MATLAB: Incorrect results using Runge Kutta 5th order method

differential equationsrunge kutta

Hello, I have written the following code for a second order differential equation, split into two first order differential equations, the code itself seems to work however after the first step all results for z and t are 0 which is incorrect, is there any errors in my code/working?
Original equation: d^2t/dx^2 = (h*p/k*ac)*(t-ta)
Split equations: dt/dx = z
dz/dx = (h*p/k*ac)*(t-ta)
5th order Runge Kutta:
ac=7.854*10^-7 %cross sectional area
k=110; % Conductive heat transfer coefficient
ta=25; % Room temperature
h=20; %Convective heat transfer coefficient
p=0.04063; %Perimeter of cylinder
s=0.001 %step size
xfinal=0.02; % final distance
N=xfinal/s;
x(1)=0; %initial values
t(1)=50;
z(1)=0;
dtdx=@(x,t,z) z;
dzdx=@(x,t,z) (h*p/k*ac)*(t-ta);
for i=1:N
k1= dtdx(x(i),t(i),z(i));
l1=dzdx(x(i),t(i),z(i));
k2=dtdx(x(i)+0.25*s,t(i)+0.25*k1*s,z(i)+0.25*l1*s)
l2=dzdx(x(i)+0.25*s,t(i)+0.25*k1*s,z(i)+0.25*l1*s)
k3=dtdx(x(i)+0.25*s,t(i)+0.125*k1*s+0.125*k2*s,z(i)+0.125*l1*s+0.125*l2*s)
l3=dzdx(x(i)+0.25*s,t(i)+0.125*k1*s+0.125*k2*s,z(i)+0.125*l1*s+0.125*l2*s)
k4=dtdx(x(i)+0.5*s,t(i)-0.5*k2*s+k3*s,z(i)-0.5*l2*s+l3*s)
l4=dzdx(x(i)+0.5*s,t(i)-0.5*k2*s+k3*s,z(i)-0.5*l2*s+l3*s)
k5=dtdx(x(i)+0.75*s,t(i)+0.1875*k1*s+0.5625*k4*s,z(i)+0.1875*l1*s+0.5625*l4*s)
l5=dzdx(x(i)+0.75*s,t(i)+0.1875*k1*s+0.5625*k4*s,z(i)+0.1875*l1*s+0.5625*l4*s)
k6=dtdx(x(i)+s,t(i)-0.4286*k1*s+0.2857*k2*s+1.7143*k3*s-1.7143*k4*s+1.1429*k5*s,z(i)-0.4286*l1*s+0.2857*l2*s+1.7143*l3*s-1.7143*l4*s+1.1429*l5*s)
l6=dzdx(x(i)+s,t(i)-0.4286*k1*s+0.2857*k2*s+1.7143*k3*s-1.7143*k4*s+1.1429*k5*s,z(i)-0.4286*l1*s+0.2857*l2*s+1.7143*l3*s-1.7143*l4*s+1.1429*l5*s)
x(i+1)=x(i)+s;
t(i+1)=t(i)*0.0111*(7*k1+32*k3+12*k4+32*k5+7*k6)*h
z(i+1)=z(i)*0.0111*(7*l1+32*l3+12*l4+32*l5+7*l6)*h
end

Best Answer

I have already answered this on your other thread, but will re-post here. The 'h' variable in your t and z updates is supposed to be stepsize, but you are confusing it with your heat parameter that is used in your derivative equations. So you probably need to use 's' here. Also, I would expect your updates to be added to the current values, not multiplied by your current values. So maybe something like this is what you need:
t(i+1) = t(i) + 0.0111*( 7*k1 + 32*k3 + 12*k4 + 32*k5 + 7*k6 )*s;
z(i+1) = z(i) + 0.0111*( 7*l1 + 32*l3 + 12*l4 + 32*l5 + 7*l6 )*s;
Related Question