Runge Kutta 4 method to solve second order ODE
109 views (last 30 days)
Show older comments
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 Comments
Answers (1)
Alan Stevens
on 15 Nov 2020
Edited: Alan Stevens
on 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 Comment
Bruce Taylor
on 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
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!