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)

deislner
function deislner  
%
% Example:  deislner
% ~~~~~~~~~~~~~~~~~~
% Solution error for simulation of cable
% motion using a second or a fourth order 
% implicit integrator.
%
% This program uses implicit second or fourth 
% order integrators to compute the dynamical 
% response of a cable which is suspended at 
% one end and is free at the other end. The 
% cable is given a uniform initial velocity. 
% A plot of the solution error is given for 
% two cases where approximate solutions are 
% generated using numerical integration rather 
% than modal response which is exact.
%
% User m functions required: 
%    mckde2i, mckde4i, cablemk, udfrevib, 
%    plterror

% Choose a model having twenty links of 
% equal length 

fprintf(...
'\nPlease wait: solution takes a while\n')
clear all; close
n=20; gravty=1.; n2=1+fix(n/2);
masses=ones(n,1)/n; lengths=ones(n,1)/n;

% First generate the exact solution by 
% modal superposition
[m,k]=cablemk(masses,lengths,gravty); 
c=zeros(size(m));
dsp=zeros(n,1); vel=ones(n,1); 
t0=0; tfin=50; ntim=126; h=(tfin-t0)/(ntim-1);

% Numbers of repetitions each solution is
% performed to get accurate cpu times for
% the chosen step sizes are shown below.
% Parameter jmr may need to be increased to
% give reliable cpu times on fast computers

jmr=500; 
j2=fix(jmr/50); J2=fix(jmr/25);
j4=fix(jmr/20); J4=fix(jmr/10);

% Loop through all solutions repeatedly to
% obtain more reliable timing values on fast
% computers
tic;  
for j=1:jmr;
  [tmr,xmr]=udfrevib(m,k,dsp,vel,t0,tfin,ntim);
end
tcpmr=toc/jmr;

% Second order implicit results
i2=10; h2=h/i2; tic;
for j=1:j2
  [t2,x2]=mckde2i(m,c,k,t0,dsp,vel,tfin,h2,i2);
end
tcp2=toc/j2; tr2=tcp2/tcpmr;

I2=5; H2=h/I2; tic;
for j=1:J2
  [T2,X2]=mckde2i(m,c,k,t0,dsp,vel,tfin,H2,I2);
end
Tcp2=toc/J2; Tr2=Tcp2/tcpmr;

% Fourth order implicit results
i4=2; h4=h/i4; tic; 
for j=1:j4
  [t4,x4]=mckde4i(m,c,k,t0,dsp,vel,tfin,h4,i4);
end
tcp4=toc/j4; tr4=tcp4/tcpmr;

I4=1; H4=h/I4; tic;
for j=1:J4
  [T4,X4]=mckde4i(m,c,k,t0,dsp,vel,tfin,H4,I4);
end
Tcp4=toc/J4; Tr4=Tcp4/tcpmr;

% Plot error measures for each solution
plterror(xmr,t2,h2,x2,T2,H2,X2,...
         t4,h4,x4,T4,H4,X4,tr2,Tr2,tr4,Tr4)

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

function [t,x,tcp] = ...
     mckde2i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
%
% [t,x,tcp]= ...
%    mckde2i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function uses a second order implicit 
% integrator % to solve the matrix differential 
% equation
%           m x'' + c x' + k x = forc(t)
% where m,c, and k are constant matrices and 
% forc is an externally defined function.
%
% Input:
% ------
% m,c,k   mass, damping and stiffness matrices
% t0      starting time
% x0,v0   initial displacement and velocity
% tmax    maximum time for solution evaluation
% h       integration stepsize
% incout  number of integration steps between 
%         successive values of output
% forc    externally defined time dependent 
%         forcing function. This parameter 
%         should be omitted if no forcing 
%         function is used.      
%
% Output:
% -------
% t       time vector going from t0 to tmax 
%         in steps of
% x       h*incout to yield a matrix of 
%         solution values such that row j 
%         is the solution vector at time t(j)
% tcp     computer time for the computation
%
% User m functions called:  none.
%----------------------------------------------

if (nargin > 9); force=1; else, force=0; end
if nargout ==3, tcp=clock; end
hbig=h*incout; 
t=(t0:hbig:tmax)'; n=length(t);
ns=(n-1)*incout; ts=t0+h*(0:ns)'; 
xnow=x0(:); vnow=v0(:); 
nvar=length(x0); 
jrow=1; jstep=0; h2=h/2;

% Form the inverse of the effective 
% stiffness matrix
mnv=h*inv(m+h2*(c+h2*k));

% Initialize the output matrix for x
x=zeros(n,nvar); x(1,:)=xnow'; 
zroforc=zeros(length(x0),1);

% Main integration loop
for j=1:ns
  tj=ts(j);tjh=tj+h2;
  if force 
    dv=feval(forc,tjh); 
  else 
    dv=zroforc; 
  end
  dv=mnv*(dv-c*vnow-k*(xnow+h2*vnow));
  vnext=vnow+dv;xnext=xnow+h2*(vnow+vnext);
  jstep=jstep+1;
  if jstep == incout
    jstep=0; jrow=jrow+1; x(jrow,:)=xnext';
  end
  xnow=xnext; vnow=vnext;
end
if nargout ==3
  tcp=etime(clock,tcp);
else
  tcp=[];
end

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

function [t,x,tcp] = ...
     mckde4i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
