Code covered by the BSD License

# Dynamical Systems Toolbox

13 Jul 2011

Bifurcation analysis of dynamical systems. Integration of AUTO bifurcation software into MATLAB.

run7d.m
```% Numerical simulation of the shimmy equations given by somieski

% Initializing the problem parameters
parameterset

% Initializing the option of the numerical integration
options=odeset('RelTol',1e-6,'AbsTol',1e-6,'events','on');

tspan=[0 10];
parameters(13)=30;

inicon=[   0.047775109049949 -42.220102984697597   0.001810029349389  -1.414095995444490   0.020328927171386];
[T,Y,TE,YE,IE]=ode45('odefile1',tspan,inicon,options,parameters);
plot(Y(:,1),Y(:,2))

data=[];
vel=[];
for v=248:0.1:250
vel=[];
parameters(13)=v;
parameters(13)
inicon=[Y(end,:)];
tspan=[0 1];
options=odeset('RelTol',1e-6,'AbsTol',1e-6,'events','off');
[T,Y]=ode45('odefile1',tspan,inicon,options,parameters);
inicon=[Y(end,:)];
tspan=[0 1];
options=odeset('RelTol',1e-6,'AbsTol',1e-6,'events','on');
[T,Y,TE,YE,IE]=ode45('odefile1',tspan,inicon,options,parameters);
vel(1:length(YE(:,1)),1)=v;
data=[data;vel IE YE];
end

% Testing Fy using a beam deflection equation
if 0
% Fy Vs deflection under pure bending moment at the free end of the beam
x=[0.00001:0.001:1];
c=1;
y=((c^2*x.^2-1).*sqrt(-1./(c^2*x.^2-1))+1)/c;
plot(x,-y,'r.')

c=[0.00001:0.001:1];
L=1;
y=((L^2*c.^2-1).*sqrt(-1./(L^2*c.^2-1))+1)./c;
plot(y,c,'r.')
% Fy vs deflection from the models of the paper by Osman
c=[0.0001:0.001:10];
y1=c.*(1-cos(0.1./c));
y2=c.*(1-sqrt((0.1./c).^2));
y3=(c/2).*((0.1./c).^2);
plot(c,y1,'r.')
hold on

plot(c,y2,'b.')
plot(c,y3,'k.')

P=1;
L=1;
r=0.5;
EI=1;
x=[0:0.001:L+r];
y=-P*(L+r).*x.^2+P.*(x.^3)/6+(P*L^2)/(2*EI)*r;
plot(x,y,'.')
alpha=[0:0.001:1];
Fyvec=[];
for i=1:length(alpha)
if abs(alpha(i)) <= 0.09
Fy=20*alpha(i);
else
Fy=20*5*pi/180*sign(alpha(i));
end
Fyvec=[Fyvec Fy];
end
plot(alpha*180/pi,Fyvec)
Fyest=20*0.1*tan(0.3*alpha)/(0.03+(3*tan(0.3*alpha))^2);
Fyest=20*0.1*tan(alpha)./(0.1+(1.3*tan(alpha)).^2);
Fyest=3*atan(7*tan(alpha)).*cos(0.95*atan(7*tan(alpha)));
end
figure
plot(T,Y(:,1),'r')
hold on
plot(T,Y(:,3),'b')
plot(T,Y(:,5),'k')
hold on