# Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd Edition

### Howard Wilson (view profile)

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

```