Runge Kutta method for coupled oscillator system.
12 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Alan Stevens
2020-11-10
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 个评论
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!