Introduction to PDEs
====================

We start with a model for air quality.

Air Quality Modeling
--------------------

Consider the following context.

* An industrial plant emits noxious fumes.

* Suppose once an hour, the plant emits fumes for a few minutes.

* One mile away is our house.

Problem: *If the wind is blowing from the plant to the house,
what will the maximum concentration of the noxious fumes be
when they reach your house?* This problems is taken
from chapter 2 of *Industrial Mathematics:
A Course in Solving Real-World Problems*
by Avner Friedman and Walter Littman, 
published by SIAM in 1994.

We assume bell shaped profiles of the fumes,
shown in :numref:`figBellShapedFumes`.

.. index:: air quality

.. _figBellShapedFumes:

.. figure:: ./figBellShapedFumes.png
   :align: center

   At the left is :math:`c_0(x)` is the profile of the fumes 
   produced at time zero, :math:`c(x, t)` is the profile 
   at place :math:`x` at time :math:`t`, at the right.

In the proposal of the solution, we make the following considerations.
Assume the wind blows at constant speed :math:`v` towards our house.
The speed multiplied by time traveled equals distance traveled.
So :math:`v~\!t` is the distance traveled 
at speed :math:`v` at time :math:`t`.
The proposed solution is then

.. math::

   c(x,t) = c_0(x - v~\!t).

The original profile of the fumes shifts in the *opposite*
direction of the wind, whence the minus sign in :math:`x - v~\!t`.

A shifted fume profile is shown in :numref:`figShiftedFumeProfile`.

.. _figShiftedFumeProfile:

.. figure:: ./figShiftedFumeProfile.png
   :align: center

   At the left is the plot of :math:`\exp(-100 x^2)`,
   and the shifted profile is :math:`\exp(-100 (x-5)^2)`,
   shown at the right.

The derivatives of the proposed solution

.. math::

   c(x,t) = c_0(x - v~\!t)

with respect to space and time are

.. math::

   \begin{array}{rcl}
   {\displaystyle \frac{\partial}{\partial x} c(x,t)}
   & = & c'_0(x - v~\!t) \cdot 1,  \\
   {\displaystyle \frac{\partial}{\partial t} c(x,t)}
   & = & c'_0(x - v~\!t) \cdot (-v).
   \end{array}

What follows behind the :math:`\cdot$ after $c'_0`
comes from the chain rule.

Eliminating :math:`c'_0(x - v~\!t)` 
yields *the advection equation*, which is:

.. index:: advection equation, transport equation

.. math::

   \frac{\partial c}{\partial t}(x,t)
   + v \frac{\partial c}{\partial x}(x,t) = 0, \quad c(x, 0) = c_0(x).

The Partial Differential Equation (PDE) we just derived
is also known as *the transport equation*.

The solution :math:`c(x,t)` is the concentration 
at location :math:`x` at time :math:`t`,
traveling at constant speed :math:`v`.
As an application, consider traffic where 
:math:`\displaystyle c(x,t)
= \frac{\mbox{number of vehicles}}{\mbox{unit mile}}`.
In :numref:`figSpeedConcentration` we see speed as a function
of the concentration.

.. _figSpeedConcentration:

.. figure:: ./figSpeedConcentration.png
   :align: center

   Speed :math:`v` in function of concentration :math:`c`,
   where :math:`v_f` is the free traffic speed, and
   :math:`c_j` is the jam density.

Traffic slows down when the number of vehicles increases.
Thus, the speed :math:`v` depends on the concentration :math:`c`:

* :math:`v_f` is the *free traffic speed*, when the road is free.

* :math:`c_j` is the *jam density*, :math:`v=0` at a traffic jam.

The application of traffic motivates the development
of the transport equation with variable speed.
The variable in the transport equation is the concentration :math:`c(x,t)`
and the partial differential equation has two terms:

1. the change in concentration :math:`c` in time :math:`t`,

2. the change in concentration :math:`c` in space :math:`x`.

With variable speed :math:`v`, the equation is then

.. math::

   \frac{\partial c}{\partial t}(x,t)
   + \frac{\partial (v~\! c)}{\partial x}(x,t) = 0, \quad c(x, 0) = c_0(x).

Pulling :math:`v` into the 
:math:`\displaystyle \frac{\partial}{\partial x}` term
agrees with :math:`c(x,t) = c_0(x - v~\!t)`, 
proposed in the traveling profile of the fumes in the air quality model.

