# Plot of three coupled oscillator.

4 views (last 30 days)
Haya Ali on 19 Dec 2022
Commented: Haya Ali on 20 Dec 2022
I am trying to built three coupled oscillator just as two coupled oscillator but i don't know what mistake I am doing. Please help me. The code for both the cases are given below
%%Two coupled
clear all; close all; clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%value of constants
a1=0.1;a2=0.2;
omega1=5;omega2=4;
G=0.01;C12=0.001;C21=0.002;
dt=0.01; %step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1(1)=0.5;
y1(1)=0.5;
x2(1)=0.5;
y2(1)=0.5;
for i=2:1000
x1(i)=x1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*x1(i-1)-omega1*y1(i-1)+G*C12*(x2(i-1)-x1(i-1)))*dt;
y1(i)=y1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*y1(i-1)+omega1*x1(i-1)+G*C12*(y2(i-1)-y1(i-1)))*dt;
x2(i)=x2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*x2(i-1)-omega2*y2(i-1)+G*C21*(x1(i-1)-x2(i-1)))*dt;
y2(i)=y2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*y2(i-1)+omega2*x2(i-1)+G*C21*(y1(i-1)-y2(i-1)))*dt;
end
figure
hold on
plot(x1,'r')
plot(x2,'g')
%% Three Coupled
close all; clear all; clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%value of constants
a1=0.1;a2=0.2;a3=0.5;
omega1=5;omega2=4;omega3=3;
G=0.01;C12=0.001;C13=0.006;
C21=0.003;C23=0.005;
C31=0.002;C32=0.004;
dt=0.01; %step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1(1)=0.5;
y1(1)=0.5;
x2(1)=1;
y2(1)=1;
x3(1)=1.5;
y3(1)=1.5;
for i=2:2000
x1(i)=x1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*x1(i-1)-omega1*y1(i-1)+G*C12*(x2(i-1)-x1(i-1)))+G*C13*(x3(i-1)-x1(i-1))*dt;
y1(i)=y1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*y1(i-1)+omega1*x1(i-1)+G*C12*(y2(i-1)-y1(i-1)))+G*C13*(y3(i-1)-y1(i-1))*dt;
x2(i)=x2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*x2(i-1)-omega2*y2(i-1)+G*C21*(x1(i-1)-x2(i-1)))+G*C23*(x3(i-1)-x2(i-1))*dt;
y2(i)=y2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*y2(i-1)+omega2*x2(i-1)+G*C21*(y1(i-1)-y2(i-1)))+G*C23*(y3(i-1)-y2(i-1))*dt;
x3(i)=x3(i-1)+((a3-x3(i-1)^2-y3(i-1)^2)*x3(i-1)-omega3*y3(i-1)+G*C31*(x1(i-1)-x3(i-1)))+G*C32*(x2(i-1)-x3(i-1))*dt;
y3(i)=y3(i-1)+((a3-x3(i-1)^2-y3(i-1)^2)*y3(i-1)+omega3*x3(i-1)+G*C31*(y1(i-1)-y3(i-1)))+G*C32*(y2(i-1)-y3(i-1))*dt;
end
figure
hold on
plot(x1,'r')
plot(x2,'g')
plot(x3,'m')

