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;
P = 2*pi/omega;
dt = P/2000;
T = 20*P;
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);
X_0 = 2;
u(1) = X_0;
v(1) = 0;
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.
