Dynamics in Mechanics and Biology
=================================

We introduce models in mechanics and biology,
using numerical methods to solve the differential equations.

The Planar Pendulum
-------------------

We introduce numerically solving ordinary differential equations
with a system of first order equations, modeling a pendulum.
Introducing an auxiliary variable :math:`v(t) = \theta'(t)`,
we transform 

.. math::

   \theta''(t) + \alpha \sin(\theta(t)) = 0, \quad \alpha = g/L,

into a system of first order differential equations,
shown in :numref:`figPendulumSystem`.
The initial conditions are :math:`\theta(0) = \pi/6`
and :math:`v(0) = 0`.

.. _figPendulumSystem:

.. figure:: ./figPendulumSystem.png
   :align: center

   A system of first order equations to model a pendulum.

To solve this problem numerically in Julia,
we use the package ``DifferentialEquations.jl``:

::

   julia> using DifferentialEquations

``DifferentialEquations.jl`` is a Julia package which exports
a numerical solver for ordinary differential equations,
described in the paper by
Christopher Rackauckas and Qing Nie:
**Differentialequations.jl -- a performant and feature-rich 
ecosystem for solving differential equations in julia**.
*Journal of Open Research Software* 5(1), 2017.

The definition of the right hand side function is done in
the function below:

::

   """
       function pendrhs!(du, u, p, t)
  
   defines the right hand side of the system of first order
   order differential equations of a model of a planar pendulum.

   The function updates the two components of du,
   which respectively define the velocity and the acceleration,
   while u[1] contains the value for the angle
   and u[2] the value of the velocity.
   The p is the array with the parameters
   and t is the independent time variable.
   """
   function pendrhs!(du, u, p, t)
       angle, velocity = u
       alpha = p
       du[1] = velocity          # instead of u[2]
       du[2] = -alpha*sin(angle) # angle is u[1]
   end

The Initial conditions, time span, and parameter are
set in the steps below:

1. The initial conditions are :math:`\theta_0 = \pi/6`, 
   :math:`v_0 = 0`.

   ::

       theta0 = pi/6  # initial angle
       v0 = 0 # initial velocity
       initconds = [theta0, v0]

2. The time span is :math:`[0, 10]`.

   ::

       timespan = (0.0, 10.0)

3. The value for the parameter :math:`\alpha = 1`.

   ::

       alpha = 1.0

Then we can define the problem:

::

    ode = ODEProblem(pendrhs!, initconds, timespan, alpha)

The solving and plotting of the solution is done by
the instructions below.

::

    sol = solve(ode)

    plot(sol, xaxis="t", label=["theta" "velocity"], layout=(2,1))


The plotted solution is shown in :numref:`figODEpendulum`.

.. _figODEpendulum:

.. figure:: ./figODEpendulum.png
   :align: center

   Displacement and velocity of a pendulum trajectory.

Damped Harmonic Motion
----------------------

A model of a suspension mechanism, using the equation for
:index:`damped harmonic motion`, written in engineering notation, is:

.. math::

   \frac{d^2}{d~\!t^2} y(t) + 2 \zeta \omega 
   \left( \frac{d}{d~\!t} y(t) \right) + \omega^2 y(t) = 0.

The symbolic solution is

.. math::

   y(t) = C_1 e^{\omega t \left( -\zeta - \sqrt{\zeta^2 - 1} \right)}
        + C_2 e^{\omega t \left( -\zeta + \sqrt{\zeta^2 - 1} \right)}.

Observe the discriminant :math:`\zeta^2 - 1 = 0`.

To determine :math:`C_1` and :math:`C_2`,
we impose the initial conditions:

.. math::

   y(0) = 1 \quad \mbox{and} \quad y'(0) = 0.

A case analysis is shown in :numref:`figUnderDamped`,
:numref:`figCriticallyDamped`, :numref:`figOverDamped`, and
:numref:`figFasterDamped`, showing respectively underdamped, critically damped,
overdamped, and faster damped motion.

.. _figUnderDamped:

.. figure:: ./figUnderDamped.png
   :align: center

   Underdamped motion, :math:`\zeta = 1/2`, :math:`\omega = 1`

.. _figCriticallyDamped:

.. figure:: ./figCriticallyDamped.png
   :align: center

   Critically damped motion, :math:`\zeta = 1`, :math:`\omega = 1`

.. _figOverDamped:

.. figure:: ./figOverDamped.png
   :align: center

   Overdamped motion, :math:`\zeta = 3/2`, :math:`\omega = 1`

.. _figFasterDamped:

