matlab - Matlab で常微分方程式 (ODE) の系を解く方法

Question

次の形式の常微分方程式系を解かなければなりません。

dx/ds = 1/x * [y* (g + s/y) - a*x*f(x^2,y^2)]
dy/ds = 1/x * [-y * (b + y) * f()] - y/s - c

ここで、x と y は調べる必要のある変数で、s は独立変数です。残りは定数です。私はこれをode45で解決しようとしましたが、これまでのところ成功していません:

y = ode45(@yprime, s, [1 1]);

function dyds = yprime(s,y)
    global g a v0 d
    dyds_1 = 1./y(1) .*(y(2) .* (g + s ./ y(2)) - a .* y(1) .* sqrt(y(1).^2 + (v0 + y(2)).^2));
    dyds_2 = - (y(2) .* (v0 + y(2)) .* sqrt(y(1).^2 + (v0 + y(2)).^2))./y(1) - y(2)./s - d;
   dyds = [dyds_1; dyds_2];
return

@yprime には連立方程式があります。次のエラー メッセージが表示されます。

YPRIME は長さ 0 のベクトルを返しますが、初期条件ベクトルの長さは 2 です。YPRIME によって返されるベクトルと初期条件ベクトルの要素数は同じでなければなりません。

何か案は？ありがとう

Answer 1

確かに、あなたの関数を見てくださいyprime。微分状態変数の数を問題と共有する単純なモデルを使用して、この例を見てください。

function dyds = yprime(s, y)
    dyds = zeros(2, 1);
    dyds(1) = y(1) + y(2);
    dyds(2) = 0.5 * y(1);
end

yprime2 つの右辺の値を保持する列ベクトルを返さなければなりません。モデルは時間に依存しないため、入力引数sは無視できます。あなたが示す例は、dy/dt = f(t, y) の形式ではないという点でやや困難です。最初のステップとして、方程式を再配置する必要があります。と に名前を変更xするy(1)と役立ちます。yy(2)

また、あなたglobal g a v0 dは空ではありませんか？これらの変数のいずれかが初期化されていない場合、状態変数に空の行列を掛けることになり、最終的に空のベクトルdydsが返されます。これはでテストできます

assert(~isempty(v0), 'v0 not initialized');

または、yprimeデバッグ ブレークポイントを使用できます。

Answer 2

ODE ソルバーの構文は次のとおりです。[s y]=ode45(@yprime, [1 10], [2 2])

そして、あなたのケースでは要素ごとの操作を行う必要はありません。つまり、.*単に使用する代わりに*