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)

...
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

Contact us