4th order runge-kutta method to solve two 1st order differential equations
17 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
%%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 Kommentare
Akzeptierte Antwort
Alan Stevens
am 13 Jul. 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')
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!