can anyone help me to correct the matlab code?

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

Accepted Answer

DGM
DGM on 17 Apr 2022
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 Comment
DGM
DGM on 18 Apr 2022
What's wrong with the coupling parameter that you had defined originally?

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!