%
%   IVPNEX1.M    Simulation of a Simple Mechanical System
%       (Example for Numerical Integration of System of ODEs (IVPs))
%
%   This is a sample problem that considers a simple mechanical system 
%   consisting of a mass, spring, and dashpot (to represent viscous
%   damping).  The system also contains a specified forcing function, f(t).  
%   A force balance on the system gives the following second order ODE:
%       my''(t) + cy'(t) + ky(t) = f(t)
%
%   For specificity, let
%       m =1,   k = 2,   c = 3,   and   f(t) = 10exp(-2t)
%   and for initial conditions, let
%       y(0) = 0   and   y'(0) = -5
%   all with appropriate units, of course.
%
%   This can be put in matrix form by letting z1 = y  and  z2 = y', giving
%       z1'  =   z2
%       z2'  =  -(k/m)z1 - (c/m)z2 + (1/m)f(t)
%   or
%       z' = A*x + B*u   (where A & B are state space matrices and u is input))
%       
%   The goal of this example is to illustrate how to solve problems of this 
%   type within Matlab using the ODE23 function (uses the RK23 method).  The 
%   numerical solution is also compared to an analytical solution derived in the 
%   course Lecture Notes for 10/24.539.
%
%   This file also uses the following two files:
%       ivpnex1a.m - function file for ode23
%       ivpnex1b.m - function file for generating value of forcing function
%
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003).
%

%
%   getting started
      clear all, close all,  nfig = 0;
%
%   setup problem specific data
      m = 1;   k = 2;   c = 3;
      yo = 0;   ypo = -5;
%
%   setup constant matrices for state space formulation
      AA = [0 1; -k/m -c/m];   BB = [0 1/m]';
%
%   time domain of interest
      to = 0;   tf = 3;   te = linspace(to,tf,61);   
%
%   evaluate and plot forcing function over interval
      ff = ivpnex1b(te);
      nfig = nfig+1;   figure(nfig)
      plot(te,ff,'r-','LineWidth',2)
      title('IVPnex1:  Forcing Function Used in Simulation')
      xlabel('time'),ylabel('f(t)'),grid
%
%   solve 2x2 system of ODEs (Initial Value Prob -- IVP)
      zo = [yo ypo]';     
      tol = 0.001;    options = odeset('RelTol',tol);
      [t,z] = ode23('ivpnex1a',[to tf],zo,options,AA,BB);
%
%   exact solution (from lecture notes)
      te = te';    ze = zeros(length(te),2);
      ze(:,1) = -(10*te + 5).*exp(-2*te) + 5*exp(-te);
      ze(:,2) =  20*te.*exp(-2*te) - 5*exp(-te);
%
%   compare exact and numerical solution 
      nfig = nfig+1;   figure(nfig)
      plot(te,ze(:,1),'-g',te,ze(:,2),'-g',t,z(:,1),'or',t,z(:,2),'or','LineWidth',2)
      title('IVPnex1:  ODE23 (pts) vs Exact Solution (line) for Mechanical System')
      xlabel('time'),ylabel('position and velocity'),grid
      gtext('z1(t)'),gtext('z2(t)')
%
      range = axis;
      xt = range(1) + 0.52*(range(2)-range(1));
      yt1 = range(3) + 0.13*(range(4)-range(3)); 
      yt2 = range(3) + 0.08*(range(4)-range(3));
      text(xt,yt1,'z1(t) - position of mass versus time')
      text(xt,yt2,'z2(t) - velocity of mass versus time')
%
%   end of problem

