Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- r= 3.2
- s= 0.03
- f= 0.8
- %% Find stationary states
- syms xi
- eqn = -(r^2 * (xi-1)*((xi/2)-1)+1+(s/(1-xi)))*xi + f ;
- %eqn2= -(r^2 * (xi-1)*((xi/2)-1)*(1-xi)+(1-xi)+(s))*xi+(f*(1-xi));
- f1 = matlabFunction(expand(eqn))
- root1=fzero(f1,[-4, 0.999999999])
- root2=fzero(f1,[1.0000001, 4])
- %% Solve ODE
- omega = 0.01;
- T = pi/omega;
- t_start = 0;
- t_end = T;
- time_span = t_start:0.01:t_end;
- xi0 = [0 0];
- [t,xi] = ode45(@odefun, time_span,xi0);
- plot(t,xi);
- hold on
- plot (t, 3*sin(omega*t));
- ylim([0,3.5])
- legend('Displacement','Force','location','NorthEast')
- title('Displacement vs Time')
- xlabel('Time, tau [sec]')
- ylabel('Displacement, \xi')
- yyaxis right
- ylabel ('Force F/F0')
- ylim([0,3.5])
- function xidot = odefun(t,xi)
- xidot = -((3.2^2).*(xi-1).*(0.5*xi-1).*xi+xi+(xi.*0.03)./(1-xi))+3*sin(0.01*t);
- end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement