C#---牛顿迭代法求解非线性方程组

NewtonMethod.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

/// <summary>
/// 非线性方程组
/// by Ted
/// 2013.6
/// </summary>
/// 
namespace NewtonIteration
{
    class NewtonMethod
    {
        private NonlinearEquations _Equations; //非线性方程组
        private NonlinearEquation[,] _JacobiEquations; //雅克比方程(_EquationNum行,_UnknowNum列)
        private int _UnknowNum; //未知数的个数
        private int _EquationNum; //方程的个数
        private int _MaxIterativeTime;// 最大迭代次数
        private double _Precision; // 精度
        private double[] _StartValue; // 初值矩阵
        private double[] _EndValue; // 计算结果

        public double[] EndValue
        {
            get { return _EndValue; }
        }

        /// <summary>
        /// 默认构造函数
        /// </summary>
        public NewtonMethod()
        {
        }

        /// <summary>
        /// 通过非线性方程组构造牛顿迭代
        /// </summary>
        /// <param name="e">原始方程组</param>
        /// <param name="u">未知数个数</param>
        /// <param name="n">方程的个数</param>
        /// <param name="m">最大迭代次数</param>
        /// <param name="p">精度</param>
        /// <param name="s">初值</param>
        public NewtonMethod(NonlinearEquations e, int u, int n, int m, double p, double[] s)
        {
            _UnknowNum = u;
            _EquationNum = n;
            _MaxIterativeTime = m;
            _Precision = p;
            _StartValue = new double[u];
            for (int i = 0; i < u; i++)
            {
                _StartValue[i] = s[i];
            }
            _Equations = new NonlinearEquations();
            _Equations.Clear();
            for (int i = 0; i < n; i++)
            {
                _Equations.Add(e[i]);
            }
            _JacobiEquations = new NonlinearEquation[_EquationNum, _UnknowNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    _JacobiEquations[i, j] = new NonlinearEquation(_UnknowNum);
                    _JacobiEquations[i, j] = PartialDifferential(_Equations[i], j);
                }
            }
            _EndValue = new double[u];
        }

        /// <summary>
        /// 非线性方程e对未知数n求偏微分
        /// </summary>
        /// <param name="e">原始方程</param>
        /// <param name="n">未知数n从0开始计数</param>
        /// <returns></returns>
        private NonlinearEquation PartialDifferential(NonlinearEquation e, int n)
        {
            NonlinearEquation result = new NonlinearEquation(e._UnknowNum);
            for (int i = 0; i < e._UnknowNum; i++)
            {
                if (i != n)
                {
                    result._Coefficient[i] = 0;
                    result._Power[i] = 0;
                }
            }

            if (e._Power[n] > 1)
            {
                result._Power[n] = e._Power[n] - 1; // 降次
                result._Coefficient[n] = e._Coefficient[n] * e._Power[n]; // 系数乘以次数
                result._Constant = 0; // 常数项归零
            }
            else
            {
                //降次后变成常数项
                result._Power[n] = 1; // 次数归1
                result._Coefficient[n] = 0; // 系数归0
                result._Constant = -e._Coefficient[n]; // 移项
            }
            return result;
        }

