y= C1*exp(m1*t)+C2*exp(m2*t)+E1*sin(omega1*t)+E2*cos(omega1*t)+E3*sin(omega2*t)+E4*cos(omega2*t);
dy=m1*C1*exp(m1*t)+m2*C2*exp(m2*t)+E1*omega1*cos(omega1*t)-E2*omega1*sin(omega1*t)+E3*omega2*cos(omega2*t)-E4*omega2*sin(omega2*t);
Y1=1.0/alpha^2*(f1-omegaf*x0-beta*x0.^3-Ga*x0.^5)*(exp(-alpha*t)+alpha*t-1);
Y2= x1/alpha*(1-exp(-alpha*t)-(omegaf*x1+3*beta*x0.^2*x1)*(t^3/6- 2*alpha*t^4/24 +3*alpha^2*t^5/120));
Y3=5*Ga*x0.^4*x1*t.^3/6;
Y4=6*beta*x0*x1.^2 +20*Ga*x0.^3*x1.^2*(t.^4/24- 4*alpha*t.^5/120);
Y5=5*Ga*x0.^9-5*x0.^4*f1*(t^4/24 - alpha*Ga*t^5/120-alpha*t^5/120);
Y6=[f1*omega1.^2-f1*omegaf+omegaf^2*x0-3*beta*f1*x0.^2+ 4*omegaf^2*beta*x0.^3+3*beta*x0.^5+omegaf*Ga*x0.^5+5*omegaf*x0.^5+ 3*beta*Ga*x0.^7+5*beta*x0.^7+ 5*alpha*Ga*x0.^4*x1+5*alpha*x0.^4*x1]*t^4/24;
Y7=[2*alpha*f1*omegaf + alpha*f1*omega1.^2+ omegaf.^2*x1 -6*beta*x1.^3 -2*alpha*omegaf.^2*x0+5*alpha*omegaf*Ga*x0 + 6*alpha*beta*f1*x0.^2- 8*alpha*beta*omegaf*x0.^3+5*alpha*omegaf*x0.^5-2*Ga*x0.^5-11*alpha*beta*Ga*x0.^7 +5*alpha*beta*x0.^7+24*beta*omegaf*x0.^2*x1+70*Ga*omegaf*x0.^4*x1+27*beta^2*x0.^4*x1-10*alpha^2*Ga*x0.^4*x1-5*alpha.^2*x0.^4*x1 -18*f1*beta*x0*x1+108*beta*Ga*x0.^6*x1-60*Ga*x0.^2*x1.^3-60*Ga*f1*x0.^3*x1+ 85*Ga^2*x0.^8*x1]*t^5/120;
x = x0 + Y1 + Y2 - Y3 - Y4 + Y5 + Y6 + Y7;
% if(j>20)
fprintf('%20.12f %20.12f %20.12f \n',t,x,y);
% end
% incrementation of the variables and of time
x0=x;
y0=y;
z0=x;
t=t0+h;
t0=t;
end

1 commentaire

Torsten
Torsten le 18 Juin 2022
Which error ? We can't run your code since you didn't provide the necessary parameters.

Connectez-vous pour commenter.

Réponses (0)

Catégories

En savoir plus sur MATLAB dans Centre d'aide et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by