Runge Kutta method for coupled oscillator system.

Question

Heya :) 2020-11-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/643100-runge-kutta-method-for-coupled-oscillator-system

评论： Heya :) 2020-11-12

I am trying to solve these equations with the help of Runge Kutta Method. I am not sure whether I am writing the code correctly or we can use this method for coupled also getting this error (mentioned at the end of the code). Please help me to refine my code.

close all; clear all; clc;
%initializing x,y,t
t(1)=0;
x1(1)=1;
y1(1)=1;
x2(1)=2;
y2(1)=2;
%value of constants
a=0.1;
b=0.3;
omega=4;
Cnp=0.2;
G=1;
h=0.1;   %step size
t=0:h:50;
%ode
p=@(t,x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1);
q=@(t,x1,y1,x2,y2) (a-x1^2-y1^2)*y1+omega*x1+G*Cnp(y2-y1);
%loop
for i=1:(length(t)-1)
    k1=p(t(i),x1(i),y1(i),x2(i),y2(i));
    l1=q(t(i),x1(i),y1(i),x2(i),y2(i));
 
    
    k2=p(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
    l2=q(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
      
    k3=p(t(i)+h/2,(x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
    l3=q(t(i)+h/2,(x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
      
    k4=p(t(i)+h,(x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    l4=q(t(i)+h,(x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    
      x(i+1) = x(i) + h*(k1+2*k2+2*k3+k4)/6;
      y(i+1) = y(i) + h*(l1+2*l2+2*l3+l4)/6;
end
%draw
plot(t,x,'r')
xlabel('X','fontsize',14,'fontweight','bold')
ylabel('y','fontsize',14,'fontweight','bold')
set(gca,'Color','k')
****Error***
Array indices must be positive integers or logical values.
Error in coupled>@(t,x1,y1,x2,y2)(a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1) (line 17)
p=@(t,x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1);
Error in coupled (line 25)
    k2=p(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2020-11-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/643100-runge-kutta-method-for-coupled-oscillator-system#answer_540085

在 MATLAB Online 中打开

Do you mean like this

%initializing x,y,t
h=0.1;   %step size
t=0:h:50;
x1 = zeros(1,numel(t)); y1 = x1; x2 = x1; y2 = x1;
x1(1)=1;
y1(1)=1;
x2(1)=2;
y2(1)=2;
%value of constants
a=0.1;
b=0.3;
omega=4;
Cnp=0.2;
G=1;
%ode
p=@(x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp*(x2-x1);
q=@(x1,y1,x2,y2) (a-x1^2-y1^2)*y1+omega*x1+G*Cnp*(y2-y1);
%loop
for i=1:(length(t)-1)
    k1=p(x1(i),y1(i),x2(i),y2(i));
    l1=q(x1(i),y1(i),x2(i),y2(i));
 
    
    k2=p((x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
    l2=q((x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
      
    k3=p((x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
    l3=q((x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
      
    k4=p((x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    l4=q((x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    
      x1(i+1) = x1(i) + h*(k1+2*k2+2*k3+k4)/6;
      y1(i+1) = y1(i) + h*(l1+2*l2+2*l3+l4)/6;
      x2(i+1) = x2(i) + h*(k1+2*k2+2*k3+k4)/6;
      y2(i+1) = y2(i) + h*(l1+2*l2+2*l3+l4)/6;
end
%draw
plot(t,x1,'r',t,x2,'y')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('x1 & x2','fontsize',14,'fontweight','bold')
set(gca,'Color','k')
legend('x1','x2','TextColor','w')
figure
plot(t,y1,'r',t,y2,'y')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('y1 & y2','fontsize',14,'fontweight','bold')
set(gca,'Color','k')
legend('y1','y2','TextColor','w')

3 个评论
显示 1更早的评论隐藏 1更早的评论

Alan Stevens 2020-11-12

That is just setting aside storage space for those variables. It makes the loop more efficient as space doesn't then have to be found for their new elements every iteration of the loop.

Heya :) 2020-11-12