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);

後退差分法

 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);
 
          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);

中心差分法

 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);