can anyone help me to correct the matlab code?
2 views (last 30 days)
Show older comments
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];
1 Comment
Accepted Answer
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
More Answers (0)
See Also
Categories
Find more on Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!