「C Sharpと数値解析 - オイラー法」の版間の差分

📢 Webサイト閉鎖と移転のお知らせ

このWebサイトは2026年9月に閉鎖いたします。

移転先はhttps://mochiu.net/index.php?title=C Sharpと数値解析 - オイラー法です。

新しい記事は移転先で追加しております。（旧サイトでは記事を追加しておりません）

← 古い編集

ビジュアルウィキテキスト

@@ 1行目: / 1行目: @@
 == 概要 ==
+オイラー法は、常微分方程式の数値解法における最も基本的な手法の1つである。<br>
+<br>
+オイラー法の考え方は、微分方程式を差分方程式で近似することである。<br>
+与えられた微分方程式 <math>\dfrac{dy}{dx} = f(x, y)</math> において、各点での関数の傾き <math>f(x, y)</math> を用いて、次の点での関数値を予測する。<br>
+<br>
+<math>y(n + 1) = y(n) + h \times f(x(n), y(n))</math><br>
+ここで、hは刻み幅、nは現在のステップ数を表す。<br>
+<br>
+オイラー法の特徴として、実装が簡単で直感的に理解しやすいというメリットがある。<br>
+<br>
+しかし、精度の面では課題があり、刻み幅hに比例する大きさの局所打ち切り誤差が生じる。<br>
+また、各ステップでの誤差が蓄積されていくため、長い区間での計算には向いていない。<br>
+<br>
+より高精度な計算が必要な場合は、ルンゲ・クッタ法等の改良された方法を使用することが推奨される。<br>
 <br><br>
@@ 7行目: / 20行目: @@
 <br>
 これは、<math>\dfrac{df(x)}{dx} = \dfrac{f(x + h) - f(x)}{h}</math> に基づいている。<br>
+<br>
+<math>\dfrac{dy}{dx_i} = f(x_i, y_i)</math> より、<br>
+<math>f(x_i, y_i) \approx \dfrac{y_{i+1} - y_i}{h} \qquad (1)</math><br>
+<br>
+上記の(1)式を変形して、<math>y_{i+1} = y_i + h f(x_i, y_i)</math> となる。<br>
 <br>
   <syntaxhighlight lang="c#">
@@ 60行目: / 78行目: @@
   var solver = new ForwardEuler(equation);
   solver.Solve(x0, y0, h, xn);
+ </syntaxhighlight>
+<br><br>
+== 前進差分法 (誤差管理機能付き) ==
+誤差管理機能付きの前進差分法は、適応的な刻み幅の制御を行うオイラー法である。<br>
+<br>
+誤差を監視しながら刻み幅を動的に調整する。<br>
+許容誤差、最小・最大刻み幅を指定することができる。<br>
+<br>
+ <syntaxhighlight lang="c#">
+ 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);
+    }
+ }
+ </syntaxhighlight>
+<br>
+ <syntaxhighlight lang="c#">
+ // 誤差管理機能付きの前進差分法による解法
+ // テスト用の微分方程式: 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}");
+ }
   </syntaxhighlight>
 <br><br>
@@ 68行目: / 223行目: @@
 <math>\dfrac{df(x)}{dx} = \dfrac{f(x) - f(x - h)}{h}</math> に基づいている。<br>
 <br>
-現在の値を使用して次のステップを1回で計算する。<br>
+<math>\dfrac{dy}{dx_i} = f(x_i, y_i)</math> より、<br>
+<math>f(x_i, y_i) \approx \dfrac{y_i - y_{i-1}}{h} \qquad (1)</math><br>
+<br>
+上記の(1)式を変形して、<math>y_{i+1} = y_i + h f(x_{i+1}, y_{i+1})</math> となる。<br>
+<br>
+後退差分法の式では、左辺と右辺の両方に <math>y_{i+1}</math> が含まれているため、今求まっている <math>f(x, y)</math> から単純な計算でを求めることはできない。<br>
+ただし、<math>f(x, y)</math> に具体的な関数を代入することにより、上記の方程式を代数的あるいは数値的に解くことで <math>y_{i+1}</math> を求めることできる。<br>
+<br>
+例えば、<math>f(x, y) = x + y</math> の場合の後退差分法の式は <math>y_{i+1} = \frac{y_i + h (x_{i} + h)}{(1 - h)}</math> となるため、<br>
+現在の <math>y_i, x_i</math> のみから計算可能である。<br>
+<br>
+後退差分法では、現在の値を使用して次のステップを1回で計算する。<br>
 ニュートン法を用いない後退差分法は、計算が単純で高速となり、メモリ使用量が少ない。<br>
 ただし、ニュートン法を用いる後退差分法と比較して、精度は低く安定性も劣る。<br>