.. figure:: ./figFasterDamped.png
   :align: center

   Faster damped motion, :math:`\zeta = 1`, :math:`\omega = 2`

Population Dynamics
-------------------

A differential equation to model the growth of a population 
was proposed by Verhulst in the 19th century.  
We use the following notation:

* Let :math:`N(t)` be the number of the population
  in function of time :math:`t`.

* The carrying capacity is :math:`K`.

* The birth rate is :math:`r`.

.. math::

   \frac{d}{d~\!t} N(t) = r \left( 1 - \frac{N(t)}{K} \right) N(t).

The rate of increase in population

* is proportional to the birth rate times the population number,

* multiplied with a term to take the carrying capacity into account.

The solution to this :math:`logistic growth` model
for the parameters :math:`K=50`, :math:`N_0 = 8`, :math:`r = 0.7`
is shown in :numref:`figVerhulst`.

.. _figVerhulst:

.. figure:: ./figVerhulst.png
   :align: center

   A model of logistic growth.

We consider next a model of growth with a :index:`predator`,
the :index:`population dynamics` of the spruce budworm.
Birds eat worms.
From the logistic growth model we subtract a predator term.
The predator term is a function of the population number: :math:`P(N(t))`.

.. math::

   \frac{d}{d~\!t} N(t) = r \left( 1 - \frac{N(t)}{K} \right) N(t)
                        - P(N(t)).

The population of the worms should reach a critical value
before the birds start noticing the worms.
Once this happens, the birds eat worms in greater numbers.
Once the population falls below another critical value,
the predation disappears.

The predator term :math:`B n^2/(A^2 + n^2)`,
for :math:`A = 30` and :math:`B = 4` is
shown in :numref:`figWormPredatorShape`.

.. _figWormPredatorShape:

.. figure:: ./figWormPredatorShape.png
   :align: center

   The shape of the predator term.

Our model

.. math::

   \frac{d}{d~\!t} N(t) = r \left( 1 - \frac{N(t)}{K} \right) N(t)
       - \frac{B (N(t))^2}{A^2 + (N(t))^2}

has too many parameters.  Applying the substitutions

.. math::

   N(t) = A u(t), \quad r = \frac{B r}{A}, \quad K = q A

leads to 

.. math::

   \frac{d}{d~\!t} A u(t) = B r \left( 1 - \frac{u(t)}{q} \right) u(t)
       - \frac{A^2 B (u(t))^2}{A^2 + A^2 (u(t))^2}.

Observe that the :math:`A^2` cancels in the predator term.

With :math:`u(t) = u(t B/A)` we can remove :math:`B`

.. math::

   \frac{d}{d~\!t} A u(t B/A)
    = B r \left( 1 - \frac{u(t B/A)}{q} \right) u(t)
    - \frac{A^2 B (u(t B/A))^2}{A^2 + A^2 (u(t B/A))^2}

because :math:`B` appears at the left and the right of the equation.
Then we replace :math:`B t/A` by :math:`t`:

.. math::

   \frac{d}{d~\!t} u(t)
    = r \left( 1 - \frac{u(t)}{q} \right) u(t)
    - \frac{(u(t))^2}{1 + (u(t))^2},

with

.. math::

   r = \frac{B}{A} r \quad \mbox{and} \quad q = \frac{K}{A}.

At the :index:`steady states` of

.. math::

   \frac{d}{d~\!t} u(t)
    = r \left( 1 - \frac{u(t)}{q} \right) u(t)
    - \frac{(u(t))^2}{1 + (u(t))^2}

the right hand side equals zero.
Replace :math:`u(t)` by :math:`x` and we find

.. math::

   r \left(1 - \frac{x}{q} \right) x - \frac{x^2}{1 + x^2} = 0.

Observe that :math:`x = 0` is a trivial factor.

The steady states for the parameters 
:math:`r = 0.4` and :math:`q = 20` are
shown in :numref:`figWormSteadyStates`.

.. _figWormSteadyStates:

.. figure:: ./figWormSteadyStates.png
   :align: center

   Steady states of a population model with a predator term.

To find the numerical roots, 
applying ``solve()``, for :math:`r = 0.4` and :math:`q = 20` gives

::

    3-element Vector{Sym}:
    0.480535939911864 - 0.e-22*I,
    2.43633300316055 + 0.e-19*I,
    17.0831310569276 - 0.e-20*I

Let us now compute a special model.
Symbolic solving of

.. math::

   \frac{d}{d~\!t} u(t)
    = 0.4 \left( 1 - \frac{u(t)}{20} \right) u(t)
    - \frac{(u(t))^2}{1 + (u(t))^2}

