No BSD License  

Highlights from
Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

image thumbnail

Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

by

Howard Wilson

 

14 Oct 2002 (Updated )

Companion Software (amamhlib)

toprun
function toprun      
% Example: toprun 
% ~~~~~~~~~~~~~~~
%
% Example that analyzes the response of a
% spinning conical top.
%
% User m functions required:
%    topde, cubrange, inputv

disp(' ');
disp(['*** Dynamics of a Homogeneous ', ...
      'Conical Top ***']); disp(' ');
disp(['Input the gravity constant and the ', ...
      'body weight (try 32.2,5)']);
[grav,wt]=inputv('? '); 
mass=wt/grav; tmp=zeros(3,1); 
disp(' '); 
disp(['Input the height and base radius ', ...
      '(try 1,.5)']);
[ht,rb]=inputv('? '); len=.75*ht;
jtrans=3*mass/20*(rb*rb+4*ht*ht); 
jaxial=3*mass*rb*rb/10;
disp(' ');
disp(['Input a vector along the initial ', ...
      'axis direction (try 0,1,0)']);
[tmp(1),tmp(2),tmp(3)]=inputv('? '); 
e3=tmp(:)/norm(tmp); r0=len*e3;
disp(' ');
disp(['Input the initial angular velocity ', ...
      '(try 0,10,2)']);
[tmp(1),tmp(2),tmp(3)]=inputv('? '); omega0=tmp; 
omegax=e3'*omega0(:); rdot0=cross(omega0,r0); 
z0=[r0(:);rdot0(:)]; uz=[0;0;1]; 
c1=wt*len^2/jtrans; c2=omegax*jaxial/jtrans;
disp(' '); 
disp(['Input tfinal,and the integration ', ...
      'tolerance (try 4.2, 1e-8)']);
[tfinl,tol]=inputv('? '); disp(' ');
fprintf( ...
  'Please wait for solution of equations.\n');

% Integrate the equations of motion
odeoptn=odeset('RelTol',tol);
[tout,zout]=ode45(@topde,[0,tfinl],z0,...
                  odeoptn,uz,c1,c2);
t=tout; x=zout(:,1); y=zout(:,2); z=zout(:,3);
vx=zout(:,4); vy=zout(:,5); vz=zout(:,6); 

% Compute total energy and angular momentum 
c3=jtrans/(len*len); taxial=jaxial/2*omegax^2;
r=zout(:,1:3)'; v=zout(:,4:6)';
etotal=(wt*r(3,:)+taxial+c3/2*sum(v.*v))'; 
h=(jaxial*omegax/len*r+c3*cross(r,v))';

% Plot the path of the gravity center
close; axis('equal'); 
axis(cubrange([x(:),y(:),z(:)])); plot3(x,y,z);
title('Path of the Top Gravity Center');
xlabel('x axis'); ylabel('y axis'); 
zlabel('z axis'); grid on; figure(gcf); 
disp(' '); disp(...
'Press [Enter] to plot error measures'), pause 
% print -deps toppath 
n=2:length(t);

% Compute energy and angular momentum error 
% quantities and plot results
et=etotal(1); enrger=abs(100*(etotal(n)-et)/et);
hzs=abs(h(1,3)); 
angmzer=abs(100*(h(n,3)-hzs)/hzs);
vec=[enrger(:);angmzer(:)]; 
minv=min(vec); maxv=max(vec);

clf; 
semilogy(t(n),enrger,'-r',t(n),angmzer,':m');
axis('normal'); xlabel('time'); 
ylabel('percent variation');
title(['Percent Variation in Total Energy ', ...
       'and z-axis Angular Momentum']);
legend(' Energy      (Upper Curve)', ...
       ' Ang. Mom. (Bottom Curve)',4);
figure(gcf), pause
% print -deps topvar

disp(' '), disp('All Done')

%=============================================

function zdot=topde(t,z,uz,c1,c2)
%
% zdot=topde(t,z,uz,c1,c2)
% ~~~~~~~~~~~~~~~
%
% This function defines the equation of motion
% for a symmetrical top. The vector z equals 
% [r(:);v(:)] which contains the Cartesian 
% components of the gravity center radius and 
% its velocity.
%
% t    - the time variable
% z    - the vector [x; y; z; vx; vy; vz]
% uz   - the vector [0;0;1]
% c1   - wt*len^2/jtrans
% c2   - omegax*jaxial/jtrans 
%
% zdot - the time derivative of z
%
% User m functions called:  none
%----------------------------------------------

z=z(:); r=z(1:3); len=norm(r); ur=r/len;

% Make certain the input velocity is 
% perpendicular to r
v=z(4:6); v=v-(ur'*v)*ur;
vdot=-c1*(uz-ur*ur(3))+c2*cross(ur,v)- ...
     ((v'*v)/len)*ur;
zdot=[v;vdot];

%=============================================

% function varargout=inputv(prompt)
% See Appendix B

% =============================================

% function range=cubrange(xyz,ovrsiz)
% See Appendix B

Contact us