MATLAB: Hi, I’m writing a Runge Kutta function to solve a higher order differential equation and I need to get the error to around 10^-11.

differential equationshomeworkoderunge kutta

My professor has given me a test code,
if true
  clear
% test y=sin(t)
% y^(4) = 2*y'-y"+2*y^(3)
t0 = 0.1; n = 100; h = 1e-2;
a = [sin(t0) cos(t0) -sin(t0) -cos(t0)]';
f = @(t,y) 2*y(2)-y(3)+2*y(4);
y = rk4_high_ode(a,t0,n,h,f);
ye = sin(t0+n*h);
error2 = abs(y-ye)
end
and my RK4 code is,
if true
  function y = rk4_high_ode(a, t0, n, h, f)
  v = a;
  for i = 1:n  
      k1 = f(t0,v);
      k2 = f(t0+.5*h,v+.5*k1*h);
      k3 = f(t0+.5*h,v+.5*k2*h);
      k4 = f(t0+   h,v+   k3*h);
      v = v+(h/6)*(k1+2*k2+2*k3+k4);
      t0 = t0 + h;
  end
  y = v(1);
end
end
Any tips would be appreciated. He told me that the term to be checked is the first element of the vector outputted by the function. However I'm getting 0.1562 as the error.

Best Answer

You need to rewrite your derivative function as a vector function to account for the higher order derivatives. You've got a 4x1 state vector, so you need a 4x1 derivative vector. E.g., this
f = @(t,y) 2*y(2)-y(3)+2*y(4);
should be written as this
f = @(t,y) [y(2); y(3); y(4); 2*y(2)-y(3)+2*y(4)];
This corresponds to the state vector (using prime notation) [ y , y' , y'' , y''' ]
y    -->  y(1)
y'   -->  y(2)
y''  -->  y(3)
y''' -->  y(4)
Related Question