C Sharpと数値解析 - オイラー法

前進差分法

前進差分法 (Forward Euler Method) は、精度は1次であるが、最も簡単な実装であり、現在の点での傾きを用いて次の点を予測する。

これは、${\dfrac {df(x)}{dx}}={\dfrac {f(x+h)-f(x)}{h}}$ に基づいている。

 using System;
 
 class ForwardEuler
 {
    // 微分方程式を表す関数の定義
    // 引数: (x, y), 戻り値: dy/dx
    private readonly Func<double, double, double> _differentialEquation;
 
    public ForwardEuler(Func<double, double, double> differentialEquation)
    {
       _differentialEquation = differentialEquation;
    }
 
    public void Solve(double x0, double y0, double h, double xn)
    {
       double x = x0;  // 初期x値
       double y = y0;  // 初期y値
 
       Console.WriteLine("前進差分法による計算結果 : ");
       Console.WriteLine($"{"x",12}{"y",12}");
       Console.WriteLine($"{x,12:F6}{y,12:F6}");
 
       // xがxnに達するまで計算を繰り返す
       while (x < xn)
       {
          // 前進差分公式 : y(n + 1) = yn + h * f(xn, yn)
          y = y + h * _differentialEquation(x, y);
          x = x + h;
 
          // 数値の右寄せ表示 (,12:F6)
          Console.WriteLine($"{x,12:F6}{y,12:F6}");
       }
    }
 }

 // 使用例
 
 // 微分方程式 dy/dx = x + y を定義
 Func<double, double, double> equation = (x, y) => x + y;
 
 // パラメータの設定
 double x0 = 0.0;    // 解析開始するxの初期値
 double y0 = 1.0;    // 解析開始するyの初期値
 double h  = 0.01;   // 刻み幅
 double xn = 1.0;    // xnまで計算
 
 // オイラー法 (前進差分法) による計算
 var solver = new ForwardEuler(equation);
 solver.Solve(x0, y0, h, xn);

後退差分法 (ニュートン法を使用しない方法)

ニュートン法を使用しない後退差分法 (Backward Euler Method) は、次式で示す点での傾きを使用して計算を行う陰的な方法である。

${\dfrac {df(x)}{dx}}={\dfrac {f(x)-f(x-h)}{h}}$ に基づいている。

現在の値を使用して次のステップを1回で計算する。
ニュートン法を用いない後退差分法は、計算が単純で高速となり、メモリ使用量が少ない。
ただし、ニュートン法を用いる後退差分法と比較して、精度は低く安定性も劣る。

 using System;
 
 class BackwardEuler
 {
    // 微分方程式を表す関数の定義
    // 引数 : (x, y), 戻り値 : dy/dx
    private readonly Func<double, double, double> _differentialEquation;
 
    public BackwardEuler(Func<double, double, double> differentialEquation)
    {
       _differentialEquation = differentialEquation;
    }
 
    public void Solve(double x0, double y0, double h, double xn)
    {
       double x = x0;  // 初期x値
       double y = y0;  // 初期y値
 
       Console.WriteLine("後退差分法による計算結果 : ");
       Console.WriteLine($"{"x",12}{"y",12}");
       Console.WriteLine($"{x,12:F6}{y,12:F6}");
 
       // xがxnに達するまで計算を繰り返す
       while (x < xn)
       {
          x = x + h;      // x値を更新
 
          // 後退差分公式: y(n + 1) = yn + h * f(x(n + 1), yn)
          y = y + h * _differentialEquation(x, y);
 
          // 数値の右寄せ表示 (,12:F6)
          Console.WriteLine($"{x,12:F6}{y,12:F6}");
       }
    }
 }

 // 使用例
 
 // 微分方程式 dy/dx = x + y を定義
 Func<double, double, double> equation = (x, y) => x + y;
 
 // パラメータの設定
 double x0 = 0.0;    // 解析開始するxの初期値
 double y0 = 1.0;    // 解析開始するyの初期値
 double h  = 0.01;   // 刻み幅
 double xn = 1.0;    // xnまで計算
 
 // オイラー法 (後退差分法) による計算
 var solver = new BackwardEuler(equation);
 solver.Solve(x0, y0, h, xn);

後退差分法 (ニュートン法を使用する方法)

ニュートン法を使用する後退差分法 (Backward Euler Method) は、次式で示す点での傾きを使用して計算を行う陰的な方法である。
ニュートン法による反復計算が必要で計算コストは高くなるが、精度および安定性は高い。

${\dfrac {df(x)}{dx}}={\dfrac {f(x)-f(x-h)}{h}}$ に基づいている。

ニュートン法では、収束判定のための許容誤差を設定して、最大反復回数による安全装置を実装する必要がある。
また、以下の例では、初期推定値として前進オイラー法の結果を使用している。

