| Mathematics | |
ODE 問題の表現
この節は、MATLAB ODE ソルバの一つを使って、初期値 ODE 問題を解く過程を記述します。van der Pol 方程式を使います。
| 注意 より詳細な情報は、ODEソルバの書式設定法 を参照してください。 |
例題：MATLABでIVP ODE問題を解く(van der Pol
方程式、ノンスティッフ)
この例題は、MATLABを使って、初期値 ODE 問題を解くために必要なステップを示し、説明します。
1 問題を1階のODE問題として書き換えます。 ODEsは、1階以上の微係数と複数の従属変数を含んでいます。MATLAB ODEソルバを使って、このような方程式をつぎの型をした一階微分方程式の等価なシステムとして表現します。
つぎのように代入して、1階の微分方程式として、記述できます。
結果は、
個の1階の ODEs の等価なシステムになります。
たとえば、van der Pol 方程式 (二階微分)を書き換えます。
ここで、
は、スカラパラメータです。そして、
を代入して、つぎの1階ODEsの結果を得ます。
2 MATLABで、1階のODEsのシステムをコード化 1階のODEsのシステムとして方程式を表現すると、MATLAB ODE ソルバが使用できる型にそれをコード化します。関数は、つぎの型になります。
dydt = odefun(t,y)
t と yが、関数の最初の2つの引数になる必要がありますが、関数は、それらを必ず使用するとは限りません。出力
dydtは、列ベクトルで、yの微係数です。
下のコードは、MATLAB 関数vdp1で表したvan
der Polシステムです。関数vdp1 は、
を仮定しています。 
と 
は、2要素ベクトルのy(1) と y(2)
になります。
function dydt = vdp1(t,y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
vdp1は、引数t と yを受け入れますが、計算では、tを使いません。
3 ソルバを問題に適用 問題を解くために使用するソルバを決定します。そして、ソルバをコールし、1階のODEsシステムを記述するために作成した関数を渡し、それと共に、問題を解くための時間間隔と初期条件ベクトルを渡します。ODE ソルバの記述に関しては、初期値問題のソルバ とODE ソルバリファレンスページ を参照してください。
。van der Pol システムに対して、初期値y(1) = 2
と y(2) = 0を使って、時間区間[0 20]で、ode45を使用します。
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
この例題は、ode45への関数ハンドルとして、@を使って、vdp1に渡します。出力結果は、時間点tを表す列ベクトルと解の配列yです。計算結果の行列yの各列は、列ベクトルtで出力される時刻に対応します。yの最初の列は、
に対応し、2番目の列は、
に対応します。
注意：
関数ハンドルに関する情報に関しては、MATLABドキュメントの中のfunction_handle (@), func2str, str2func
のリファレンスページとM-ファイルプログラミングの節を参照してください。
|
4 ソルバ出力の可視化 plot
コマンドを使って、ソルバ出力を可視化することができます。
plot(t,y(:,1),'-',t,y(:,2),'--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2')
また、ソルバの出力関数を使って、出力を処理することもできます。ソルバは、連続的な時間ステップの各ステップの後で、積分プロパティOutputFcnの中に設定し、関数をコールします。odesetを使って、希望する関数をOutputFcnに設定します。詳細は、OutputFcn を参照してください。
付加的なパラメータをODE 関数に渡す
関数は、options引数の後に任意の入力引数を設定して、ODE関数やユーザがoptionsに設定した関数に渡します。たとえば、
muを明示的に設定するのではなく、muパラメータを渡すことにより、van
der Pol 関数を一般化します。 function dydt = vdp1(t,y,mu) dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
options引数の後に、設定することで、パラメータmuを関数vdp1に渡します。この例題では、options
= []を使っています。[t,y] = ode45(@vdp1,tspan,y0,[],mu)
これらの関数をベースにした例題は、vdpode を参照してください。
例題：van der Pol 方程式、 µ = 1000 (スティッフ)
| 注意 MATLABの中で、初期値 ODE 問題を解くための詳細な例は、例題：MATLAB の中で、IVP ODE を解く (van der Pol 方程式、ノンスティッフ)を参照してください。 |
この例題は、スティッフな問題を取り扱います。スティッフな問題では、シミュレーション時間に比べ、非常に短い時間スケールで解は変化しますが、興味のある解(適切な精度を保った数値解)は、はるかに大きい時間スケールでしか変化しません。スティッフ問題用に設計されていない方法では、最も速い変化を解くために十分な細かさの時間間隔を使うために、時間変化のゆっくりしている部分では、効率の悪いものになります。
を大きくして1000にすると、van
der Pol方程式の解は、急激に変化し、より長い時間スケールで振動を示します。初期値問題の解を近似することは、非常に難しい問題です。この特別な問題はスティッフなので、ode45のようなスティッフに対応しないソルバは、非現実的な効率の悪いものになります。スティッフなソルバode15sは、このような問題に対して作られたものです。
つぎのコードvdp1000は、
= 1000のとき、ODEファイルの中にvan del
Polシステムを、どのように表現するかを示したものです。
function dydt = vdp1000(t,y) dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
注意 この例題は、ODE関数の中に、 をハードコードしています。vdpode
デモは、同じ問題を解いていますが、ODE
関数への付加的な引数として、ユーザ定義の  を渡しています。詳細は、 付加的な ODE ソルバ引数をとしてください。
|
さて、関数ode15sを使って、vdp1000を解きましょう。[2;
0]の初期条件ベクトルはそのままですが、計算時間を[0
3000]にします。スケーリングの目的(積分ステップ間隔とシステムの振幅周期を調べる)で、y(t)の最初の成分をプロットしましょう。
[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]);
plot(t,y(:,1),'-');
title('Solution of van der Pol Equation, \mu = 1000');
xlabel('time t');
ylabel('solution y_1');
| |
初期値問題ソルバ | ODE ソルバ性能の改良 | |