        /// <summary>
        /// 根据初值构造雅克比矩阵
        /// </summary>
        /// <param name="x">初值矩阵</param>
        private double[,] JacobiMatrix(double[] x)
        {
            double[,] JacobiM = new double[_EquationNum, _UnknowNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    JacobiM[i, j] = _JacobiEquations[i, j].EquationValue(x);
                }
            }
            return JacobiM;
        }

        /// <summary>
        /// 返回方程组e的系数矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[,] CoefficientMatrix(NonlinearEquations e)
        {
            double[,] c = new double[_EquationNum, _UnknowNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    c[i, j] = e[i]._Coefficient[j];
                }
            }
            return c;
        }

        /// <summary>
        /// 返回方程组e的系数矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[,] CoefficientMatrix(NonlinearEquation[] e)
        {
            double[,] c = new double[_EquationNum, _UnknowNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    c[i, j] = e[i]._Coefficient[j];
                }
            }
            return c;
        }

        /// <summary>
        /// 返回方程组e的常数项矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[] ConstantMatrix(NonlinearEquations e)
        {
            double[] y = new double[_EquationNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                y[i] = e[i]._Constant;
            }
            return y;
        }

        /// <summary>
        /// 返回方程组e的常数项矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[] ConstantMatrix(NonlinearEquation[] e)
        {
            double[] y = new double[_EquationNum];
            for (int i = 0; i < _EquationNum; i++)
            {
                y[i] = e[i]._Constant;
            }
            return y;
        }

        /// <summary>
        /// 返回方程组e的增广矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[,] AugmentationMatrix(NonlinearEquations e)
        {
            double[,] a = new double[_EquationNum, _UnknowNum + 1];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    a[i, j] = e[i]._Coefficient[j];
                }
                a[i, _UnknowNum] = e[i]._Constant;
            }
            return a;
        }

        /// <summary>
        /// 返回方程组e的增广矩阵
        /// </summary>
        /// <param name="e">方程组</param>
        private double[,] AugmentationMatrix(NonlinearEquation[] e)
        {
            double[,] a = new double[_EquationNum, _UnknowNum + 1];
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    a[i, j] = e[i]._Coefficient[j];
                }
                a[i, _UnknowNum] = e[i]._Constant;
            }
            return a;
        }

        /// <summary>
        /// 判断一个浮点数是否等于0
        /// </summary>
        private bool IsDoubleZero(double v)
        {
            if (v < 0.0001 && v > -0.0001)
            {
                return true;
            }
            return false;
        }

        /// <summary>
        /// 判断方程组是否为线性
        /// </summary>
        private bool IsEquationsLinear()
        {
            for (int i = 0; i < _EquationNum; i++)
            {
                for (int j = 0; j < _UnknowNum; j++)
                {
                    if (_Equations[i]._Power[j] >= 2) { return false; }
                }
            }
            return true;
        }

        /// <summary>
        /// 根据系数矩阵c和常数项y求解线性方程组
        /// </summary>
        /// <param name="c">系数矩阵</param>
        /// <param name="y">常数项矩阵</param>
        /// <returns >value="0" 线性方程组无解</returns>
        /// <returns >value="1" 线性方程组有唯一解</returns>
        /// <returns >value="2" 线性方程组有无穷解</returns>
        private int LinearEquations(double[,] c, double[] y, ref double[] x)
        {
            // 行_EquationNum
            // 列_UnknowNum
            x = new double[_UnknowNum];
            int ResultForm = 1;

            #region 奇次线性方程组
            bool ZeroY = true;
            for (int i = 0; i < _EquationNum; i++)
            {
                if (!IsDoubleZero(y[i])) { ZeroY = false; break; }
            }
            #endregion

            #region 交换列主
            for (int i = 0; i < _UnknowNum; i++)
            {
                for (int j = i + 1; j < _EquationNum; j++)
                {
                    double MainValue = Math.Abs(c[i, i]); //列主的绝对值大小
                    int MainLine = i; //列主的行号
                    if (MainValue < Math.Abs(c[j, i]))
                    {
                        MainLine = j;
                        MainValue = Math.Abs(c[j, i]);
                        double temp;
                        for (int k = 0; k < _UnknowNum; k++)
                        {
                            temp = c[i, k]; c[i, k] = c[MainLine, k]; c[MainLine, k] = temp;
                        }
                        temp = y[i]; y[i] = y[MainLine]; y[MainLine] = temp;
                    }
                }
            }
            #endregion

            #region 消元
            for (int i = 0; i < _UnknowNum; i++)
            {
                if (!IsDoubleZero(c[i, i]))
                {
                    for (int j = i + 1; j < _EquationNum; j++)
                    {
                        double ratio = c[j, i] / c[i, i];
                        for (int k = 0; k < i + 1; k++)
                        {
                            c[j, k] = 0;
                        }
                        for (int k = i + 1; k < _UnknowNum; k++)
                        {
                            c[j, k] = c[j, k] - ratio * c[i, k];
                        }
                        y[j] = y[j] - ratio * y[i];
                    }
                }
            }
            #endregion

            #region 判断解的性质
            int FirstZeroColumn = 0; // 第一次系数全为0的行号
            for (int i = 0; i < _EquationNum; i++)
            {
                bool IsColumnZero = true;
                for (int j = 0; j < _UnknowNum; j++)
                {
                    // 这一行的元素全部为0
                    if (!IsDoubleZero(c[i, j])) { IsColumnZero = false; break; }
                }
                if (IsColumnZero) { break; }
                FirstZeroColumn++;
            }
            if (FirstZeroColumn < _EquationNum)
            {
                // 有系数全为0的行
                if (!IsDoubleZero(y[FirstZeroColumn])) { return 0; } //无解
                double CoefficientMatrixRank = FirstZeroColumn; //系数矩阵的秩
                if (CoefficientMatrixRank == _UnknowNum) { ResultForm = 1; } // 唯一解
                else if (CoefficientMatrixRank < _UnknowNum) { ResultForm = 2; } // 无穷解
                else { ResultForm = 2; }
            }
            else { ResultForm = 1; }
            #endregion

            #region 回代
            if (!ZeroY)
            {
                //回代
                for (int i = 0; i < _UnknowNum; i++)
                {
                    for (int j = _UnknowNum - 1; j >= 0; j--)
                    {
                        double Value = y[j];
                        for (int k = j + 1; k < _UnknowNum; k++)
                        {
                            Value = Value - c[j, k] * x[k];
                        }
                        x[j] = Value / c[j, j];
                    }
                }
            }
            else
            {
                // 齐次线性方程组,有一组零解
                for (int i = 0; i < _UnknowNum; i++) { x[i] = 0; }
            }
            #endregion

            return ResultForm;
        }

        /// <summary>
        /// 使用牛顿迭代法求非线性方程组
        /// </summary>
        /// <returns >value="0" 线性方程组无解</returns>
        /// <returns >value="1" 线性方程组有唯一解</returns>
        /// <returns >value="2" 线性方程组有无穷解</returns>
        /// <returns >value="3" 非线性方程组求解失败</returns>
        /// <returns >value="4" 非线性方程组求解成功</returns>
        /// <returns >value="5" 非线性方程组达到最大迭代次数</returns>
        public int NewtonValue()
        {
            if (IsEquationsLinear())
            {
                //是线性方程组
                double[,] c = new double[_EquationNum, _UnknowNum];
                double[] y = new double[_EquationNum];
                c = CoefficientMatrix(_Equations);
                y = ConstantMatrix(_Equations);
                return LinearEquations(c, y, ref _EndValue);
            }
            else
            {
                //非线性方程组,牛顿迭代法
                bool IsCirculate = false;
                int IterativeTime = 0;
                double[] x = new double[_UnknowNum];
                for (int j = 0; j < _UnknowNum; j++) { x[j] = _StartValue[j]; }
                do
                {
                    IsCirculate = false;
                    double[,] c = JacobiMatrix(x);
                    double[] y = new double[_EquationNum];
                    for (int j = 0; j < _EquationNum; j++) { y[j] = -_Equations[j].EquationValue(x); }
                    double[] detx = new double[_UnknowNum];
                    if (LinearEquations(c, y, ref detx) == 0) return 3;
                    for (int j = 0; j < _UnknowNum; j++) { x[j] += detx[j]; }
                    // 当差值小于初值时则迭代终止
                    for (int j = 0; j < _UnknowNum; j++)
                    {
                        if (Math.Abs(detx[j]) >= _Precision)
                        {
                            IsCirculate = true;
                            break;
                        }
                    }
                    IterativeTime++;
                    if (IterativeTime >= _MaxIterativeTime) { break; } // 达到最大迭代次数
                } while (IsCirculate);
                for (int j = 0; j < _UnknowNum; j++) { _EndValue[j] = x[j]; }
                if (!IsCirculate) { return 4; }
                if (IterativeTime >= _MaxIterativeTime) { return 5; }
                return 3;
            }
        }
    }
}