We study traffic flow in greater detail in a later lecture.

In moving to three dimensions, operators notations are useful.
:math:`c(x_1, x_2, x_3, t)` is the concentration at the point 
:math:`(x_1, x_2, x_3)`, at time :math:`t`.
The velocity vector of the transport is

.. math::

   \vec{v}(x_1, x_2, x_3, t) = 
   v_1 {\bf i} + v_2 {\bf j} + v_3 {\bf k},

where :math:`v_1`, :math:`v_2`, :math:`v_3`
are functions of :math:`(x_1, x_2, x_3, t)`.

.. topic:: Definition of the del operator

   The *del operator* is
   :math:`\displaystyle \nabla = 
   {\bf i} \frac{\partial}{\partial x}
   + {\bf j} \frac{\partial}{\partial y}
   + {\bf k} \frac{\partial}{\partial z}`.

The three dimensional advection equation is then

.. math::

   \frac{\partial c}{\partial t} + \nabla \cdot ( \vec{v}~\!c ) = 0,

formulated with the :index:`del operator` and the :index:`dot product`

.. math::

   \nabla \cdot \vec{v} = \frac{\partial v_1}{\partial x}
                        + \frac{\partial v_2}{\partial y}

Let :math:`(x_0, y_0)` be the origin of the fumes.
Consider the radius :math:`R`, with

.. math::

   R^2 = (x - x_0)^2 + (y - y_0)^2.

Then the initial concentration