is not possible.  We use ``DifferentialEquations.jl``
and define the right hand size of the system of differential equations
in the function below.

::

   """
       function wormrhs!(du, u, p, t)

   defines the right hand side for a special worm model.
   """
   function wormrhs!(du, u, p, t)
       growth = 0.4*(1 - u[1]/20)*u[1]
       predator = - (u[1]^2)/(u[1]^2 + 1)
       du[1] = growth + predator
   end

The solution of this special model is shown
in :numref:`figWormSpecial`.

.. _figWormSpecial:

.. figure:: ./figWormSpecial.png
   :align: center

   The solution of a special population model with predator term.

Observe that we started at 2.4,
very close to a steady state by ended up at 0.4.
Consider again :numref:`figWormSteadyStates`.
As :math:`r` increases from 0.4 to 0.5 while :math:`q=20` remains fixed,
the two steady states become closer to each other and then disappear.
As :math:`r` increases, we get a catastrophic outbreak of the worm population.

.. index:: catastrophic jump, hysteresis

Image now that we would start at :math:`r=0.6`, at a large steady state.
The return to a smaller value of the population will only be possible
when the large steady state, jointly with the unstable middle steady
state disappears as the predator line sinks below the growth line.
The catastrophic jumps happen at different values of the parameter 
:math:`r` and are depending on the direction of the parameter change.
This is called *hysteresis*.

A catastrophe in visual perception is shown
in :numref:`figCatastrophePerception`.

.. _figCatastrophePerception:

.. figure:: /figCatastrophePerception.png
   :align: center

   Look at the pictures below from left to right
   and then from right to left.
   Where the sudden jump in perception happens depends on the order,
   whether we look left-to-right or right-to-left.

Proposal for a Topic of a Project
---------------------------------

One deficiency of population models is that they do not take
the age structure into account.

Read section 1.7 of 
*Mathematical Biology I.  An Introduction*
by James D. Murray, Springer-Verlag 2002.

Consider the following.

* Work through section 1.7 computationally and
  build a numerical model to reproduce the plots in that section.

* What is the effect of an aging population on its growth?

Exercises
---------

1. An extended model of the pendulum adds damping:

   .. math::

      \frac{d^2 \theta}{d t^2} + A \frac{d \theta}{d t} 
      + \alpha \sin(\theta) = 0, \quad A = 1/10, 

   with initial values :math:`\theta(0) = \pi/6`,
   :math:`\theta'(0) = \pi/6`.

   * Define the ODE to solve with ``DifferentialEquations.jl``.

   * Solve the ODE and plot the trajectories of angles and velocities.

2. Another extended model of the pendulum adds damping and forcing:

   .. math::

      \frac{d^2 \theta}{d t^2} + A \frac{d \theta}{d t} 
      + \alpha \sin(\theta) = B \sin(t), \quad A = 1/10, \quad B = 10,

   with initial values :math:`\theta(0) = \pi/6`,
   :math:`\theta'(0) = \pi/6`.

   * Define the ODE to solve with ``DifferentialEquations.jl``.

   * Solve the ODE and plot the trajectories of angles and velocities.

3. Consider

   .. math::

      \frac{d^2 y}{d~\!t ^2} + \frac{1}{10} \frac{d~\!y}{d~\!t} + 4 y = 0,
      \quad y(0) = 1, \quad y'(0) = 0.

   * Compute the solution trajectory for this problem.

   * How long before :math:`y(t) \leq 0.01`?

4. Consider :math:`u(t)` as the solution of

   .. math::

      \frac{d~\!u}{d~\!t} = r~\!u \left( 1 - \frac{u}{q} \right)
      - \frac{u^2}{1 + u^2},

   where :math:`r` and :math:`q` are parameters with positive values.

   1. Show that the regions for which there is 1 positive steady state,
      or 3 positive steady states, are described by

      .. math::

         r = \frac{2 a^3}{(1 + a^2)^2} \quad \mbox{and} \quad
         q = \frac{2 a^3}{a^2 - 1},

      for some parameter :math:`a`.

   2. Solve the problem for :math:`a = \sqrt{3}` and make a plot.

   3. Interpret the solution, for when :math:`a \rightarrow \infty`?

   4. Interpret the solution, for when :math:`a \rightarrow 1`?

Bibliography
------------

1. Charles R. MacCluer, middle of Chapter 9 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.

2. James D. Murray:
   *Mathematical Biology I.  An Introduction*.
   Third Edition, Springer-Verlag 2002.