NonlinearEquations.cs

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;

/// <summary>
/// 非线性方程组
/// by Ted
/// 2013.6
/// </summary>

namespace NewtonIteration
{
    class NonlinearEquations : CollectionBase
    {
        /// <summary>
        /// 默认构造函数
        /// </summary>
        public NonlinearEquations()
        {
        }

        /// <summary>
        /// 构造函数
        /// </summary>
        /// <param name="value"></param>
        public NonlinearEquations(NonlinearEquations Equations)
        {
            AddRange(Equations);
        }

        /// <summary>
        /// 构造函数
        /// </summary>
        /// <param name="value"></param>
        public NonlinearEquations(NonlinearEquation[] Equations)
        {
            AddRange(Equations);
        }

        /// <summary>
        /// 索引符
        /// </summary>
        /// <param name="index"></param>
        public NonlinearEquation this[int index]
        {
            get { return (NonlinearEquation)List[index]; }
            set { List[index] = value; }
        }

        /// <summary>
        /// 添加一个非线性方程
        /// </summary>
        /// <param name="newEquation"></param>
        public void Add(NonlinearEquation newEquation)
        {
            List.Add(newEquation);
        }

        /// <summary>
        /// 移除一个非线性方程
        /// </summary>
        /// <param name="oldEquation"></param>
        public void Remove(NonlinearEquation oldEquation)
        {
            List.Remove(oldEquation);
        }

