時変のパラメータを持​つ微分方程式を解くに​はどうすればよいです​か?

26 ビュー (過去 30 日間)
MathWorks Support Team
MathWorks Support Team 2009 年 6 月 29 日
下記の微分方程式で、時間ベクトルとパラメータ値が1対1で定義されているパラメータp(t)を含む、微分方程式があります。
y'(t) + p(t)*y(t) = 0, y(0) = 1
 
t p(t)
-------------------
0 0
: :
π/2 1
このような時変のパラメータを持つ微分方程式を解くにはどうすればよいですか?

採用された回答

MathWorks Support Team
MathWorks Support Team 2009 年 6 月 29 日
入れ子関数または無名関数を用いて時間ベクトルとパラメータ値を導関数へ引き渡します。
導関数を定義している関数内では、interp1関数を用いて時刻 t におけるパラメータ値を補間します。
入れ子関数、無名関数に関する説明は、下記にある
  「関連ソリューション:関数の入力引数として指定する関数にパラメータ値を与えるにはどのようにすればよいですか。」
をご覧下さい。
なお、以下のサンプルプログラムでは、時間ベクトルと対になるパラメータ値を定義してから、解析を行っています。
1. 入れ子関数
入れ子関数を使用すると、導関数を定義している関数内で時間ベクトルとパラメータ変数をそのまま使用することができます。
function [T1,Y1] = ode_param_sample1()
TSPAN = [0 pi/2]; % 計算開始・終了時間(TSPAN)、初期状態(IC)の定義
IC = 1;
tv = linspace(0,pi/2,10); % 時間ベクトル(tv)と対応するパラメータ値(ftv)を定義
ftv = sin(tv);
[T1 Y1] = ode45(@myODE1, TSPAN, IC); % 常微分方程式を解く
% 導関数の定義(入れ子関数)
function dydt1 = myODE1(t1,y1)
ft = interp1(tv, ftv, t1,'spline'); % 時刻 t におけるパラメータの値を補間
dydt1 = - ft *y1; % 時刻 t における導関数の値を計算
end
end
なお、上記をMATLABファイルとして保存し、下記コマンドを実行することで、解析結果が得られます。
>> [T1,Y1] = ode_param_sample1;
2. 無名関数
無名関数の形式を用いて、ODE関数に時間ベクトルとパラメータ値の引渡しを行ないます。
また、導関数を定義している関数では、時間ベクトルとパラメータ値の両方が入力引数となるように定義します。
function [T2,Y2] = ode_param_sample2()
TSPAN = [0 pi/2]; % 計算開始・終了時間(TSPAN)、初期状態(IC)の定義
IC = 1;
tv = linspace(0,pi/2,10); % 時間ベクトル(tv)と対応するパラメータ値(ftv)を定義
ftv = sin(tv);
% 常微分方程式を解く
[T2 Y2] = ode45(@(t,y) myODE2(t,y,tv,ftv), TSPAN, IC); % 無名関数の記述
function dydt2 = myODE2(t2,y2,tv,ftv) % 導関数
ft = interp1(tv, ftv, t2,'spline'); % 時刻 t におけるパラメータの値を補間
dydt2 = - ft *y2; % 時刻 t における導関数の値を計算
なお、上記をMATLABファイルとして保存し、下記コマンドを実行することで、解析結果が得られます。
>> [T2,Y2] = ode_param_sample2;

その他の回答 (0 件)

カテゴリ

Help Center および File Exchange数値積分と微分方程式 についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!