0

微分方程式をプロットするのに助けが必要です...それはすべてファンキーになり続け、グラフは本来あるべき姿とは異なります。

function [dydt] = diff(y,t)

dydt = (-3*y)+(t*(exp(-3*t)));

end

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.5;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = diff(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

これは私のコードですが、グラフは本来あるべき姿とはまったく異なります。forステートメントのインデックスのようなものが欠けていますか?

編集

エラーが発生します：

Error using diff
Difference order N must be a positive integer scalar.

Error in diff3 (line 12)
    yPrime = diff(t(k-1),y(k-1));
4

1 に答える 1

1

Danil Asotsky と horchler がコメントで指摘したエラーを修正した後:

  1. 組み込み関数「diff」との名前の競合を回避する
  2. 引数の順序を に変更しt,yます。
  3. dt時間ステップを 0.1 に減らす
  4. ODE 右辺を無名関数に変換する

(そして関数定義で不要な括弧を削除します)、コードは次のようになります。

F = @(t,y) -3*y+t*exp(-3*t);

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.1;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = F(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

これは期待どおりに機能します:

出力

于 2015-01-05T21:13:19.937 に答える