Plot of three coupled oscillator.

3 次查看（过去 30 天）
Haya Ali
Haya Ali 2022-12-19
评论： Haya Ali 2022-12-20
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
Sam Chak 2022-12-19
Ran in:
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')
  4 个评论
显示 2更早的评论
Sam Chak
Sam Chak 2022-12-19
Ran in:
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
Haya Ali 2022-12-20
Thank you for the efforts Sam :)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by