C Sharpと数値解析 - オイラー法
ナビゲーションに移動
検索に移動
概要
前進差分法
前進差分法 (Forward Euler Method) は、精度は1次であるが、最も簡単な実装であり、現在の点での傾きを用いて次の点を予測する。
これは、 に基づいている。
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);
後退差分法 (ニュートン法を使用しない方法)
ニュートン法を使用しない後退差分法 (Backward Euler Method) は、次式で示す点での傾きを使用して計算を行う陰的な方法である。
に基づいている。
現在の値を使用して次のステップを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) は、次式で示す点での傾きを使用して計算を行う陰的な方法である。
ニュートン法による反復計算が必要で計算コストは高くなるが、精度および安定性は高い。
に基づいている。
ニュートン法では、収束判定のための許容誤差を設定して、最大反復回数による安全装置を実装する必要がある。
また、以下の例では、初期推定値として前進オイラー法の結果を使用している。
ただし、収束性は初期推定値に依存する場合があることに注意する。
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);
中心差分法
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);