Runge Kutta method for coupled oscillator system.

12 views (last 30 days)
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));

Accepted Answer

Alan Stevens
Alan Stevens on 10 Nov 2020
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 Comments
Alan Stevens
Alan Stevens on 12 Nov 2020
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Automotive in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!