        /// <summary>
        /// 根据值确定该方程在方程组中的索引
        /// </summary>
        /// <param name="currentEquation"></param>
        public int IndexOf(NonlinearEquation currentEquation)
        {
            return List.IndexOf(currentEquation);
        }

        /// <summary>
        /// 查看该方程是否在方程组中
        /// </summary>
        /// <param name="currentEquation"></param>
        public bool Contains(NonlinearEquation currentEquation)
        {
            return List.Contains(currentEquation);
        }

        /// <summary>
        /// 插入一个新的非线性方程
        /// </summary>
        /// <param name="currentEquation"></param>
        public void Insert(int index, NonlinearEquation currentEquation)
        {
            List.Insert(index, currentEquation);
        }

        /// <summary>
        /// 添加非线性方程数组到当前集合
        /// </summary>
        /// <param name="Equations"></param>
        public void AddRange(NonlinearEquation[] Equations)
        {
            for (int i = 0; i < Equations.Length; i++)
            {
                this.Add(Equations[i]);
            }
        }

        /// <summary>
        /// 添加非线性方程集合到当前集合
        /// </summary>
        /// <param name="Equations"></param>
        public void AddRange(NonlinearEquations Equations)
        {
            for (int i = 0; i < Equations.Capacity; i++)
            {
                Add((NonlinearEquation)Equations.List[i]);
            }
        }
    }
}

NonlinearEquation.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace NewtonIteration
{
    class NonlinearEquation
    {
        public int _UnknowNum; // 未知数的个数
        public double[] _Coefficient; // 未知数的系数
        public int[] _Power; // 未知数的次数(权),次数为0表示常数
        public double _Constant; //常数项

        /// <summary>
        /// 默认构造函数
        /// </summary>
        public NonlinearEquation()
        {
            _UnknowNum = 0;
        }

        /// <summary>
        /// 构造一个带有n个未知数的非线性方程
        /// </summary>
        /// <param name="n">未知数的个数</param>
        public NonlinearEquation(int n)
        {
            _UnknowNum = n;
            _Coefficient = new double[n];
            _Power = new int[n];
        }

        /// <summary>
        /// 默认
        /// </summary>
        public void InitEquation()
        {
            for (int i = 0; i < _UnknowNum; i++)
            {
                _Coefficient[i] = 0;
                _Power[i] = 1;
            }
            _Constant = 0;
        }

        /// <summary>
        /// 把方程转化成字符串输出
        /// </summary>
        public string EquationToString()
        {
            if (_UnknowNum <= 0) return null;
            string equation = "";

            bool IsCoefficientAllZero = true;
            for (int i = 0; i < _UnknowNum; i++)
            {
                if (!IsDoubleZero(_Coefficient[i])) { IsCoefficientAllZero = false; break; }
            }
            if (IsCoefficientAllZero)
            {
                // 系数全为0
                equation += "0";
            }
            else
            {
                bool IsAddFirst = true;
                for (int i = 0; i < _UnknowNum; i++)
                {
                    if (_Coefficient[i] < 0)
                    {
                        if (i > 0 && !IsAddFirst)
                            equation += " + ";
                        // 系数小于0,系数加括号
                        equation += "(";
                        equation += _Coefficient[i].ToString();
                        equation += ")*X";
                        equation += (i + 1).ToString();
                        if (_Power[i] > 1)
                        {
                            equation += "^";
                            equation += _Power[i].ToString();
                        }
                        if (IsAddFirst) IsAddFirst = false;
                    }
                    else if (_Coefficient[i] > 0)
                    {
                        if (i > 0 && !IsAddFirst)
                            equation += " + ";
                        // 系数大于0,系数不加括号
                        equation += _Coefficient[i].ToString();
                        equation += "*X";
                        equation += (i + 1).ToString();
                        if (_Power[i] > 1)
                        {
                            equation += "^";
                            equation += _Power[i].ToString();
                        }
                        if (IsAddFirst) IsAddFirst = false;
                    }
                    // 系数等于0,不输出
                }
            }
            equation += " = ";
            equation += _Constant.ToString();
            return equation;
        }

        /// <summary>
        /// 设置方程未知数的值并输出
        /// </summary>
        public double EquationValue(double[] x)
        {
            double y = 0;
            for (int i = 0; i < _UnknowNum; i++)
            {
                double result = 0;
                if (_Coefficient[i] != 0)
                {
                    result = 1;
                    for (int j = 0; j < _Power[i]; j++)
                    {
                        // 次数
                        result *= x[i];
                    }
                    result *= _Coefficient[i]; // 系数
                }
                y += result;
            }
            y -= _Constant; // 常数项移项
            return y;
        }

        /// <summary>
        /// 判断一个浮点数是否等于0
        /// </summary>
        private bool IsDoubleZero(double v)
        {
            if (v < 0.0001 && v > -0.0001)
            {
                return true;
            }
            return false;
        }
    }
}

