3

常微分方程式系..初期値問題..時間または独立変数に依存するパラメータをどのように解くのですか? 私が持っている方程式を言う

 Dy(1)/dt=a(t)*y(1)+b(t)*y(2);
 Dy(2)/dt=-a(t)*y(3)+b(t)*y(1);
 Dy(3)/dt=a(t)*y(2);

ここで、a(t) はベクトルで、b(t) =c*a(t); ここで、a と b の値は単調ではなく時間ステップごとに時間とともに変化します。この投稿を使用してこれを解決しようとしました....しかし、同じ原則を適用すると...エラーメッセージが表示されました

「griddedInterpolant を使用したエラー ポイントの座標が厳密な単調な順序で配列されていません。」

誰かが私を助けてくれますか？

4

1 に答える 1

4

答えの最初の部分または2番目の部分があなたに関連しているかどうかを確認するために最後まで読んでください：

パート1：最初に、計算を説明する関数と.mを与える関数を 含むファイルを作成します。例：次のコードを含むというファイルを作成します。abfun_name.m

function Dy = fun_name(t,y)

Dy=[ a(t)*y(1)+b(t)*y(2); ...
    -a(t)*y(3)+b(t)*y(1); ...
     a(t)*y(2)] ;
end

    function fa=a(t);
        fa=cos(t); % or place whatever you want to place for a(t)..
    end

    function fb=b(t);
        fb=sin(t);  % or place whatever you want to place for b(t)..
    end

次に、次のコードを含む2番目のファイルを使用します。

t_values=linspace(0,10,101); % the time vector you want to use, or use tspan type vector, [0 10]
initial_cond=[1 ; 0 ; 0];  
[tv,Yv]=ode45('fun_name',t_values,initial_cond);
plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o');
legend('y1','y2','y3');

もちろん、fun_name.m私が書いた場合、a(t)とのサブ関数を使用する必要はありません。可能であればb(t)、明示的な関数形式を使用できます（など）。Dycos(t)

パート2 ：、が（パート1のように）の関数として表現できない、たまたま持っている数のベクトルである場合 a(t)、それぞれが発生する時間ベクトルも必要になります。これは次のようになります。もちろん、ODEに使用するのと同じ時間ですが、補間が機能する限り、そうである必要はありません。期間や解像度が異なる場合は、一般的なケースを扱います。次に、次のことを実行して、ファイルを作成します。b(t)tfun_name.m

function Dy = fun_name(t, y, at, a, bt, b) 

a = interp1(at, a, t); % Interpolate the data set (at, a) at times t
b = interp1(at, b, t); % Interpolate the data set (bt, b) at times t

Dy=[ a*y(1)+b*y(2); ...
    -a*y(3)+b*y(1); ...
     a*y(2)] ;

これを使用するには、次のスクリプトを参照してください。

%generate bogus `a` ad `b` function vectors with different time vectors `at` and `bt`

at= linspace(-1, 11, 74); % Generate t for a in a generic case where their time span and sampling can be different
bt= linspace(-3, 33, 122); % Generate t for b    
a=rand(numel(at,1));
b=rand(numel(bt,1));
% or use those you have, but you also need to pass their time info...

t_values=linspace(0,10,101); % the time vector you want to use 
initial_cond=[1 ; 0 ; 0];  
[tv,Yv]= ode45(@(t,y) fun_name(t, y, at, a, bt, b), t_values, initial_cond); % 
plot(tv,Yv(:,1),'+',tv,Yv(:,2),'x',tv,Yv(:,3),'o');
legend('y1','y2','y3');
于 2013-01-15T22:58:27.843 に答える