ODE45で常微分方程式の初期値問題を解く：無償ソフトで技術計算しよう【シミュレーション基礎編】（3）（2/2 ページ）

[伊藤孝宏，MONOist]

無名関数を用いる方法

　無名関数を用いる方法では、任意の関数名＝＠（独立変数、従属変数）の後に式を記載します。連立微分方程式のように、従属変数が複数ある場合、従属変数を配列とします。ode45には定義した関数名を、関数名に相当する引数の位置に記載すれば、定義された関数が渡されます。例えば、下記のように記載すると、cos(x)がode45に渡されます。

--> f=@(x) cos(x);
--> [x,y]=ode45(f,[0,pi],0);

　無名関数を用いる方法では、ode45内で直接定義することが可能です。この場合、関数名を除いた部分を、ode45の関数名に相当する引数の位置に記載します。例えば、下記のように記載すると、ode45にcos(x)が渡されます。ヘルプファイルの説明が、この方法に該当します。

--> [x,y]=ode45(@(x) cos(x),[0,pi],0);

関数mファイルを用いる方法

　係数を別の行に記載したい場合、複数の行に渡って記載することになります。そうした場合、関数mファイルを用います。

　1行目をfunction　変数＝関数名（独立変数、従属変数）として、2行目に変数＝式を記載します。次に、関数名を名称として、mファイルに保存します。ode45に渡すには、関数名に相当する引数の位置に＠関数名として記載します。関数mファイルでは、関数を入れ子構造にできるため、例では、ex403という関数mファイルの中に、funtion f(x,y)としてcos(x)を定義しています。実行させるには、コマンドウィンドウでex403と入力します。

function ex403
    [x,y]=ode45(@f,[0,pi],0);
    plot(x,y,'o-');grid('on');
function xy=f(x)
    xy=cos(x);

ex403.m

＞＞「ex403.m」ダウンロード

　ここで、ex403.mをFreeMatで計算した結果とMATLABで計算した結果を参考までに下記に示します。図3がFreeMatで解いた結果で、図4がMATLABで解いた結果です。MATLABの方が分割を細かくし、より滑らかな結果が出るようにしていることが分かります。

図3：FreeMatで解いた結果

図4：MATLABで解いた結果

連立微分方程式の解法

　ode45は連立微分方程式も解くことが可能です。例として、下記の連立微分方程式を初期条件x(0)=1、y(0)=1で解いてみます。

　参考までに解析解を求めてみます。MATLAB（Mupad）で解いてみると、下記に示すように、結構複雑な式になります。このように、微分方程式を解析的に解くことは困難を伴う場合が多いです。一方、工学分野では、結果が数値で得られれば十分というものも多く、そういった場合、数値解法が威力を発揮します。次に、数値的に解いてみます。

　無名関数をode45内で定義すれば、下記のように1行で計算できます。この場合、関数の与え方は、fに結果が返るとすると、xに相当するのはf(1)、yに相当するのはf(2)と、変数をfの配列にします。次に、式をセミコロン;でつなぎ、行方向に並べた配列とします。同じく初期条件も[1;1]と行方向に並べた配列とします。

--> [t,f]=ode45(@(t,f) [f(1)+2*f(2);3*f(1)+4*f(2)],[0,1],[1;1]);

という式を解いた結果が、tは行方向の配列で、fはn行2列の配列として返ってきます。fの1列目にはx、fの2列目にはyが格納されています。下記のように入力してグラフ化すれば、図5の結果が得られます。

--> plot(t,f);grid('on');xlim([0,1.4]);
--> legend('t_x','t_y');
--> xlabel('t');ylabel('x,y');

図5：dx/dt=x+2*y、dy/dt=3*x+4*yの数値解法結果

　「legend」はプロットごとに名称を付けるコマンドで、プロットされた順にプロットの種類とlegend内に記載された名称が並んで表示されます。

　数値として結果を取り出すには、x=f(:,1);y=f(:,2);などとします。

---

　次回は、応用編としてode45を使ってさまざまな微分方程式を解いてみます。

参考文献

• 「MATLABハンドブック」小林一行著、秀和システム刊
• 「はじめてのFreeMat」赤間世紀著、工学社刊

無償ソフトで技術計算しよう

無償ソフトで技術計算しよう

無償の工学計算ソフトでも、かなり高度な計算ができる！

筆者紹介

伊藤孝宏（いとう・たかひろ）

1960年生。小型モーターメーカーのエンジニア。博士（工学）。専門は流体工学、音・振動工学。現在は、LabVIEWを使って、音不良の計測・診断ソフト、特性自動検査装置などの開発を行っている。