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

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')

Answer 1

[Paul]

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.

[Sam Chak]

I just change one thing:

h = 0.01;