4th order runge-kutta method to solve two 1st order differential equations
26 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
%%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.]
0 commentaires
Réponse acceptée
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'])
disp([t' y'])
plot(t,y,'o-',t,z,'*-'),grid
xlabel('t'), ylabel('y & z')
legend('y','z')
Plus de réponses (0)
Voir également
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!