Plot of three coupled oscillator.

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')

Answer 1

[Sam Chak]

Hi @Haya Ali

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')

[Haya Ali]

Thanks a lot!

[Sam Chak]

Hi @Haya Ali, don't mention it as I'm also a learner. Sometimes, we tend to commit typing errors in very long equations that contain many brackets. I guess your was due to expanding to 3D from the copied 2D equations.

[Sam Chak]

Hi @Haya Ali

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]

Thank you for the efforts Sam :)