%
%   PENDULUM_2.M  Evaluate and plot the dynamics of a simple linear pendulum
%                  via the FD method (compared to the analytical solution)
%
%   This example evaluates and plots the angular position of a pendulum with a point 
%   mass at the end.  The focus is on the solution of IVPs using the finite difference
%   (FD) method.  For comparison, we also plot the analytical solution.  The equations
%   for the FD method are developed in the Lecture Notes following Lesson 4 (see 
%   fd_intro.pdf), with the basic mathematical model and analytical solution given in
%   the Lesson 1 Lecture Notes (also see pendulum_dynamics.pdf).  
%
%   A similar program, pendulum_1.m, focuses only on evaluating and plotting the 
%   analytical solution.
%
%   File prepared by J. R. White, UMass-Lowell (last update: Oct. 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   define problem parameters
      m = 1;                    % pendulum mass (kg)
      c = 2;                    % friction coeff (kg/s)
      L = 1;                    % length of pendulum arm (m)
      g = 10;                   % gravitational acceleration (m/s^2)
      poso = 30*pi/180;         % initial position (degrees converted to radians)
      velo = 0;                 % initial velocity (radians/s)
%
%       *************************
%       *  ANALYTICAL SOLUTION  *
%       *************************
%
%   evaluate roots of characteristic equation and compute equation coefficients
      a = -c/(2*m);   b = sqrt(g/L - a^2);
      c1 = poso;          c2 = (velo - a*c1)/b;
      d1 = b*c2 + a*c1;   d2 = a*c2 - b*c1;
%
%   define discrete time domain variable (with time in seconds)
      Nt = 101;   to = 0;   tf = 5;   t = linspace(to,tf,Nt);
%
%   evaluate analytical expressions for angular position 
      posa = exp(a*t).*(c1*cos(b*t) + c2*sin(b*t));
%
%       *****************
%       *  FD SOLUTION  *
%       *****************
%
%   evaluate equation coefficients and starting points
%   (can vary Nt to see the sensitivity of the numerical soln to value of dt)
      N = Nt-1;  dt = (tf-to)/N;  % dt is consistent with the time vector from above
      a = 1 + c*dt/(2*m);   b = 2 - g*dt*dt/L;   d = c*dt/(2*m) - 1;
      posn = zeros(1,Nt);   posn(1) = poso;    posn(2) = posn(1) + velo*dt;
%
%   now loop over all remaining discrete time points
      for i = 2:N 
        posn(i+1) = b*posn(i)/a + d*posn(i-1)/a;
      end
%
%       *************************
%       *  PLOT BOTH SOLUTIONS  *
%       *************************
%
      nfig = nfig+1;   figure(nfig)
      plot(t,posa*180/pi,'r-',t,posn*180/pi,'b:','LineWidth',2),grid    
      title('PENDULUM\_2: Pendulum Dynamics  --  \theta(t) vs t')
      xlabel('Time (sec)'),ylabel('Angular Position (degrees)')
      legend('Analytical Solution','FD Solution')
%
%   end of problem
