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

ナビゲーションに移動 検索に移動

C Sharpと数値解析 - オイラー法 (ソースを閲覧)

2025年1月22日 (水) 14:16時点における版

4,761 バイト追加 、 2025年1月22日 (水)
→‎前進差分法
60行目: 60行目:
  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>
「https://mochiuwiki.one/index.php/特別:携帯機器差分/13245」から取得

案内メニュー