MochiuWiki : SUSE, EC, PCB

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

📢 Webサイト閉鎖と移転のお知らせ
このWebサイトは2026年9月に閉鎖いたします。
新しい記事は移転先で追加しております。（旧サイトでは記事を追加しておりません）

← 古い編集
ビジュアルウィキテキスト
r
 
(同じ利用者による、間の5版が非表示)
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>
<br><br>


7行目: 20行目:
<br>
<br>
これは、<math>\dfrac{df(x)}{dx} = \dfrac{f(x + h) - f(x)}{h}</math> に基づいている。<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>
<br>
  <syntaxhighlight lang="c#">
  <syntaxhighlight lang="c#">
60行目: 78行目:
  var solver = new ForwardEuler(equation);
  var solver = new ForwardEuler(equation);
  solver.Solve(x0, y0, h, xn);
  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>
  </syntaxhighlight>
<br><br>
<br><br>
68行目: 223行目:
<math>\dfrac{df(x)}{dx} = \dfrac{f(x) - f(x - h)}{h}</math> に基づいている。<br>
<math>\dfrac{df(x)}{dx} = \dfrac{f(x) - f(x - h)}{h}</math> に基づいている。<br>
<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>
ただし、ニュートン法を用いる後退差分法と比較して、精度は低く安定性も劣る。<br>
ただし、ニュートン法を用いる後退差分法と比較して、精度は低く安定性も劣る。<br>
242行目: 408行目:


== 中心差分法 ==
== 中心差分法 ==
中心差分法 (Central Euler Method) は、前後の点を使用して中心での傾きを計算する。<br>
<br>
<math>\dfrac{df(x)}{dx} = \dfrac{f(x + h) - f(x - h)}{2h}</math> に基づいている。<br>
<br>
2次の精度を持ち、前進差分法より精度が高いが、初期ステップの処理が必要となる。<br>
<br>
  <syntaxhighlight lang="c#">
  <syntaxhighlight lang="c#">
  using System;
  using System;
279行目: 451行目:
           x = x + h;
           x = x + h;
   
   
           y_prev = y;     // 現在の値を前の値として保存
           y_prev = y; // 現在の値を前の値として保存
           y = y_next;     // 次の値を現在の値に更新
           y = y_next; // 次の値を現在の値に更新
   
   
          // 数値の右寄せ表示 (,12:F6)
           Console.WriteLine($"{x,12:F6}{y,12:F6}");
           Console.WriteLine($"{x,12:F6}{y,12:F6}");
       }
       }
「http://mochiuwiki.e2.valueserver.jp/index.php?title=C_Sharpと数値解析_-_オイラー法」から取得