필터 지우기
필터 지우기

can anyone help me to correct the matlab code?

조회 수: 3 (최근 30일)
prajith samraj
prajith samraj 2022년 4월 17일
댓글: DGM 2022년 4월 18일
function vdpd_qo_syn
clear all
clear all
figure(1)
%Computational Length and step size
n=100; %total data legnth (n/h)
h=0.01;
%%%%%%%%%%%%%%%%%%%system 1%%%%%%%%%%%%%%%%
alpha= 0.01;
beta = 9.74;
omegaf= 0.3;
gamma= 0.2;
W1= 0.3;
ff1 =0.4999;
%% Initial condition
x10 = 0.11; x20 = 0.2;
%%%%%%%%%%%%%%%system 2%%%%%%%%%%%%%%%%
e = 0.22;
a= 0.956;
b = 0.149;
ff2 =.59;
W2= 0.3;
%% Initial condition
x30 = 0.11; x40 = 0.2;
%%%%%%%%%%%%%%%%%%coupling %%%%%%%%%%%%%%%%
cou=0.2;
%time step and initial condition
tspan = 0:h:n;
y0 = [x10; x20; x30; x40];
%[t,y] = ode45(@(t,x) f(t,x,a,b,beta,rho,ff1,W1,W2,ff2,cou),tspan,y0);
%[t,y] = ode45(@(t,x) f(t,x,alpha,omegaf,beta,gamma,ff1,ff2,W1,W2,epsilon,alpha1,beta1,cou),tspan,y0);
[t,y] = ode45(@(t,x) f(t,x,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou),tspan,y0);
%% transient%%%%%%%%%%%%%%%%remove%%%%%%%%%%%%%%%
k=n/h;
k1=(k/2:k);
size(k1);
x1=y(k1,1); x2=y(k1,2); x3=y(k1,3); x4=y(k1,4); ti=t(k1);
%plot the variable
%plot(x1,x2)
subplot(2,2,4);plot(t,y(:,1)); xlabel('t'); ylabel('x'); title('time series of qo oscillaor');
subplot(2,2,1); plot(x1,x2); xlabel('x'); ylabel('y'); title('Quintic oscillator ');
subplot(2,2,2); plot(x3,x4); xlabel('x^P'); ylabel('y^P'); title('VdpD oscillator');
subplot(2,2,3); plot(x1,x3); xlabel('x'); ylabel('y'); title('Synchronization');
figure(2)
subplot(2,1,1);plot(ti,x1,'r',ti,x3,'b'); xlabel('Time(sec)'); ylabel('x1,y1');
subplot(2,1,2);plot(ti,x1,'r',ti,x3,'b'); xlabel('Time(sec)'); ylabel('x1,y1');
%figure(2)
% plot(ti,x2,'r',ti,x4,'b'); xlabel('Time(sec)'); ylabel('x2,y2');
fprintf('Total length %d and the code taken transient after %d', k,k/2)
function dy = f(t,y,alpha,beta,omegaf,gamma,ff1,W1,e,a,b,ff2,W2,cou)
x1 = y(1); x2 = y(2); x3 = y(3); x4 = y(4);
%e1 = x3-x1;
%e2= x4-x2;
%sys-I
dx1=x2+cou*(x3-x1);
dx2=-alpha*x2-omegaf*x1-beta*x1.^3- gamma*x1.^5+ff1*cos(W1*t);
%Sys II
u1 = -x4+x2-x3+x1;
u2 = -(omegaf-alpha +beta*x1.^2 +gamma*x1.^4)*x1 -b*x2.^3-(e+alpha)*x3 +e*x2.^2*x4 -f2*sin(W2*t)-f1*cos(W1*t);
dx3=x4+u1;
dx4=e*(1-x3.^2)*x4-a*x3+b*x3.^2+ff2*cos(W2*t)+u2;
dy = [dx1; dx2; dx3; dx4];

채택된 답변

DGM
DGM 2022년 4월 17일
Just a guess, but f1 and f2 are undefined.
If i change them to ff1 and ff2, the code runs:
u2 = -(omegaf-alpha +beta*x1.^2 +gamma*x1.^4)*x1 -b*x2.^3-(e+alpha)*x3 +e*x2.^2*x4 -ff2*sin(W2*t)-ff1*cos(W1*t);
You'll have to decide if that substitution makes sense. I only made that change based on the context. I have no insight into what this code is doing.
  댓글 수: 1
DGM
DGM 2022년 4월 18일
What's wrong with the coupling parameter that you had defined originally?

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Cartesian Coordinate System Conversion에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by