can anyone help me to correct the matlab code?
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Akzeptierte Antwort
DGM
am 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 Kommentar
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Filter Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!