Divergent 4th-order Runge-Kutta method
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Carlos Rueda
el 13 de Feb. de 2021
Comentada: Carlos Rueda
el 14 de Feb. de 2021
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/518407/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/518412/image.png)
When running and plotting the 2nd-order RK algorithm I get following plot:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/518417/image.png)
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.
0 comentarios
Respuesta aceptada
Alan Stevens
el 13 de Feb. de 2021
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
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!