Matlabで微分方程式の計算練習

今回は簡単にMatlabで微分方程式を解く練習をします。

対象とする微分方程式は以下とします。

$$ u'' + u = 0 $$

この理論解の１つは以下のようになり、u-u'平面でのグラフは円になります。 $$ u = \cos{t} $$ $$ u'=\sin{t} $$

この方程式をMatlabのODEソルバーで解くには、まず以下のように方程式を関数で定義してやります。

function dydt = odefcn(t,y)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -y(1);

このファイルをodefcn.mという名前で保存します。

そして別に以下のようなスクリプトを書きます。

tspan = [0 10];
y0 = [1 0];
[t1,y1] = ode45(@(t,y) odefcn(t,y), tspan, y0);
plot(y1(:,1),y1(:,2),'-o')

このファイル（例えばodetest.m)を実行してやると以下のようなグラフができます。

MatlabのODEソルバーについては別の機会に説明しますが、ode45は精度が高く、例えばタイムステップを10から1000に増やしても、円周軌道を維持し続けます。

tspan = [0 1000];
y0 = [1 0];
[t1,y1] = ode45(@(t,y) odefcn(t,y), tspan, y0);
plot(y1(:,1),y1(:,2),'-o')

一方、例えば精度がode45ほどではないode23では、1000ステップ計算すると軌道がどんどん円周から外れていってしまいます。

tspan = [0 1000];
y0 = [1 0];
[t1,y1] = ode23(@(t,y) odefcn(t,y), tspan, y0);
plot(y1(:,1),y1(:,2),'-o')

modelicaで計算する際は特に精度に特に気を遣っていないのですが、本当は注意しなくてはいけませんね。今回はここまでにします。