【C#】ニュートン法を使って連立非線形方程式の解を求める - 旅行好きなソフトエンジニアの備忘録

ニュートン法を使って連立非線形方程式の解を求めます。ニュートン法を使って連立非線形方程式の解を求める場合、ヤコビアンの計算が必要となり多くの計算量が必要となってしまいます。このため、ここではニュートン法の変種を使って連立非線形方程式の解を求めるクラスを実装します。ニュートン法の変種を使った場合、計算回数は増えますが、一回当たりの計算量が減ることにより、解が求まるまでの時間が短縮される可能性が高まります。連立非線形方程式の解を求めるクラスは以下になります。このクラスではここ（【C#】偏微分を計算するクラス - 旅行好きなソフトエンジニアの備忘録）で作成したPartialDerivativeといクラスを利用しています。

public class NewtonMethodMD
{
    /// <summary>
    /// ニュートン法（の変種）により非線形連立方程式の解を求めるメソッド
    /// </summary>
    /// <param name="fs">非線形方程式群</param>
    /// <param name="initialX">解の初期値</param>
    /// <param name="iteration">計算回数</param>
    /// <param name="acceptableError">許容誤差</param>
    /// <param name="xn">計算終了後の解</param>
    /// <returns>許容誤差以内で解が求まった場合にtrue</returns>
    public static bool Compute(Func<double[], double>[] fs, double[] initialX, int iteration, double acceptableError, out double[] xn)
    {
        ValidateArguments(fs, initialX, iteration, acceptableError);

        bool found = false;
        xn = new double[initialX.Length];
        Array.Copy(initialX, xn, xn.Length);
        for (int i = 0; i < iteration; ++i)
        {
            double[] xnCopy = xn;
            double error = fs.Select(f => Math.Abs(f(xnCopy))).Sum();
            if (error <= acceptableError)
            {
                found = true;
                break;
            }

            // 一変数ずつ解を更新する
            for (int j = 0; j < fs.Length; ++j)
            {
                double fx = fs[j](xn);
                double fdx = PartialDerivative.Compute(fs[j], xn, j);
                xn[j] -= fx/fdx;
            }
        }

        return found;
    }

    private static void ValidateArguments(Func<double[], double>[] fs, double[] initialX, int iteration, double acceptableError)
    {
        if (fs == null)
        {
            throw new ArgumentNullException("fs");
        }
        if (initialX == null)
        {
            throw new ArgumentNullException("initialX");
        }
        if (fs.Length != initialX.Length)
        {
            throw new ArgumentException("fs count and initialX length must be the same");
        }
        if (fs.Length <= 1)
        {
            throw new ArgumentException("fs length must be more than 1");
        }
        if (iteration <= 0)
        {
            throw new ArgumentOutOfRangeException("iteration", "iteration must be positive");
        }
        if (acceptableError <= 0.0)
        {
            throw new ArgumentOutOfRangeException("acceptableError", "acceptableError must be positive");
        }
    }
}

次にNewtonMethodMDクラスを使って以下の連立非線形方程式の解を求めます。
$f_0(x_0,x_1)=-0.1x_0^2+x_0-0.2x_1^2-0.7=0$
$f_1(x_0,x_1)=-0.1x_0-0.2x_0x_1^2+x_1-0.7=0$

// f0(x0,x1)
Func<double[], double> f0 = x => -0.1*x[0]*x[0] + x[0] - 0.2*x[1]*x[1] - 0.7;
// f1(x0,x1)
Func<double[], double> f1 = x => -0.1*x[0] - 0.2*x[0]*x[1]*x[1] + x[1] - 0.7;
var fs = new Func<double[], double>[] { f0, f1 };
// (x0,x1)の初期値を(0.5,0.5)とする
var initialX = new double[] { 0.5, 0.5 };
// 計算回数は30
int iteration = 30;
// 許容誤差は0.001
double acceptableError = 0.001;
double[] xn;
bool found = NewtonMethodMD.Compute(fs, initialX, iteration, acceptableError, out xn);