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

概要

オイラー法は、常微分方程式の数値解法における最も基本的な手法の1つである。

オイラー法の考え方は、微分方程式を差分方程式で近似することである。
与えられた微分方程式 ${\dfrac {dy}{dx}}=f(x,y)$ において、各点での関数の傾き $f(x,y)$ を用いて、次の点での関数値を予測する。

$y(n+1)=y(n)+h\times f(x(n),y(n))$
ここで、hは刻み幅、nは現在のステップ数を表す。

オイラー法の特徴として、実装が簡単で直感的に理解しやすいというメリットがある。

しかし、精度の面では課題があり、刻み幅hに比例する大きさの局所打ち切り誤差が生じる。
また、各ステップでの誤差が蓄積されていくため、長い区間での計算には向いていない。

より高精度な計算が必要な場合は、ルンゲ・クッタ法等の改良された方法を使用することが推奨される。

前進差分法

前進差分法 (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);

前進差分法 (誤差管理機能付き)

誤差管理機能付きの前進差分法は、適応的な刻み幅の制御を行うオイラー法である。

誤差を監視しながら刻み幅を動的に調整する。
許容誤差、最小・最大刻み幅を指定することができる。

 using System;
 using System.Collections.Generic;
 
 /// <summary>
 /// 適応的な刻み幅の制御を備えたオイラー法による微分方程式ソルバ
 /// </summary>
 public class AdaptiveEuler
 {
    private readonly Func<double, double, double> _function;
    private readonly double _initialStepSize;
    private readonly double _tolerance;
    private readonly double _minStepSize;
    private readonly double _maxStepSize;
 
    /// <summary>
    /// コンストラクタ
    /// </summary>
    /// <param name="function">微分方程式の右辺を表す関数 f(t, y)</param>
    /// <param name="initialStepSize">初期刻み幅</param>
    /// <param name="tolerance">許容誤差</param>
    /// <param name="minStepSize">最小刻み幅</param>
    /// <param name="maxStepSize">最大刻み幅</param>
    public AdaptiveEuler(Func<double, double, double> function, double initialStepSize = 0.01,
                         double tolerance = 1e-6, double minStepSize = 1e-8, double maxStepSize = 0.1)
    {
       _function = function;
       _initialStepSize = initialStepSize;
       _tolerance = tolerance;
       _minStepSize = minStepSize;
       _maxStepSize = maxStepSize;
    }
 
    /// <summary>
    /// 指定された区間で微分方程式を解く
    /// </summary>
    /// <param name="t0">開始時刻</param>
    /// <param name="y0">初期値</param>
    /// <param name="tEnd">終了時刻</param>
    /// <returns>計算結果の時刻と解のペアのリスト</returns>
    public List<(double Time, double Value, double StepSize, double Error)> Solve(double t0, double y0, double tEnd)
    {
       var result = new List<(double Time, double Value, double StepSize, double Error)>();
       var t = t0;
       var y = y0;
       var h = _initialStepSize;
 
       // 初期値を結果リストに追加
       result.Add((t, y, h, 0.0));
 
       while (t < tEnd)
       {
          // 現在の刻み幅が終了時刻を超えないように調整
          if (t + h > tEnd)
          {
             h = tEnd - t;
          }
 
          // 2つの異なる刻み幅で解を計算して、誤差を推定
          var y1 = SingleStep(t, y, h);                    // 1ステップで計算
          var y2_1 = SingleStep(t, y, h/2);               // 半分のステップで1回目
          var y2 = SingleStep(t + h/2, y2_1, h/2);        // 半分のステップで2回目
 
          // 2つの解の差から誤差を推定
          var error = Math.Abs(y2 - y1);
 
          // 誤差が許容範囲内かどうかを確認
          if (error <= _tolerance || h <= _minStepSize)
          {
             // 解を更新
             t += h;
             y = y2;  // より精度の高い解を採用
 
             // 結果を保存
             result.Add((t, y, h, error));
 
             // 次の刻み幅を調整 (誤差が小さすぎる場合は大きく、大きすぎる場合は小さく)
             if (error < _tolerance / 10 && h < _maxStepSize)
             {
                h = Math.Min(h * 2, _maxStepSize);
             }
          }
          else
          {
             // 誤差が大きすぎる場合は刻み幅を半分にして再試行
             h = Math.Max(h / 2, _minStepSize);
          }
       }
 
       return result;
    }
 
    /// <summary>
    /// オイラー法による1ステップの計算
    /// </summary>
    /// <param name="t">現在の時刻</param>
    /// <param name="y">現在の値</param>
    /// <param name="h">刻み幅</param>
    /// <returns>次のステップでの値</returns>
    private double SingleStep(double t, double y, double h)
    {
       return y + h * _function(t, y);
    }
 }

 // 誤差管理機能付きの前進差分法による解法
 
 // テスト用の微分方程式: dy/dt = -y
 // 解析解 : y = y0 * e^(-t)
 Func<double, double, double> f = (t, y) => -y;
 
 var solver = new AdaptiveEuler(
    function        : f
    initialStepSize : 0.1,
    tolerance       : 1e-6,
    minStepSize     : 1e-8,
    maxStepSize     : 0.5
 );
 
 var result = solver.Solve(t0: 0, y0: 1, tEnd: 5);
 
 // 結果の出力
 foreach (var (time, value, stepSize, error) in result)
 {
    Console.WriteLine($"t: {time:F6}, y: {value:F6}, h: {stepSize:E3}, error: {error:E3}");
 }

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

ニュートン法を使用しない後退差分法 (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);

中心差分法

中心差分法 (Central Euler Method) は、前後の点を使用して中心での傾きを計算する。

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

2次の精度を持ち、前進差分法より精度が高いが、初期ステップの処理が必要となる。

 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;  // 次の値を現在の値に更新
 
          // 数値の右寄せ表示 (,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 CentralEuler(equation);
 solver.Solve(x0, y0, h, xn);