Runge Kutta 4 method to solve second order ODE
105 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Please help. I have been stuck at it for a while:
I am trying to solve a second order differential equation where U_dot= V and V_dot = d*U-c*U^3-b*V+a*sin(w*t)
Now I made the code (analytical part is meant for double checking myself, ignore it). My code gives an error saying "index exceeds array bounds". The loop refuses to accept anything other than for i = 1:1. How can I run this loop with Runge Kutta calculations?
Thank you!
w = 1.3;
dt = 2*pi/(w*100);
a= 0.25;
b=0.1;
c=1;
d = 1;
x0=1;
y0=0;
t=0:dt:5000;
transient = 4250;
tran_str= 300;
% %analytical
% w0=sqrt(-d);
% A=sqrt(x0^2+(y0/w0)^2);
% phi = atan(x0*w0/y0);
% xan = A*sin(w0*t+phi);
% yan = A*w0*cos(w0*t+phi);
%RK4
f1 = @(t,x,y) y;
f2 = @(t,x,y) d*x-c*x^3-b*y+a*sin(w*t);
h=dt
for i=1:length(t-1)
t(i+1) = t(i)+i*h;
k1x = f1(t(i),x0(i),y0(i));
k2x = f2(t(i)+0.5*h,x0(i)+0.5*k1y*h, y0(i)+0.5*k1y*h);
k3x = f2(t(i)+0.5*h,x0(i)+0.5*k2y*h, y0(i)+0.5*k2y*h);
k4x = f2(t(i)+0.5*h,x0(i)+0.5*k3y*h, y0(i)+0.5*k3y*h);
k1y = f2(t(i),x0(i),y0(i));
k2y = f2(t(i)+0.5*h,x0(i)+0.5*k1y*h, y0(i)+0.5*k1y*h);
k3y = f2(t(i)+0.5*h,x0(i)+0.5*k2y*h, y0(i)+0.5*k2y*h);
k4y = f2(t(i)+0.5*h,x0(i)+0.5*k3y*h, y0(i)+0.5*k3y*h);
y(i+1) = y0(i)+1/6*(k1y+2*k2y+2*k3y+k4y);
end
3 commentaires
Réponses (1)
Alan Stevens
le 15 Nov 2020
Modifié(e) : Alan Stevens
le 15 Nov 2020
Your integration loop is mightily scrambled. It should be like this
x(1) = x0;
y(1) = y0;
for i=1:length(t)-1
t(i+1) = i*h;
k1x = f1(t(i),x(i),y(i));
k1y = f2(t(i),x(i),y(i));
k2x = f1(t(i)+0.5*h,x(i)+0.5*k1x*h, y(i)+0.5*k1y*h);
k2y = f2(t(i)+0.5*h,x(i)+0.5*k1x*h, y(i)+0.5*k1y*h);
k3x = f1(t(i)+0.5*h,x(i)+0.5*k2x*h, y(i)+0.5*k2y*h);
k3y = f2(t(i)+0.5*h,x(i)+0.5*k2x*h, y(i)+0.5*k2y*h);
k4x = f1(t(i)+h,x(i)+k3x*h, y(i)+k3y*h);
k4y = f2(t(i)+h,x(i)+k3x*h, y(i)+k3y*h);
x(i+1) = x(i)+1/6*(k1x+2*k2x+2*k3x+k4x)*h;
y(i+1) = y(i)+1/6*(k1y+2*k2y+2*k3y+k4y)*h;
end
(However, I don't think your analytical solution is the solution to your equations.)
1 commentaire
Bruce Taylor
le 22 Oct 2023
This is a beautiful solution. I modified it to analyze a series R-L-C circuit and compared the result to the exact solution...perfect match.
Bruce Taylor
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!