can anyone help me to correct the matlab code?

3 次查看(过去 30 天)
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.

更多回答(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