はじめに
matlabで今回、オイラー法による微分方程式の解法を実装しなさいとの課題が大学院で出ていたので、同じような課題の人もいるんじゃないかと思ってn番煎じだとは思いますが記載します。
コードは下記に挙げているのでよかったらどうぞ!
そもそも、オイラー法ってなに?
オイラー法:初期値問題の数値解法
オイラー法は初期値問題(特定の時点での値とその変化率を表す微分方程式が与えられている問題)を数値的に解く最も単純な手法の一つです。
オイラー法は以下のように機能します:
- 初期設定:解を探している微分方程式と初期値(特定の時点での値)を入力します。
- 時間ステップの設定:解を見つけたい時間範囲を小さな時間ステップに分割します。
- 値の更新:各時間ステップにおいて、微分方程式の右辺(つまり、変化率)を使用して、その時点での値から次の時点での値を予測します。
数学的には、オイラー法は以下の反復公式によって定義されます:
y_{n+1} = y_n + h*f(t_n, y_n)
ここで、
-
y_{n+1}
は次のステップでの解の近似値, -
y_n
は現在のステップでの解の近似値, -
h
は時間ステップの大きさ(一定), -
f(t_n, y_n)
は微分方程式の右辺を表す関数。
オイラー法は直感的で実装が容易な一方で、誤差が累積しやすく、精度がそれほど高くないという欠点があります。したがって、より複雑な問題に対しては、より高度な手法(例えばルンゲ・クッタ法)がしばしば用いられます。
オイラー法: 微分方程式の数値解法
オイラー法は初期値問題を解くための最も基本的な数値解法です。初期値問題とは、微分方程式とその解の初期値が与えられ、それに基づいて解を見つける問題です。
オイラー法の基本概念
オイラー法は、微分方程式の解を近似するための逐次的な手法です。与えられた微分方程式と初期値から始めて、小さな時間ステップで進むことで、解の近似値を生成します。
具体的には、時間ステップ h
と初期値 y(t0)=y0
が与えられたとき、次のステップ t0+h
における解 y(t0+h)
の近似値 y1
は次の式で計算されます:
y1 = y0 + h*f(t0, y0)
ここで、f(t, y)
は微分方程式の右辺を表します。
オイラー法のサンプルプログラム(MATLAB)
以下に、オイラー法を使用して微分方程式 dy/dx = y - t^2 + 1
を解くMATLABプログラムを示します。この微分方程式の解は (t+1)^2 - 0.5*exp(t)
です。
% ODEの右辺の関数と真の解の関数の定義
f = @(t, y) y - t^2 + 1.0;
g = @(t) (t+1)^2 - 0.5*exp(t);
% オイラー法の実行
[t, w, error] = eulerODE2(f, g);
% 近似解と真の解のプロット
plot(t, w, 'r', t, g(t), 'b');
legend('Euler Approximation', 'True solution');
xlabel('t');
ylabel('y');
% 誤差のプロット
figure;
plot(t, error);
xlabel('t');
ylabel('Error');
function [t, w, error] = eulerODE2(f, g)
% 入力:
% f - ODEの右辺の関数
% g - 真の解の関数
% 出力:
% t - 時間ベクトル
% w - 近似解ベクトル
% error - 真の解との誤差ベクトル
% ユーザーからの入力を受け取る
a = input('Enter the start of the interval, a: ');
b = input('Enter the end of the interval, b: ');
alpha = input('Enter the initial value at a, alpha: ');
N = input('Enter the number of steps, N: ');
h = (b-a)/N; % ステップ幅
t = a:h:b; % 時間ベクトル
w = zeros(1, N+1); % 近似解を格納するベクトル
w(1) = alpha; % 初期値の設定
% オイラー法
for i = 1:N
w(i+1) = w(i) + h*f(t(i), w(i));
end
% 真の解との誤差
true_sol = g(t); % 真の解ベクトル
error = abs(true_sol - w); % 誤差ベクトル
end
これはオイラー法の基本的な説明と実装例です。より複雑な微分方程式や他の数値解法についても同様の方法で解くことができます。
出力はこんな感じ
近似解と真の解のプロット
最後に
matlab初めて初心者で綺麗なコードではないでしょうが、誰かの役に立てば幸いです。
pythonでの実装は下記の方のがわかりやすくいいかと思います。