ODEソルバの使い方

ODEソルバの使い方
担当：柳澤優香
常微分方程式(ODE)の数値解法
•
単独一階常微分方程式
𝑑𝑦
= 𝑓 𝑡, 𝑦
𝑑𝑡
が,f(t,y)が点(t0,y0)を含むある二次元領域で有界かつ連続であれば,
その近傍で初期条件y(t0)=y0を満たす解を少なくとも１つ持つ.
解の存在と一意性が保証されても, その厳密解を求めることは難しく, 適当な
分点t1,t2,…におけるy(t1),y(t2)…を求めることになる．
⇨ これを数値的に解くという
<代表的な数値解法>
• オイラー(Euler)法, 修正オイラー法
• ルンゲ・クッタ(Runge-Kutta)法
1
オイラー(Euler)法
• 刻み幅 h (h>0)を設定し下の計算を繰り返す.
注)精度を上げようと, hを小さくしすぎると丸め誤差の影響を受けて誤差がかえって増大する
• 初期値を t0, y0 とすると
のように、t の値を t0 から h ずつ増やしながら
y0 から yn の値を数値的に求める.
2
実行例:Euler法
dy/dt = y, t0=0, y0=1
t0=0; t_end=5; N=100;
t=linspace(t0, t_end, N+1); % 刻む数を指定（t0とt_endの間を100刻み）
h = (t_end - t0)/N;
%刻み幅
y(1) = 1;
%yの初期値
for n=1:N, y(n+1)= y(n) + h*f(y(n),t(n)); end % オイラー法
plot(t,y,'-.r',t,exp(t),'-b')
functiondy =expo(y,t)
dy=y;
end
解析解 y=et と比較
3
ODEソルバ(solver)
• MATLABでは、微分方程式の解を求める組み込み関
数として,odeソルバが用意されている
4次のルンゲ・クッタ法を使用
注）各ODEソルバで使われるアルゴリズムは、精度と設計されたシステムのタイプ（stiffの度合い）により異なる。
堅い方程式(stiffequation)とは、解の緩やかな変化と急激な変化が混在する方程式で、刻み幅を小さく設定しな
ければならない場合がある。
4
ルンゲ・クッタ(Runge-Kutta)法
• 微分方程式の数値解法として、4次のルンゲ・クッタ法を使用。
𝑡'() = 𝑡' + ℎ
𝑘) = ℎ . 𝑓(𝑡' , 𝑦' )
ℎ
𝑘)
𝑘- = ℎ . 𝑓(𝑡' + , 𝑦' + )
2
2
ℎ
𝑘𝑘2 = ℎ . 𝑓(𝑡' + , 𝑦' + )
2
2
𝑘3 = ℎ . 𝑓(𝑡' + ℎ, 𝑦' + 𝑘2)
𝑦'()
1
= 𝑦' + (𝑘) + 2𝑘- + 2𝑘2 +𝑘3 )
6
5
実行例:Runge-Kutta法
functiondy=rungef(t,y)
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=y(3);
dy(3)=-4*y(3)-5*y(2)-2*y(1);
end
tspan=[03];%積分範囲
y0=[5;-8;16];%初期値を列ベクトルで与える。
[t,y]=ode45('rungef',tspan,y0);
plot (t,y(:,1));
6

Download Report