Runge Kutta 4 method to solve second order ODE

17 次查看（过去 30 天）
Alina Li
Alina Li 2020-11-14
评论： Bruce Taylor 2023-10-22
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 个评论
显示 1更早的评论
Alina Li
Alina Li 2020-11-14
Thank you for the question, it actually made me realize that I should have included it. I edited the code above. Unfrotunately, the same error persists.
Haya M
Haya M 2020-11-15
What do you mean by k1y in k2x and check the rest as well.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Alan Stevens
Alan Stevens 2020-11-15
编辑：Alan Stevens 2020-11-15
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 个评论
Bruce Taylor
Bruce Taylor 2023-10-22
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by