=====Math 445 lecture 21; Ordinary Differential Equations=====

Lecture involves 

  * physics => ordinary differential equation (ODE)
  * manipulating ODE until it fits form Matlab can integrate
  * Matlab's ode23, ode45 functions, plotting results
  * perhaps some idea of how Matlab integrates ODE 

=====Physics=====

Pendulum physics, F = ma with linear air resistance

<latex>
 m l \theta'' = -mg sin \theta - b \theta'     

 \theta'' + c \theta' + g/l \sin \theta = 0    
</latex>

This is a 2nd-order nonlinear differential equation in one variable, theta.
2nd order b/c $\theta''$ term. Nonlinear b/c $\sin \theta$. Linearization for small 
theta would approximate sin theta = theta for small $\theta$. No big thing

=====Get ODE into Matlab form=====

Matlab can integrate 1st-order nonlinear ODEs in n variables, of form

  dx/dt = f(t,x)

for vector x. Fortunately we can always convert an nth-order ODE in 1 variable
to a 1st order ODE in n variables, as follows. Let

<latex>
  x1 = theta
  x2 = theta'
</latex>

Then

  dx1/dt = d theta/dt = x2

  dx2/dt = d theta'/dt = theta'' 
         = - c theta' - g/l sin theta
         = -c x2 - g/l sin x1

So for vector x = (x1, x2),

  dx/dt = f(t,x) = (x2,  -c x2 - g/l sin x1)

=====Integrating ODE in Matlab =====

Let's code this as an anonymous function in Matlab, setting constants prior

  c = 0;    % no air resistance
  g = 9.8;  % gravity in mks (m/s^2)x
  l = 1.0;  % one meter pendulum 

  f = @(t,x) [x(2); -c*x(2) - g/l * sin(x(1))];

  % integrate system of ODEs from t=0 to t=100, from initial condition
  % [x1; x2] = [theta; theta'] = [pi/10; 0], using Matlab's ode45

  [t, x] = ode45(f, [0 20], [pi/10; 0])

  % outputs are an 
  % N-vector t of timesteps
  % N x 2 matrix x whose cols are x1, x2 = theta, theta' at the N timesteps
  % e.g. x(:,1) is N-vector of theta at times corresponding to N-vector t. 

  % can graph theta versus t as follows

  plot(t, x(:,1), 'b')