.. math::

   c_0(x,y) = 
   \left\{
   \begin{array}{cc}
   {\displaystyle 50 
   \left( 1 + \cos \left( \frac{\pi R}{4} \right) \right)}
   & \mbox{if } R < 4, \\
   0 & {\mbox{otherwise,}}
   \end{array}
   \right.

is referred to as a :index:`cosine hill`.

Diffusion
---------

The model in :numref:`figShiftedFumeProfile` is not realistic,
because it does not incorporate the change in concentration.
Consider the evolution of the profile of fumes
in :numref:`figFumesDiffusion`.

.. _figFumesDiffusion:

.. figure:: ./figFumesDiffusion.png
   :align: center

   The diffusion of the profile of fumes.

The transition from high to low concentration is *diffusion*.
To introduce a plausible mathematical model for :index:`diffusion`,
consider the following symbolic calculation:

::

   >>> from sympy import var, diff, exp
   >>> x, t = var("x, t")
   >>> c = exp(-t)*exp(I*x)
   >>> c
   exp(-t)*exp(I*x)
   >>> d1t = diff(c,t)
   >>> d1t
   -exp(-t)*exp(I*x)
   >>> d1x = diff(c,x)
   >>> d1x
   I*exp(-t)*exp(I*x)
   >>> d2x = diff(d1x, x)
   >>> d2x
   -exp(-t)*exp(I*x)
   >>> d1t - d2x
   0

A plausible model works with a decaying exponential in time:

.. math::

   c(x,t) = e^{-t} e^{I~\!x}, \quad I = \sqrt{-1}.

Then we have

.. math::

   \frac{\partial c}{\partial t} = - e^{-t} e^{I~\!x}

and

.. math::

   \frac{\partial^2 c}{\partial x^2} = - e^{-t} e^{I~\!x},
   \quad \mbox{as } I^2 = -1.

Thus:

.. index:: heat equation

.. math::

   \frac{\partial c}{\partial t} = \frac{\partial^2 c}{\partial x^2}

is an example of the *heat equation*.

.. topic:: Definition of the Laplacian operator

   The *Laplacian operator* :math:`\nabla^2` 
   applied to :math:`f = f(x,y,z)` is :math:`\nabla \cdot \nabla f`.

.. math::

   \nabla f = \frac{\partial f}{\partial x} {\bf i}
            + \frac{\partial f}{\partial y} {\bf j}
            + \frac{\partial f}{\partial z} {\bf k}
   \mbox{ is the gradient vector of } f.

.. math::

   \begin{array}{rcl}
   \nabla \cdot \nabla f
   & = & {\displaystyle
   \left( {\bf i} \frac{\partial}{\partial x}
        + {\bf j} \frac{\partial}{\partial y}
        + {\bf k} \frac{\partial}{\partial z} \right) \cdot \nabla f} \\
   & = & {\displaystyle
    \frac{\partial^2 f}{\partial x^2}
    + \frac{\partial^2 f}{\partial y^2}
    + \frac{\partial^2 f}{\partial z^2}}
   \end{array}

The :index:`heat equation`, or the equation of :index:`diffusion` is then

.. math::

   \frac{\partial u}{\partial t} = \nabla^2 u, \quad u = u(x,y,z,t).

Consider the evolution of the profile of fumes
in :numref:`figFumesEvolution`.

.. _figFumesEvolution:

.. figure:: ./figFumesEvolution.png
   :align: center

   The evolution of the profile of fumes with diffusion.

The transition from high to low concentration is modeled by
the *advection diffusion equation*:

.. math::

   \frac{\partial c}{\partial t}(x,t)
   + \frac{\partial (v~\! c)}{\partial x}(x,t)
   = k \frac{\partial^2 c}{\partial x^2}, \quad k > 0,
   \quad c(x, 0) = c_0(x).

For :math:`c = c(x_1, x_2, x_3, t)`, 
the :index:`advection diffusion equation` is

.. math::

   \frac{\partial c}{\partial t} + \nabla \cdot ( \vec{v}~\!c )
   = \sum_{i=1}^3 \frac{\partial}{\partial x_i}
   \sum_{j=1}^3 \left( k_{i,j} \frac{\partial c}{\partial x_j} \right),

where 

* :math:`k = (k_{i,j})` is a positive definite matrix, and

* :math:`c(x_1, x_2, x_3, 0) = c_0(x_1, x_2, x_3)`
  is the initial concentration.

Proposals for Topics of a Project
---------------------------------

1. numerically solve the advection diffusion equation

   The advection diffusion equation is solved in
   chapter 2 of *Industrial Mathematics:
   A Course in Solving Real-World Problems*
   by Avner Friedman and Walter Littman, 
   published by SIAM in 1994.

   Study chapter 2 of the Friedman-Littman book.

   Review section 8.1
   of *Numerical Analysis* of T. Sauer (MCS 471).

   Do the following:

   * Describe the solution using your own words.

   * Do the computational experiments to examine the stability.

   * Write a report on your computational experiments.

2. Fick's laws of diffusion

   The heat equation was derived by Adolf Fick in 1855.

   Study the wikipedia page on the derivation of Fick's second law
   and consult the references.

   Do the following:

   * Describe the derivation of the heat equation.

   * Reproduce the illustration on the wikipedia page
     on Fick's laws of diffusion.

   * Write a report on your computational experiments.

3. population density of age segment

   We are interested to model the population of an age segment.

   Let :math:`u = u(x,t)` be the population density defined as

   .. math::

      P(b, t) - P(a, t) = \int_a^b u(x, t) dx

   the population at time :math:`t` between 
   the ages :math:`a` and :math:`b`.
   Let :math:`q(x)` be the probability of death less net immigration
   at age :math:`x` per year.

   Do the following:

   * Show that :math:`u_t + u_x + q u = 0`.

   * Solve analytically (see page 227 of the text book).

   * Simulate the model numerically and interpret the results.

Exercises
---------

1. Let :math:`f(x)` be continuously differentiable and consider

   .. math::

      \frac{d~\!x}{d~\!t} = f(x), \quad t > 0, \quad x(0) = x_0.

   Its solution :math:`x(t)` depends on :math:`x_0`,
   so we write :math:`x(t, x_0)`
   as the unique curve that passes through the point :math:`(x_0, 0)`.

   Show that :math:`x(t, x_0)` is differentiable and that

   .. math::

      z(t) = \frac{\partial x(t, x_0)}{\partial x_0}

   satisfies

   .. math::

      \frac{d~\!z}{d~\!t} = f_x(x(t, x_0)) z, 
      \quad f_x = \frac{d~\!f}{d~\!x},
      \quad z(0) = 1.

2. Consider :math:`f(x,y)` as the value of the gray scale in a picture.
   The problem is to detect the edges.

   Use divided differences for the second derivatives in

   .. math::

      \nabla^2 f =
      \frac{\partial^2 f}{\partial x^2}
      + \frac{\partial^2 f}{\partial y^2}.

   Why is the application of this operator useful to detect edges?

3. What is the Laplacian operator in polar coordinates?

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

1. Avner Friedman and Walter Littman, chapter 2 of
   *Industrial Mathematics: 
   A Course in Solving Real-World Problems*, SIAM 1994.

2. Charles R. MacCluer, start of Chapter 11 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.