Having issues with RK4 on small angle pendulum
显示 更早的评论
I am trying to implemnt and compaire some numerical solvers, but having issues when it comes to the RK4 as its not behaving as i would expect.
I would expect it to be stable and follow the analytic soultion, but it becomes unstable. And it cant figure out whats wrong in my rk4 code.
equations to be solved:
omega_dot = -(g/l)*theta
theta_dot=omega;
clf
clear all
%Declearing stuff
l = 1;
m = 1;
g = 9.81;
theta(1) = 0.2;
omega(1) = 0;
delta_t = 0.01;
time = 10;
steps = time/delta_t;
%Analytic answer
t = [0:delta_t:time];
theta_a = theta(1)*cos(sqrt(g/l)*t);
%Euler solver
for i = 1:steps
E_E(i)= (1/2)*m*l^2*omega(i)^2+m*g*l*(1-cos(theta(i))); %#ok<*SAGROW>
omega(i+1) = omega(i)-(g/l)*theta(i)*delta_t;
theta(i+1) = theta(i)+omega(i)*delta_t;
end
%Last step of total Energy
E_E (i+1)= (1/2)*m*l^2*omega(i+1)^2+m*g*l*(1-cos(theta(i+1)));
figure(1)
hold on;
% plot(t,theta)
plot(t,theta_a,'.')
figure(2)
hold on;
plot(t,E_E)
%Euler-Cromer
for i = 1:steps
E_EC (i)= (1/2)*m*l^2*omega(i)^2+m*g*l*(1-cos(theta(i))); %#ok<*SAGROW>
omega(i+1) = omega(i)-(g/l)*theta(i)*delta_t;
theta(i+1) = theta(i)+omega(i+1)*delta_t;
end
%Last step of total Energy
E_EC (i+1)= (1/2)*m*l^2*omega(i+1)^2+m*g*l*(1-cos(theta(i+1)));
%
% figure(1)
% hold on;
% plot(t,theta)
% figure(2)
% plot(t,E_EC)
for i = 1:steps
k1 = -(g/l)*theta(i);
k2 = -(g/l)*(theta(i)+delta_t*0.5*k1);
k3 = -(g/l)*(theta(i)+delta_t*0.5*k2);
k4 = -(g/l)*(theta(i)+delta_t*k3);
omega(i+1) = omega(i)+((k1+2*k2+2*k3+k4)/6)*delta_t;
k1 = omega(i);
k2 = omega(i)+delta_t*0.5*k1;
k3 = omega(i)+delta_t*0.5*k2;
k4 = omega(i)+delta_t*k3;
theta(i+1) = theta(i)+((k1+2*k2+2*k3+k4)/6)*delta_t;
end
figure(1)
hold on;
plot(t,theta,'x')
采纳的回答
darova
2019-9-30
0 个投票
Please use button fro code inserting
2 个评论
Truls Kjøsnes Olsen
2019-9-30
darova
2019-9-30
I think each step of RK4 should be connected
for i = 1:steps-1
k1o = -(g/l)*theta(i);
k1t = omega(i);
k2o = -(g/l)*(theta(i)+delta_t*0.5*k1t);
k2t = omega(i) + delta_t*0.5*k2o;
k3o = -(g/l)*(theta(i)+delta_t*0.5*k2t);
k3t = omega(i) + delta_t*0.5*k3o;
k4o = -(g/l)*(theta(i)+delta_t*k3t);
k4t = omega(i) + delta_t*k4o;
omega(i+1) = omega(i)+((k1o+2*k2o+2*k3o+k4o)/6)*delta_t;
theta(i+1) = theta(i)+((k1t+2*k2t+2*k3t+k4t)/6)*delta_t;
end
Dimensions must agree
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
