微软公司内部培训程序员资料---求解非线性方程组的类

/*
 * 求解非线性方程组的类 NLEquations
 * 
 * 周长发编制
 */

using System;

namespace MSAlgorithm
{
    public abstract class NLEquations
    {
        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">变量</param>
        /// <returns>double型,函数值</returns>
        protected virtual double Func(double x)
        {
            return 0;
        }

        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">变量值数组</param>
        /// <returns>double型,函数值</returns>
        protected virtual double Func(double[] x)
        {
            return 0;
        }

        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">变量</param>
        /// <param name="y">函数值数组</param>
        protected virtual void Func(double x, double[] y)
        {
        }

        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">二元函数的变量</param>
        /// <param name="y">二元函数的变量</param>
        /// <returns>double型,函数值</returns>
        protected virtual double Func(double x, double y)
        {
            return 0;
        }

        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">二元函数的变量值数组</param>
        /// <param name="y">二元函数的变量值数组</param>
        /// <returns>double型,函数值</returns>
        protected virtual double Func(double[] x, double[] y)
        {
            return 0;
        }

        /// <summary>
        /// 虚函数：计算方程左端函数值，必须在引申类中覆盖该类函数
        /// </summary>
        /// <param name="x">已知变量值数组</param>
        /// <param name="p">已知函数值数组</param>
        protected virtual void FuncMJ(double[] x, double[] p)
        {
        }

        /// <summary>
        /// 基本构造函数
        /// </summary>
        public NLEquations()
        {
        }

        /// <summary>
        /// 求非线性方程实根的对分法
        /// 调用时，须覆盖计算方程左端函数f(x)值的虚函数
        /// double Func(double x)
        /// </summary>
        /// <param name="nNumRoots">在[xStart, xEnd]内实根个数的预估值</param>
        /// <param name="x">一维数组，长度为m。返回在区间[xStart, xEnd]内搜索到的实根，实根个数由函数值返回</param>
        /// <param name="xStart">求根区间的左端点</param>
        /// <param name="xEnd">求根区间的右端点</param>
        /// <param name="dblStep">搜索求根时采用的步长</param>
        /// <param name="eps">精度控制参数</param>
        /// <returns>int型，求得的实根的数目</returns>
        public int GetRootBisect(int nNumRoots, double[] x, double xStart, double xEnd, double dblStep, double eps)
        {
            int n, js;
            double z, y, z1, y1, z0, y0;

            //根的个数清0
            n = 0;

            //从左端点开始搜索
            z = xStart;
            y = Func(z);

            //循环求解
            while ((z <= xEnd + dblStep / 2.0) && (n != nNumRoots))
            {
                if (Math.Abs(y) < eps)
                {
                    n = n + 1;
                    x[n - 1] = z;
                    z = z + dblStep / 2.0;
                    y = Func(z);
                }
                else
                {
                    z1 = z + dblStep;
                    y1 = Func(z1);

                    if (Math.Abs(y1) < eps)
                    {
                        n = n + 1;
                        x[n - 1] = z1;
                        z = z1 + dblStep / 2.0;
                        y = Func(z);
                    }
                    else if (y * y1 > 0.0)
                    {
                        y = y1;
                        z = z1;
                    }
                    else
                    {
                        js = 0;
                        while (js == 0)
                        {
                            if (Math.Abs(z1 - z) < eps)
                            {
                                n = n + 1;
                                x[n - 1] = (z1 + z) / 2.0;
                                z = z1 + dblStep / 2.0; y = Func(z);
                                js = 1;
                            }
                            else
                            {
                                z0 = (z1 + z) / 2.0;
                                y0 = Func(z0);
                                if (Math.Abs(y0) < eps)
                                {
                                    x[n] = z0;
                                    n = n + 1;
                                    js = 1;
                                    z = z0 + dblStep / 2.0;
                                    y = Func(z);
                                }
                                else if ((y * y0) < 0.0)
                                {
                                    z1 = z0;
                                    y1 = y0;
                                }
                                else
                                {
                                    z = z0;
                                    y = y0;
                                }
                            }
                        }
                    }
                }
            }

            //返回实根的数目
            return (n);
        }

        /// <summary>
        /// 求非线性方程一个实根的牛顿法
        /// 调用时，须覆盖计算方程左端函数f(x)及其一阶导数f'(x)值的虚函数:
        /// void Func(double x, double[] y)
        /// y(0) 返回f(x)的值
        /// y(1) 返回f'(x)的值
        /// </summary>
        /// <param name="x">传入迭代初值（猜测解），返回在区间求得的一个实根</param>
        /// <param name="nMaxIt">递归次数</param>
        /// <param name="eps">精度控制参数</param>
        /// <returns>bool 型，求解是否成功</returns>
        public bool GetRootNewton(ref double x, int nMaxIt, double eps)
        {
            int l;
            double d, p, x0, x1 = 0.0;
            double[] y = new double[2];

            //条件值
            l = nMaxIt;
            x0 = x;
            Func(x0, y);

            //求解，控制精度
            d = eps + 1.0;
            while ((d >= eps) && (l != 0))
            {
                if (y[1] == 0.0)
                    return false;

                x1 = x0 - y[0] / y[1];
                Func(x1, y);

                d = Math.Abs(x1 - x0);
                p = Math.Abs(y[0]);
                if (p > d)
                    d = p;
                x0 = x1;
                l = l - 1;
            }

            x = x1;

            return true;
        }

        /// <summary>
        /// 求非线性方程一个实根的埃特金迭代法
        /// 调用时，须覆盖计算方程左端函数f(x)值的虚函数
        /// double Func(double x)
        /// </summary>
        /// <param name="x">传入迭代初值（猜测解），返回在区间求得的一个实根</param>
        /// <param name="nMaxIt">递归次数</param>
        /// <param name="eps">精度控制参数</param>
        /// <returns>bool 型，求解是否成功</returns>
        public bool GetRootAitken(ref double x, int nMaxIt, double eps)
        {
            int flag, l;
            double u, v, x0;

            //求解条件
            l = 0;
            x0 = x;
            flag = 0;

            //迭代求解
            while ((flag == 0) && (l != nMaxIt))
            {
                l = l + 1;
                u = Func(x0);
                v = Func(u);
                if (Math.Abs(u - v) < eps)
                {
                    x0 = v;
                    flag = 1;
                }
                else
                    x0 = v - (v - u) * (v - u) / (v - 2.0 * u + x0);
            }

            x = x0;

            //是否在指定的迭代次数内达到求解精度
            return (nMaxIt > l);
        }

        /// <summary>
        /// 求非线性方程一个实根的连分式解法
        /// 调用时，须覆盖计算方程左端函数f(x)值的虚函数
        /// double Func(double x)
        /// </summary>
        /// <param name="x">传入迭代初值（猜测解），返回在区间求得的一个实根</param>
        /// <param name="eps">精度控制参数</param>
        /// <returns>bool 型，求解是否成功</returns>
        public bool GetRootPq(ref double x, double eps)
        {
            int i, j, m, it = 0, l;
            double z, h, x0, q;
            double[] a = new double[10];
            double[] y = new double[10];

            //求解条件
            l = 10;
            q = 1.0e+35;
            x0 = x;
            h = 0.0;

            //连分式求解
            while (l != 0)
            {
                l = l - 1;
                j = 0;
                it = l;
                while (j <= 7)
                {
                    if (j <= 2)
                        z = x0 + 0.1 * j;
                    else
                        z = h;

                    y[j] = Func(z);
                    h = z;
                    if (j == 0)
                        a[0]