連載
» 2015年02月20日 10時00分 公開

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

[伊藤孝宏,MONOist]
前のページへ 1|2       

無名関数を用いる方法

 無名関数を用いる方法では、任意の関数名=@(独立変数、従属変数)の後に式を記載します。連立微分方程式のように、従属変数が複数ある場合、従属変数を配列とします。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」赤間世紀著、工学社刊

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

FreeMat

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


筆者紹介

伊藤孝宏(いとう・たかひろ)

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



前のページへ 1|2       

Copyright © ITmedia, Inc. All Rights Reserved.