4th order runge-kutta method to solve two 1st order differential equations

4 ビュー (過去 30 日間)
Jue
Jue 2022 年 7 月 13 日
コメント済み: Jue 2022 年 7 月 13 日
%%4th order runge-kutta for Project 2
fy=@(t,y,z) z;
fz=@(t,y,z) (36000/(1500-40*t))-9.81-((0.1962*z^2)/(1500-40*t));
%boundary conditions
t(1) =0;
z(1) =0;
y(1)=0;
j=0;
h=1; %%step size
tfinal=10;
N=tfinal/h;
%Update loop
for j=1:N
t(j+1)=t(j)+h;
k1y=fy(t(j),y(j),z(j));
k1z=fz(t(j),y(j),z(j));
k2y=fy(t(j)+0.5*h , y(j)+0.5*h*k1y , z(j)+0.5*h*k1z);
k2z=fz(t(j)+0.5*h,y(j)+0.5*h*k1y,z(j)+0.5*h*k1z);
k3y=fy(t(j)+0.5*h,y(j)+0.5*h*k2y,z(j)+0.5*h*k2z);
k3z=fz(t(j)+h/2,y(j)+0.5*h*k2y,z(j)+0.5*h*k2z);
k4y=fy(t(j)+h,y(j)+h*k3y,z(j)+h*k3z);
k4z=fz(t(j)+h,y(j)+h*k3y,z(j)+h*k3z);
y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
fprintf('t = %6.4f, y = %6.4f\n', t, y);
so i need to display y values from t=0 to t=10 but the t output is not as expected. and y range shouldnt be a zero in it. please help
[SL changed the last line from a comment in the code to text in a text region.]

採用された回答

Alan Stevens
Alan Stevens 2022 年 7 月 13 日
Like this?
%%4th order runge-kutta for Project 2
fy=@(t,y,z) z;
fz=@(t,y,z) (36000/(1500-40*t))-9.81-((0.1962*z^2)/(1500-40*t));
%boundary conditions
z(1) =0;
y(1)=0;
h=1; %%step size
tfinal=10;
N=tfinal/h;
%Update loop
t = 0:h:N*h;
for j=1:N
k1y=fy(t(j),y(j),z(j));
k1z=fz(t(j),y(j),z(j));
k2y=fy(t(j)+0.5*h, y(j)+0.5*h*k1y, z(j)+0.5*h*k1z);
k2z=fz(t(j)+0.5*h, y(j)+0.5*h*k1y, z(j)+0.5*h*k1z);
k3y=fy(t(j)+0.5*h, y(j)+0.5*h*k2y, z(j)+0.5*h*k2z);
k3z=fz(t(j)+0.5*h, y(j)+0.5*h*k2y, z(j)+0.5*h*k2z);
k4y=fy(t(j)+h, y(j)+h*k3y, z(j)+h*k3z);
k4z=fz(t(j)+h, y(j)+h*k3y, z(j)+h*k3z);
y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
disp([' t y'])
t y
disp([t' y'])
0 0 1.0000 7.2008 2.0000 29.2186 3.0000 66.6541 4.0000 120.0703 5.0000 189.9819 6.0000 276.8443 7.0000 381.0419 8.0000 502.8774 9.0000 642.5612 10.0000 800.2017
plot(t,y,'o-',t,z,'*-'),grid
xlabel('t'), ylabel('y & z')
legend('y','z')
  1 件のコメント
Jue
Jue 2022 年 7 月 13 日
yes thank you :)

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by