Divergent 4th-order Runge-Kutta method

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.

Answer 1

[Alan Stevens]

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

[Carlos Rueda]

Thank you! That totally worked. The funny thing is, I had already tried a notation without function handles with the same wrong plot...