MATLAB: Divergent 4th-order Runge-Kutta method

MATLABnumerical integrationrunge kutta

I am trying to implement a numeric integration 4th-order Runge-Kutta method applied to simple harmonic motion with following ecuations, where u :
v(n+1) = v(n) − Δt*ω^2*u(n),
u(n+1) = u(n) + Δt*v(n+1).
v is the first derivative of u and describes the acceleration whereas u describes the velocity of the oscilating system. I tried to implement the generic algorithm found in a number of posts in the forum but when I try it I found some divergency in the plot, which doesn't happen with the 2nd-order RK method nor with the Euler-Cromer method.
I used the following code introducing function handles for the equations above to avoid errors when evaluating u(n) and v(n). I based my code on following derivation of the RK4 method:
And this is the code I created based on the 2nd-order RK method:
function osc_RK4()
omega = 2; % frequency of the oscilating system
P = 2*pi/omega; % period of the oscilating system
dt = P/2000; % step size
T = 20*P; % interval considered for the study  
N_t = floor(T/dt);
t = linspace(0, N_t*dt, N_t+1);
u = zeros(N_t+1,1);
v = zeros(N_t+1,1);
% Initial condition
X_0 = 2;
u(1) = X_0;
v(1) = 0;
% functions to compute armonic movement
u_prime = @(v) v;
v_prime = @ (u,w) -w^2*u;
for n = 1:N_t
  
  v_circumflex = u_prime(v(n) + 0.5*dt*u_prime(v(n)));
  v_tilde = u_prime(v(n) + 0.5*dt*u_prime(v_circumflex));
  v_bar = u_prime(v(n) + dt*u_prime(v_tilde));
  
  u_circumflex =v_prime(u(n) + 0.5*dt*v_prime(u(n),omega), omega);
  u_tilde = v_prime(u(n) + 0.5*dt*v_prime(u_circumflex,omega), omega);
  u_bar = v_prime(u(n) + dt*v_prime(u_tilde,omega), omega);
  
  u(n+1) = u(n) + (dt/6)*(v(n) + 2*v_circumflex + 2*v_tilde + v_bar); 
  v(n+1) = v(n) + (dt/6)*(- omega^2*u(n)+ 2*u_circumflex + 2*u_tilde + u_bar); 
end
plot(t,u,'b+',t, X_0*cos(omega*t), 'k--','LineWidth',2)
I get following plot:
When running and plotting the 2nd-order RK algorithm I get following plot:
As you can see, the numeric solution matches perfectly the analyical one using the same parameters. If anyone has a clue about what am I coding wrong, thanks in advance.

Best Answer

I think your rather cumbersome notation has got you confused between the values u and v and their rates of change with time! Try replacing your loop with:
for n = 1:N_t
  
  k1v = -omega^2*u(n); 
  k1u = v(n);
  v1 = v(n) + 0.5*dt*k1v;
  u1 = u(n) + 0.5*dt*k1u;
  
  k2v = -omega^2*u1;
  k2u = v1;
  v2 = v(n) + 0.5*dt*k2v;
  u2 = u(n) + 0.5*dt*k2u;
  
  k3v = -omega^2*u2;
  k3u = v2;
  v3 = v(n) + dt*k3v;
  u3 = u(n) + dt*k3u;
  
  k4v = -omega^2*u3;
  k4u = v3;
   
  u(n+1) = u(n) + (dt/6)*(k1u + 2*k2u + 2*k3u + k4u); 
  v(n+1) = v(n) + (dt/6)*(k1v+ 2*k2v + 2*k3v + k4v); 
end
Related Question