Frequency vs Amplitude using ODE45 of system of equations?
1 view (last 30 days)
Show older comments
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 Comments
Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!