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

26 vues (au cours des 30 derniers jours)
Jue
Jue le 13 Juil 2022
Commenté : Jue le 13 Juil 2022
%%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.]

Réponse acceptée

Alan Stevens
Alan Stevens le 13 Juil 2022
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')

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by