ただし、収束性は初期推定値に依存する場合があることに注意する。

 using System;
 
 class BackwardEulerNewton
 {
    // 微分方程式を表す関数の定義
    private readonly Func<double, double, double> _differentialEquation;
 
    // 微分方程式のy偏導関数を表す関数の定義
    private readonly Func<double, double, double> _partialDerivative;
 
    // ニュートン法の収束判定のための許容誤差
    private const double Tolerance = 1e-10;
 
    // ニュートン法の最大反復回数
    private const int MaxIterations = 100;
 
    public BackwardEulerNewton(Func<double, double, double> differentialEquation, Func<double, double, double> partialDerivative)
    {
       _differentialEquation = differentialEquation;
       _partialDerivative = partialDerivative;
    }
 
    // ニュートン法による非線形方程式の求解
    private double NewtonSolver(double x_next, double y_current, double h, double initial_guess)
    {
       double y = initial_guess;
 
       for (int i = 0; i < MaxIterations; i++)
       {
          // F(y) = y - y_current - h * f(x_next, y)
          double F = y - y_current - h * _differentialEquation(x_next, y);
 
          // F'(y) = 1 - h * df/dy(x_next, y)
          double dF = 1 - h * _partialDerivative(x_next, y);
 
          // ニュートン法による更新
          double y_new = y - F / dF;
 
          // 収束判定
          if (Math.Abs(y_new - y) < Tolerance)
          {
             return y_new;
          }
 
          y = y_new;
       }
 
       // 収束しない場合は警告を出して最後の値を返す
       Console.WriteLine("Warning: Newton method did not converge");
 
       return y;
    }
 
    public void Solve(double x0, double y0, double h, double xn)
    {
       double x = x0;  // 初期x値
       double y = y0;  // 初期y値
 
       Console.WriteLine("後退差分法 (ニュートン法使用) による計算結果 :");
       Console.WriteLine($"{"x",12}{"y",12}");
       Console.WriteLine($"{x,12:F6}{y,12:F6}");
 
       // xがxnに達するまで計算を繰り返す
       while (x < xn)
       {
          double x_next = x + h;
 
          // 初期推定値として前進オイラー法の結果を使用
          double initial_guess = y + h * _differentialEquation(x, y);
 
          // ニュートン法で次のステップの値を求解
          y = NewtonSolver(x_next, y, h, initial_guess);
          x = x_next;
 
          // 数値の右寄せ表示 (,12:F6)           
          Console.WriteLine($"{x,12:F6}{y,12:F6}");
       }
    }
 }

 // 微分方程式 dy/dx = x + y を定義
 Func<double, double, double> equation = (x, y) => x + y;
 
 // 微分方程式のy偏導関数 ∂f/∂y = 1 を定義
 Func<double, double, double> derivative = (x, y) => 1.0;
 
 // パラメータの設定
 // パラメータの設定
 double x0 = 0.0;    // 解析開始するxの初期値
 double y0 = 1.0;    // 解析開始するyの初期値
 double h  = 0.01;   // 刻み幅
 double xn = 1.0;    // xnまで計算
 
 // オイラー法 (後退差分法) による計算
 var solver = new BackwardEulerNewton(equation, derivative);
 solver.Solve(x0, y0, h, xn);

中心差分法

 using System;
 
 class CentralEuler
 {
    // 微分方程式を表す関数の定義
    // 引数 : (x, y), 戻り値 : dy/dx
    private readonly Func<double, double, double> _differentialEquation;
 
    public CentralEuler(Func<double, double, double> differentialEquation)
    {
       _differentialEquation = differentialEquation;
    }
 
    public void Solve(double x0, double y0, double h, double xn)
    {
       double x = x0;          // 初期x値
       double y = y0;          // 初期y値
       double y_prev = y0;     // 前のステップのy値
       double y_next;          // 次のステップのy値
 
       Console.WriteLine("中心差分法による計算結果 :");
       Console.WriteLine($"{"x",12}{"y",12}");
       Console.WriteLine($"{x,12:F6}{y,12:F6}");
 
       // 最初のステップは前進差分で計算
       y = y + h * _differentialEquation(x, y);
       x = x + h;
       Console.WriteLine($"{x,12:F6}{y,12:F6}");
 
       // xがxnに達するまで計算を繰り返す
       while (x < xn)
       {
          // 中心差分公式 : y(n + 1) = y(n - 1) + 2h * f(xn, yn)
          y_next = y_prev + 2 * h * _differentialEquation(x, y);
          x = x + h;
 
          y_prev = y;     // 現在の値を前の値として保存
          y = y_next;     // 次の値を現在の値に更新
 
          Console.WriteLine($"{x,12:F6}{y,12:F6}");
       }
    }
 }

 // 使用例
 
 // 微分方程式 dy/dx = x + y を定義
 Func<double, double, double> equation = (x, y) => x + y;
 
 // パラメータの設定
 double x0 = 0.0;    // 解析開始するxの初期値
 double y0 = 1.0;    // 解析開始するyの初期値
 double h  = 0.01;   // 刻み幅
 double xn = 1.0;    // xnまで計算
 
 // オイラー法 (中心差分法) による計算
 var solver = new CentralEuler(equation);
 solver.Solve(x0, y0, h, xn);