Sam Chak on 19 Dec 2022
I don't know what oscillators are since you didn't provide the math equations to compare with. However, on the programming side, it appears that the location of parenthesis caused the problem. It is possible to randomly shuffle the locations of the brackets by some algorithms until the system responses are stable.
Check if these are the expected responses.
% value of constants
a1 = 0.1;
a2 = 0.2;
a3 = 0.5;
omega1 = 5;
omega2 = 4;
omega3 = 3;
G = 0.01;
C12 = 0.001;
C13 = 0.006;
C21 = 0.003;
C23 = 0.005;
C31 = 0.002;
C32 = 0.004;
dt = 0.01; % step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1(1) = 0.5;
y1(1) = 0.5;
x2(1) = 1;
y2(1) = 1;
x3(1) = 1.5;
y3(1) = 1.5;
for i = 2:1000
% x1(i) = x1(i-1) + ((a1 - x1(i-1)^2 - y1(i-1)^2)*x1(i-1) - omega1*y1(i-1) + G*C12*(x2(i-1) - x1(i-1))) + G*C13*(x3(i-1) - x1(i-1))*dt;
% y1(i) = y1(i-1) + ((a1 - x1(i-1)^2 - y1(i-1)^2)*y1(i-1) + omega1*x1(i-1) + G*C12*(y2(i-1) - y1(i-1))) + G*C13*(y3(i-1) - y1(i-1))*dt;
% x2(i) = x2(i-1) + ((a2 - x2(i-1)^2 - y2(i-1)^2)*x2(i-1) - omega2*y2(i-1) + G*C21*(x1(i-1) - x2(i-1))) + G*C23*(x3(i-1) - x2(i-1))*dt;
% y2(i) = y2(i-1) + ((a2 - x2(i-1)^2 - y2(i-1)^2)*y2(i-1) + omega2*x2(i-1) + G*C21*(y1(i-1) - y2(i-1))) + G*C23*(y3(i-1) - y2(i-1))*dt;
% x3(i) = x3(i-1) + ((a3 - x3(i-1)^2 - y3(i-1)^2)*x3(i-1) - omega3*y3(i-1) + G*C31*(x1(i-1) - x3(i-1))) + G*C32*(x2(i-1) - x3(i-1))*dt;
% y3(i) = y3(i-1) + ((a3 - x3(i-1)^2 - y3(i-1)^2)*y3(i-1) + omega3*x3(i-1) + G*C31*(y1(i-1) - y3(i-1))) + G*C32*(y2(i-1) - y3(i-1))*dt;
x1(i) = x1(i-1) + ( ( a1 - x1(i-1)^2 - y1(i-1)^2 )*x1(i-1) - omega1*y1(i-1) + G*C12*( x2(i-1) - x1(i-1) ) + G*C13*( x3(i-1) - x1(i-1) ) )*dt;
y1(i) = y1(i-1) + ( ( a1 - x1(i-1)^2 - y1(i-1)^2 )*y1(i-1) + omega1*x1(i-1) + G*C12*( y2(i-1) - y1(i-1) ) + G*C13*( y3(i-1) - y1(i-1) ) )*dt;
x2(i) = x2(i-1) + ( ( a2 - x2(i-1)^2 - y2(i-1)^2 )*x2(i-1) - omega2*y2(i-1) + G*C21*( x1(i-1) - x2(i-1) ) + G*C23*( x3(i-1) - x2(i-1) ) )*dt;
y2(i) = y2(i-1) + ( ( a2 - x2(i-1)^2 - y2(i-1)^2 )*y2(i-1) + omega2*x2(i-1) + G*C21*( y1(i-1) - y2(i-1) ) + G*C23*( y3(i-1) - y2(i-1) ) )*dt;
x3(i) = x3(i-1) + ( ( a3 - x3(i-1)^2 - y3(i-1)^2 )*x3(i-1) - omega3*y3(i-1) + G*C31*( x1(i-1) - x3(i-1) ) + G*C32*( x2(i-1) - x3(i-1) ) )*dt;
y3(i) = y3(i-1) + ( ( a3 - x3(i-1)^2 - y3(i-1)^2 )*y3(i-1) + omega3*x3(i-1) + G*C31*( y1(i-1) - y3(i-1) ) + G*C32*( y2(i-1) - y3(i-1) ) )*dt;
end
figure
hold on
plot(x1)
plot(x2)
plot(x3)
grid on, legend('x', 'y', 'z')
Sam Chak on 19 Dec 2022
This is another version in 3-D.
% value of constants
a1 = 0.1;
a2 = 0.2;
a3 = 0.5;
omega1 = 5;
omega2 = 4;
omega3 = 3;
G = 0.01;
C12 = 0.001;
C13 = 0.006;
C21 = 0.003;
C23 = 0.005;
C31 = 0.002;
C32 = 0.004;
dt = 0.01; % step size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1(1) = 0.5;
y1(1) = 0.5;
x2(1) = 1;
y2(1) = 1;
x3(1) = 1.5;
y3(1) = 1.5;
for i = 2:6000
x1(i) = x1(i-1) + ( ( a1 - x1(i-1)^2 - y1(i-1)^2 )*x1(i-1) - omega1*y1(i-1) + G*C12*( x2(i-1) - x1(i-1) ) + G*C13*( x3(i-1) - x1(i-1) ) )*dt;
y1(i) = y1(i-1) + ( ( a1 - x1(i-1)^2 - y1(i-1)^2 )*y1(i-1) + omega1*x1(i-1) + G*C12*( y2(i-1) - y1(i-1) ) + G*C13*( y3(i-1) - y1(i-1) ) )*dt;
x2(i) = x2(i-1) + ( ( a2 - x2(i-1)^2 - y2(i-1)^2 )*x2(i-1) - omega2*y2(i-1) + G*C21*( x1(i-1) - x2(i-1) ) + G*C23*( x3(i-1) - x2(i-1) ) )*dt;
y2(i) = y2(i-1) + ( ( a2 - x2(i-1)^2 - y2(i-1)^2 )*y2(i-1) + omega2*x2(i-1) + G*C21*( y1(i-1) - y2(i-1) ) + G*C23*( y3(i-1) - y2(i-1) ) )*dt;
x3(i) = x3(i-1) + ( ( a3 - x3(i-1)^2 - y3(i-1)^2 )*x3(i-1) - omega3*y3(i-1) + G*C31*( x1(i-1) - x3(i-1) ) + G*C32*( x2(i-1) - x3(i-1) ) )*dt;
y3(i) = y3(i-1) + ( ( a3 - x3(i-1)^2 - y3(i-1)^2 )*y3(i-1) + omega3*x3(i-1) + G*C31*( y1(i-1) - y3(i-1) ) + G*C32*( y2(i-1) - y3(i-1) ) )*dt;
end
figure
hold on
plot3(x1, x2, x3)
view(45, 30)
Haya Ali on 20 Dec 2022
Thank you for the efforts Sam :)