please help to fix my code for 4th order Runge Kutta method for second order ODE.

2 views (last 30 days)
Siti Nafisah Bekti Lelana
Siti Nafisah Bekti Lelana on 11 Jun 2022
Edited: Paul on 11 Jun 2022
Hi, im trying to solve second order ODE by using 4th order Runge Kutta method.
I start by converting it into two first order ode whereby dy/dt=z and dz/dt= (21285 + 392.4t - 0.1962z^2) / (1500 - 40*t)
with initial condition t0=0, y0=0, z0=0
step size, h=1.
i dont know why i never get the code to produce the correct output.
this is my coding:
% Define function handles
fy=@(t,y,z) z;
fz=@(t,y,z) (21285 + (392.4*t) - (0.1962*(z^2))) / (1500 - 40*t);
%Initial conditions
t(1)=0;
y(1)=0;
z(1)=0;
%Step size
h=1;
tfinal=10;
N=tfinal/h;
%Update loop
for i=1:N
%Update time
t(i+1)=t(i)+h;
%Update Runge Kutta
k1y=fy(t(i),y(i),z(i));
k1z=fz(t(i),y(i),z(i));
k2y=fy(t(i)+h/2,y(i)+h/2*k1y,z(i)+h/2*k1z);
k2z=fz(t(i)+h/2,y(i)+h/2*k1y,z(i)+h/2*k1z);
k3y=fy(t(i)+h/2,y(i)+h/2*k2y,z(i)+h/2*k2z);
k3z=fz(t(i)+h/2,y(i)+h/2*k2y,z(i)+h/2*k2z);
k4y=fy(t(i)+h,y(i)+h*k3y,z(i)+h/k3z);
k4z=fz(t(i)+h,y(i)+h*k3y,z(i)+h/k3z);
y(i+1)=y(i)+h/6*(k1y + 2*k2y +2*k3y + k4y);
z(i+1)=z(i)+h/6*(k1z + 2*k2z +2*k3z + k4z);
end
%Plot solution
plot (t,y)
xlabel ('t')
ylabel ('y')
legend ('Time')

Answers (1)

Paul
Paul on 11 Jun 2022
Edited: Paul on 11 Jun 2022
Hi Siti Nafisah Bekti Lelana,
Why do you think it's the wrong answer? I didn't check the equations, but the code seems to yield a result very close to that from ode45, and I suspect it would be closer if the step size, h, is reduced. Having said that, the code could be simplified by combining y and z into a single vector variable, e.g., x, and then just applying the ode4 equations to x.
fy=@(t,y,z) z;
fz=@(t,y,z) (21285 + (392.4*t) - (0.1962*(z^2))) / (1500 - 40*t);
%Initial conditions
t(1)=0;
y(1)=0;
z(1)=0;
%Step size
h=1;
tfinal=10;
N=tfinal/h;
%Update loop
for i=1:N
%Update time
t(i+1)=t(i)+h;
%Update Runge Kutta
k1y=fy(t(i),y(i),z(i));
k1z=fz(t(i),y(i),z(i));
k2y=fy(t(i)+h/2,y(i)+h/2*k1y,z(i)+h/2*k1z);
k2z=fz(t(i)+h/2,y(i)+h/2*k1y,z(i)+h/2*k1z);
k3y=fy(t(i)+h/2,y(i)+h/2*k2y,z(i)+h/2*k2z);
k3z=fz(t(i)+h/2,y(i)+h/2*k2y,z(i)+h/2*k2z);
% k4y=fy(t(i)+h,y(i)+h*k3y,z(i)+h/k3z);
% k4z=fz(t(i)+h,y(i)+h*k3y,z(i)+h/k3z);
k4y=fy(t(i)+h,y(i)+h*k3y,z(i)+h*k3z);
k4z=fz(t(i)+h,y(i)+h*k3y,z(i)+h*k3z);
y(i+1)=y(i)+h/6*(k1y + 2*k2y +2*k3y + k4y);
z(i+1)=z(i)+h/6*(k1z + 2*k2z +2*k3z + k4z);
end
%Plot solution
plot (t,y,t,z)
xlabel ('t')
ylabel ('y and z')
legend ('y','z')
[tt,xx] = ode45(@(t,x)([fy(t,x(1),x(2));fz(t,x(1),x(2))]),[0 10],[0 0]);
plot(tt,xx)
Edit: I just noticed that k4 equations appear to have an error. Shoudln't the they both have h*k3z instead of h/k3z? See in-code edit above.

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!