Frequency vs Amplitude using ODE45 of system of equations?
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have a system of equations and I want to plot the frequency vs amplitude curves for the system.
this is my function
function xdot= multifunction(t,x)
global alfa Meo v eta zeta P vt psi
xdot = zeros(8,1);
xdot(1)=x(2);
xdot(2)= - ((1/2)*pi^4 *(zeta+1)*x(1)) - Meo *x(2) - (1/8)* eta * (-x(6)*x(1)^2 + 4*x(6)*x(1)*x(5) + x(6)*x(3)*(x(3)+2*x(7)) ...
+ 2*x(8)*(-x(1)*x(3) + x(3)*x(5) + 2*x(1)*x(7)) + 2*x(4)*(x(1)*(2*x(3)-x(7))+ x(5)*(x(3)+x(7))) + x(2)*(3*x(1)^2 - 2*x(1)*x(5) ...
+ 2*(x(3)^2 + x(5)^2 - x(3)*x(7) + x(7)^2))) ...
- (1/8)*pi^4*x(1)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) ...
- (1/2)*alfa*v*pi*(-2*x(2)*x(3) + 2*x(6)*x(3) + x(4)*(x(1) - 3*x(5)) + 3*x(8)*x(5) - 4*x(6)*x(7)) + ((1/2)*pi^2 * x(1)*(alfa^2 * v^2 - P)) ...
- (psi*vt)*((4/15)*(5*x(3) + 2*x(7)));
xdot(3)=x(4);
xdot(4)= - ((8)*pi^4 *(zeta+1)*x(3)) - Meo *x(4) - (1/8)* eta * (x(8)*(-x(1)^2 + 2*x(1)*x(5) + x(5)^2 + 4*x(3)*x(7)) ...
+2*(x(2)*x(1)*(2*x(3) - x(7)) + x(4)*x(1)*(x(3) + x(7)) + x(2)*x(5)*(x(3) + x(7))+ x(6)*x(5)*(2*x(3) + x(7))) ...
+ x(4)*(2*x(1)^2 + 2*x(1)*x(5) + 3*x(3)^2 + 2*x(5)^2 + 2*x(7)^2)) - (1/2)*pi^4*x(3)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) ...
- (1/2)*alfa*v*pi*( x(7)*x(1) + 2*x(8)*x(3) + x(2)*(x(1) - 3*x(5)) - 4*x(4)*x(7)) + (2*pi^2 * x(3)*(alfa^2 * v^2 - P)) ...
- (psi*vt)*((4/15)*(5*x(1) - 9*x(5)));
xdot(5)=x(6);
xdot(6)= - ((81/2)*pi^4 *(zeta+1)*x(5)) - Meo *x(6) - (1/8)* eta * (2*x(6)*x(1)^2 + 3*x(6)*x(5)^2 + 2*x(8)*x(3)*(x(1)+x(5)) ...
+ 2*x(4)*x(3)*(x(1)+2*x(5)) + 4*x(8)*x(5)*x(7) +2*x(4)*x(7)*(x(1)+x(5)) + 2*x(6)*(x(3)^2 + x(3)*x(7)+x(7)^2) + x(2)*(-x(1)^2 + 4*x(1)*x(5) + x(3)*(x(3)+2*x(7)))) ...
- (9/8)*pi^4*x(5)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) - (1/2)*alfa*v*pi*( x(4)*x(1) + x(8)*x(1) + 2*x(2)*(x(3) - 2*x(7)));
+ ((9/2)*pi^2 * x(5)*(alfa^2 * v^2 - P)) - (psi*vt)*((12/35)*(7*x(3) - 10*x(7)));
xdot(7)=x(8);
xdot(8)= - ((128)*pi^4 *(zeta+1)*x(7)) - Meo *x(8) - (1/8)* eta * ((2*(x(2)*x(3))*(-x(1) + x(5)) + x(6)*x(3)*(x(1)+x(5)) + x(8)*(x(1)^2 + x(3)^2 +x(5)^2)) ...
+ 4*x(7)*(x(2)*x(1) + x(6)*x(5)) + 3* x(8)*x(7)^2 + x(4)*(-x(1)^2 + 2*x(1)*x(5) + x(5)^2 + 4*x(3)*x(7))) ...
- 2*pi^4*x(7)*(x(1)^2 + 4*x(3)^2 + 9*x(5)^2 + 16*x(7)^2) - (1/2)*alfa*v*pi*(x(6)*x(1) + 2*x(4)*x(3) + 3*x(2)*x(5));
+ (8*pi^2 * x(7)*(alfa^2 * v^2 - P)) - (psi*vt)*((8/105)*(7*x(1) + 45*x(5)));
end
all material properties and functions are specified and coded in another script. I only need the ODE45 code to execute and get the frequency vs amplitude curves
The script I wrote is to plot the Phase portrait and the wave form.
clc;clear all;close all;format long;
global alfa Meo v vt eta zeta P psi
E=1.1e12;
row=1.3e3;
G=280e9;
L=60e-9;
h=1e-9;
w= 10*h;
gap= 5*h;
I=(1/12)*w*(h^3);
Ar=w*L;
v= 2;
vt = 1;
Meo1=200e-8;
Po = 1e5;
Pa = 1e7;
lambda=64e-9;
lam = (lambda*Po)/(Pa*gap);
cpp = (Meo1*(L^2)*(w^3))/((gap^3)*sqrt(E*I*Ar*row));
Pr=100;
ssf = h;
eta1 = 2e7;
r = sqrt(1/12)*h;
zeta = ((L^2)*G*ssf)/(E*(r^2));
alfa = r/L;
P = (Pr*L^2)/(E*(r^2));
eta = (eta1*(L^2)*sqrt(row/E)/(Ar*row));
Meo = (cpp*(L^2)*sqrt(row/E)/(Ar*row*r));
psi = L*sqrt((row*Ar)/(E*I));
%% ODE Solver
t0=0;
Tt=2*pi/81.8
Tf=500*Tt
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,x]=ode45(@multifunction,0:Tt/100:Tf,[.01 0 0 0 0 0 0 0],options);
n=length(t);
as=fix(0.98*n);
t=t(as:n);
x1=x(as:n,1);
x2=x(as:n,2);
%% Results
figure(1)
set(gcf,'Color',[1,1,1])
subplot(1,2,1);hold on;
plot(t,x1,'k','linewidth',0.1)
grid on
title (' Wave Form','FontSize',12)
xlabel('Time','FontSize',12)
ylabel('$q1$', 'interpreter', 'latex','FontSize',20)
%%
subplot(1,2,2);hold on;
plot(x1,x2,'k','linewidth',0.1)
grid on
title ('Phase Portrait','FontSize',12)
xlabel('$q1$', 'interpreter', 'latex','FontSize',20)
ylabel('$\dot q1$', 'interpreter', 'latex','FontSize',20)
Please help!
0 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!