Program.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace NewtonIteration
{
    class Program
    {
        static void Main(string[] args)
        {
            int _UnknowNum = 2;                  // 未知数的个数
            int _EquationNum = 2;                // 方程的个数
            int _MaxIterativeTime = 100;          // 最大迭代次数
            double _Precision = 0.001;           // 精度
            double[] _StartValue = {1, 1};       // 初值矩阵
            double[] _EndValue;                  // 计算结果
            NonlinearEquations _Equations;       // 非线性方程组
            NewtonMethod _Newton;                // 牛顿迭代式

            // 赋值
            _Equations = new NonlinearEquations();
            NonlinearEquation[] newEquations = new NonlinearEquation[_EquationNum];
            double [,]coefficient = new double[2,2]{{1, 2},{2, 1}};
            int [,]power = new int[2, 2] { { 1, 1 }, { 2, 2 } };
            double []constant = new double[]{3, 5};

            for (int i = 0; i < _EquationNum; i++)
            {
                newEquations[i] = new NonlinearEquation(_UnknowNum);
                for (int j = 0; j < _UnknowNum; j++)
                {
                    newEquations[i]._Coefficient[j] = coefficient[i,j];
                    newEquations[i]._Power[j] = power[i,j];
                }
                newEquations[i]._Constant = constant[i];
            }
            _Equations.Clear();
            _Equations.AddRange(newEquations);

            // 调用牛顿迭代法
            if (_UnknowNum > 0 || _EquationNum > 0 || _MaxIterativeTime > 0 || _Precision > 0)
            {
                _Newton = new NewtonMethod(_Equations, _UnknowNum, _EquationNum, _MaxIterativeTime, _Precision, _StartValue);
                int CalculateResult = _Newton.NewtonValue();
                _EndValue = _Newton.EndValue;
                string mesg = "";
                // 将结果显示在屏幕上
                switch (CalculateResult)
                {
                    case 0:
                        mesg = "线性方程组求解失败！线性方程组无解。";
                        break;
                    case 1:
                        mesg = "线性方程组求解成功！有且仅有一组解。";
                        break;
                    case 2:
                        mesg = "线性方程组求解失败！有无穷多解。";
                        break;
                    case 3:
                        mesg = "非线性方程组求解失败！";
                        break;
                    case 4:
                        mesg = "非线性方程组求解成功！";
                        break;
                    case 5:
                        mesg = "达到最大迭代次数！";
                        break;
                }
                Console.WriteLine(mesg);

                if (CalculateResult == 1 || CalculateResult == 4 || CalculateResult == 5)
                {
                    for (int i = 1; i < _UnknowNum + 1; i++)
                    {
                        string x = "x" + i.ToString() + "=" + _EndValue[i - 1].ToString();
                        Console.WriteLine(x);
                    }
                }
            }
            else
            {
                Console.WriteLine("参数不全！");
                return;
            }

            Console.ReadLine();
        }
    }
}