%
%   HW14_DEMO1.M  Compute numerical + symbolic solution to HW14_DEMO1
%
%   This problem was derived from Prob. 12 pg. 196 of Edwards and Penny (Prentice 
%   Hall, 1996).  The first part addresses the NUMERICAL solution of a particular 
%   mechanical oscillation problem and the solution is compared to an analytical  
%   solution developed in a previous class demo.  ODE23 is used as the ODE solver.
%
%   The force balance equation for this problem is
%       x'' + 6x' + 13x = 10*sin(5t)
%   with   x(0) = 0   and x'(0) = v(0) = 0
%   This can be written as two first order equations and evaluated as an IVP with 
%   a standard ODE solver.
%
%   Calls HW14_DEMO1_EQN.M to evaluate desired vector function, f(z,t).  Here z is a
%   2-element vector with z(1) = x and z(2) = x' = v.  With this definition, the 
%   above equation can be written as
%       dz(1)/dt = z(2)
%       dz(2)/dt = -13*z(1)-6*z(2)+10*sin(5*t)
%
%   The second part of the problem gives the SYMBOLIC solution using Matlab's 
%   dsolve command.  This symbolic solution is simply printed to the screen.
%
%   File prepared by J. R. White, UMass-Lowell (April 2005)
%
      
%
%   getting started
      clear all;   close all;   nfig = 0;
%
%   NUMERICAL Solution
%
%   define integration variables and initial conditions
      to = 0;   tf = 5;     % time in seconds
      xo = 0;   vo = 0;     % initial position (m) and speed (m/s)
      zo = [xo vo]';        % initial conditions in vector form
%
%   Evaluate x(t) and v(t) (ie. components of vector z) numerically
      [t,z] = ode23('hw14_demo1_eqn',[to tf],zo);
%
%   Evaluate x(t) analytically (amplitude and phase angle form)
      te = linspace(to,tf,200);
      xe = 0.7737*exp(-3*te).*cos(2*te-1.1903) - 0.3095*cos(5*te-0.3805);
%
%   Plot x(t) profiles
      nfig = nfig+1;   figure(nfig)
      plot(t,z(:,1),'ro',te,xe,'g-','LineWidth',2),grid    
      ylabel('Position (meters)'),xlabel('Time (seconds)')
      title('HW14\_Demo1:  Numerical and Analytical Solution Profiles')
      legend('ODE23 Solution','Analytical Solution')
%
%   SYMBOLIC Solution
%
      S = dsolve('D2x + 6*Dx + 13*x = 10*sin(5*t)', 'x(0) = 0','Dx(0) = 0','t');
      disp('Symbolic Solution to HW14_Demo1: '), pretty(simple(S)), disp(' ')
%
%   end of problem