%
% [t,x,tcp]= ...
%    mckde4i(m,c,k,t0,x0,v0,tmax,h,incout,forc)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function uses a fourth order implicit 
% integrator with fixed stepsize to solve the 
% matrix differential equation
%           m x'' + c x' + k x = forc(t)
% where m,c, and k are constant matrices and 
% forc is an externally defined function.
%
% Input:
% ------
% m,c,k   mass, damping and stiffness matrices
% t0      starting time
% x0,v0   initial displacement and velocity
% tmax    maximum time for solution evaluation
% h       integration stepsize
% incout  number of integration steps between 
%         successive values of output
% forc    externally defined time dependent 
%         forcing function. This parameter 
%         should be omitted if no forcing 
%         function is used.      
%
% Output:
% -------
% t       time vector going from t0 to tmax 
%         in steps of h*incout
% x       matrix of solution values such 
%         that row j is the solution vector 
%         at time t(j)
% tcp     computer time for the computation
%
% User m functions called:  none.
%----------------------------------------------

if nargin > 9, force=1; else, force=0; end 
if nargout ==3, tcp=clock; end
hbig=h*incout; t=(t0:hbig:tmax)';
n=length(t); ns=(n-1)*incout; nvar=length(x0);
jrow=1; jstep=0; h2=h/2; h12=h*h/12;

% Form the inverse of the effective stiffness 
% matrix for later use.

m12=m-h12*k;
mnv=inv([[(-h2*m-h12*c),m12]; 
        [m12,(c+h2*k)]]);

% The forcing function is integrated using a 
% 2 point Gauss rule
r3=sqrt(3); b1=h*(3-r3)/6; b2=h*(3+r3)/6;

% Initialize output matrix for x and other 
% variables
xnow=x0(:); vnow=v0(:); 
tnow=t0; zroforc=zeros(length(x0),1);

if force 
  fnow=feval(forc,tnow); 
else 
  fnow=zroforc;
end
x=zeros(n,nvar); x(1,:)=xnow'; fnext=fnow;

% Main integration loop
for j=1:ns
  tnow=t0+(j-1)*h; tnext=tnow+h;
  if force
    fnext=feval(forc,tnext); 
    di1=h12*(fnow-fnext);
    di2=h2*(feval(forc,tnow+b1)+ ...
            feval(forc,tnow+b2));
    z=mnv*[(di1+m*(h*vnow)); (di2-k*(h*xnow))]; 
    fnow=fnext;
  else
    z=mnv*[m*(h*vnow); -k*(h*xnow)];
  end
  vnext=vnow + z(1:nvar); 
  xnext=xnow + z((nvar+1):2*nvar);
  jstep=jstep+1;

  % Save results every incout steps
  if jstep == incout
    jstep=0; jrow=jrow+1; x(jrow,:)=xnext';
  end

  % Update quantities for next step
  xnow=xnext; vnow=vnext; fnow=fnext;
end
if nargout==3
  tcp=etime(clock,tcp);
else 
  tcp=[]; 
end

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

function [m,k]=cablemk(masses,lngths,gravty)
%
% [m,k]=cablemk(masses,lngths,gravty)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Form the mass and stiffness matrices for 
% the cable.
%
% masses     - vector of masses
% lngths     - vector of link lengths
% gravty     - gravity constant
% m,k        - mass and stiffness matrices
%
% User m functions called:  none.
%----------------------------------------------

m=diag(masses);
b=flipud(cumsum(flipud(masses(:))))* ...
  gravty./lngths;
n=length(masses); k=zeros(n,n); k(n,n)=b(n);
for i=1:n-1
  k(i,i)=b(i)+b(i+1); k(i,i+1)=-b(i+1);
  k(i+1,i)=k(i,i+1);
end

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

function plterror(xmr,t2,h2,x2,T2,H2,X2,...
         t4,h4,x4,T4,H4,X4,tr2,Tr2,tr4,Tr4)
% plterror(xmr,t2,h2,x2,T2,H2,X2,...
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%          t4,h4,x4,T4,H4,X4,tr2,Tr2,tr4,Tr4)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Plots error measures showing how different
% integrators and time steps compare with 
% the exact solution using modal response.
%
% User m functions called:  none
%----------------------------------------------

% Compare the maximum error in any component
% at each time with the largest deflection
% occurring during the complete time history 
maxd=max(abs(xmr(:)));
er2=max(abs(x2-xmr)')/maxd;
Er2=max(abs(X2-xmr)')/maxd;
er4=max(abs(x4-xmr)')/maxd;
Er4=max(abs(X4-xmr)')/maxd;

plot(t2,er2,'-',T2,Er2,'--');
title(['Solution Error For Implicit ',...
       '2nd Order Integrator']);
xlabel('time');
ylabel('solution error measure');
lg1=['h= ', num2str(h2),  ...
     ', relative cputime= ', num2str(tr2)];
lg2=['h= ', num2str(H2),  ...
     ', relative cputime= ', num2str(Tr2)];
legend(lg1,lg2,2); figure(gcf); 
disp('Press [Enter] to continue'); pause
% print -deps deislne2

plot(t4,er4,'-',T4,Er4,'--');
title(['Solution Error For Implicit ',...
       '4th Order Integrator']);
xlabel('time');
ylabel('solution error measure');
lg1=['h= ', num2str(h4),  ...
     ', relative cputime= ', num2str(tr4)];
lg2=['h= ', num2str(H4),  ...
     ', relative cputime= ', num2str(Tr4)];
legend(lg1,lg2,2); figure(gcf); 
% print -deps deislne4 
disp(' '), disp('All Done')

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

% function [t,u,mdvc,natfrq]=...
%              udfrevib(m,k,u0,v0,tmin,tmax,nt)
% See Appendix B

Contact us