Separation of Variables
=======================

A partial differential equation can be reduced
to ordinary differential equations via the method
of the separation of variables.

Heat Diffusion
--------------

Consider the cooling an insulated tube at the ends.
Imagine a tube filled with hot water.
The tube is insulated everywhere, except at the ends.
*At the uninsulated ends of the tube cooling happens,
what is the evolution of the temperature distribution?*

To solve this problem, we develop
a mathematical model for :index:`heat diffusion`.

Let :math:`u(x, t)` be the temperature 
at position :math:`x \in [0,1]` at time :math:`t \geq 0`.

The heat diffusion is modeled by a partial differential equation

.. math::

   \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2},

where

* :math:`u(x, 0) = 1` is the initial condition, temperature is everywhere 1,

* :math:`u(0, t) = 0` is the left boundary condition, temperature at zero,

* :math:`u(1, t) = 0` is the right boundary condition, temperature at zero.

The geometric interpretation is as follows:
*The rate at which the temperature decreases is proportional
to the concavity of the temperature distribution.*

A plausible temperature distribution is shown
in :numref:`figSineProfile`.

.. _figSineProfile:

.. figure:: ./figSineProfile.png
   :align: center

   A sine profile of the temperature distribution.

Decoupling Space from Time
--------------------------

Given is a partial differential equation in :math:`u = u(x,t)`,
in two variables:

* :math:`x` specifies the position, *where* :math:`u` is,

* :math:`t` is a moment in time, *when* :math:`u` exists.

In an attempt to find a solution, we simplify the problem,
assuming position and time to be independent from each other.
This assumption implies that we can write :math:`u` as a product of

1. a function in :math:`x`, and

2. a function in :math:`t`.

This method reduces the given partial differential equation in two
decoupled ordinary differential equations.
This method is called *separation of variables*.

Assuming that :math:`x` and :math:`t` are independent from each other:

.. math::

   U(x, t) = T(t) X(x),

implies that we propose :math:`U`,
a product of two functions :math:`T` and :math:`X`,
as the solution :math:`u` to the partial differential equation
:math:`\displaystyle \frac{\partial u}{\partial t}
= \frac{\partial^2 u}{\partial x^2}`.

Now we impose the partial differential equation on :math:`U(x, t)`:

.. math::

   \frac{\partial U}{\partial t} = \frac{\partial T}{\partial t} X(x)
   \quad \mbox{and} \quad
   \frac{\partial^2 U}{\partial x^2} = T(t) \frac{\partial^2 X}{\partial x^2}.

Then the original partial differential equation becomes

.. math::

   \frac{\partial T}{\partial t} X(x)
   = T(t) \frac{\partial^2 X}{\partial x^2}.

Now comes the decoupling into ordinary differential equations.
The above equation is rewritten as

.. math::

   \frac{\partial T}{\partial t} \frac{1}{T(t)}
   = \frac{1}{X(x)} \frac{\partial^2 X}{\partial x^2}.

The assumption that space and time are independent from each other
implies the decoupling of the equation,
wedging in a constant :math:`\lambda`:

.. math::

   \frac{\partial T}{\partial t} \frac{1}{T(t)}
   = \lambda =
   \frac{1}{X(x)} \frac{\partial^2 X}{\partial x^2}.

The constant :math:`\lambda` is independent
of :math:`x` and :math:`t`.
This then leads to two ordinary differential equations:

.. math::

   \frac{\partial T}{\partial t}
   = \lambda T(t) \quad \Rightarrow \quad T(t) = T(0) e^{\lambda t}

and

.. math::

   \frac{\partial^2 X}{\partial x^2} = \lambda X(x).

Because we model a :index:`cooling off` process, we want decay:

.. math::

   \lambda = - \omega^2 < 0.

The choice of symbol for the separation constant :math:`\lambda`
suggest a connection with the :index:`eigenvalue problem`.
Using the :index:`Laplacian` :math:`\nabla^2`, we abbreviate
:math:`\displaystyle \frac{\partial^2 X}{\partial x^2} =  \lambda X(x)`
as :math:`\nabla^2 X = \lambda X`.

The function :math:`X(x)` is an *eigenfunction* of :math:`\nabla^2`
corresponding to :math:`\lambda`.

As a solution to

.. math::

   \frac{\partial^2 X}{\partial x^2} = - \omega^2 X(x),

consider a linear combination of a cosine and a sine function:

.. math::

   \begin{array}{rcl}
   X(x) & = & a \cos(\omega x) + b \sin(\omega x) \\
   X'(x) & = & -a~\!\omega \sin(\omega x) + b~\!\omega \cos(\omega x) \\
   X''(x) & = & -a~\!\omega^2 \cos(\omega x) - b~\!\omega^2 \sin(\omega x) \\
          & = & -\omega^2 X(x).
   \end{array}

Therefore, for some constants :math:`a` and :math:`b`,

.. math::

   X(x) = a \cos(\omega x) + b \sin(\omega x)

is a solution to the second order differential equation.

To arrive at the general form of a solution, combining

.. math::

   T(t) = T(0) e^{- \omega^2 t}

and

.. math::

   X(x) = a \cos(\omega x) + b \sin(\omega x)

leads to

.. math::

   U(x, t) = T(t) X(t)
           = T(0) e^{-\omega^2 t}
   \left( \vphantom{e^{-\omega^2 t}} a 
          \cos(\omega x) + b \sin(\omega x) \right).

This solution depends on :math:`T(0)`, :math:`\omega`,
:math:`a`, and :math:`b`.

Let :math:`T(0)a = A` and :math:`T(0)b = B`, 
then the general form becomes

.. math::

   U(x, t) = e^{-\omega^2 t}
   \left( \vphantom{e^{-\omega^2 t}} A 
          \cos(\omega x) + B \sin(\omega x) \right).

On the general form of a solution
we impose the boundary conditions

.. math::

   u(0, t) = 0 \quad \mbox{and} \quad u(1, t) = 0.

The condition :math:`u(0, t) = 0` leads to

.. math::

   U(0, t) = \underbrace{e^{-\omega^2 t}}_{\not= 0}
   \left( \vphantom{e^{-\omega^2 t}} A \cdot 1
          + B \cdot 0 \right)
   \quad \Rightarrow \quad A = 0.

So we simplified the form into
:math:`\displaystyle
U(x, t) = e^{-\omega^2 t}
\left( \vphantom{e^{-\omega^2 t}} B \sin(\omega x) \right)`.

To find values for :math:\omega`,
imposing :math:`u(1,t) = 0` leads to

.. math::

   U(1, t) = \underbrace{e^{-\omega^2 t}}_{\not= 0}
   \left( \vphantom{e^{-\omega^2 t}} B \sin(\omega) \right)
   \quad \Rightarrow \quad B = 0 \quad \mbox{or} \quad \sin(\omega) = 0.

Setting :math:`B` to zero sets :math:`U` to zero,
contradicting :math:`u(x, 0) = 1`, thus :math:`B \not= 0`.

.. math::

   \sin(\omega) = 0 \quad \Rightarrow \quad 
   \omega = n \pi, \quad n = 1,2, \ldots.

As :math:`B \not= 0`, for each multiple :math:`n` of :math:`\pi`,
we set :math:`B = B_n`.  So, we have

.. math::

   U_n(x, t) = e^{- n^2 \pi^2 t} 
   \left( \vphantom{e^{-\omega^2 t}} B_n \sin(n \pi x) \right),
   \quad n = 1,2, \ldots,

and thus infinitely many solutions to the heat equation
:math:`\displaystyle \frac{\partial u}{\partial t}
= \frac{\partial^2 u}{\partial x^2}`.

The operators :math:`\displaystyle \frac{\partial}{\partial t}`
and :math:`\displaystyle \frac{\partial^2}{\partial x^2}` are linear.
Therefore, consider

.. math::

   U_n(x, t) = e^{- n^2 \pi^2 t} 
   \left( \vphantom{e^{-\omega^2 t}} B_n \sin(n \pi x) \right),
   \quad n = 1,2, \ldots,

as a basis of functions.  Our new general form of a solution is then

.. math::

   U(x, t) = \sum_{n=1}^\infty B_n
   \left( e^{- n^2 \pi^2 t} \sin(n \pi x) \right).

The values :math:`B_n` are the coefficients of the solutions in the basis
of sine functions with decaying amplitudes.

Fourier Series
--------------

Imposing the initial condition :math:`u(x, 0) = 1` on

.. math::

   U(x, t) = \sum_{n=1}^\infty B_n
   \left( e^{- n^2 \pi^2 t} \sin(n \pi x) \right)

leads to

.. math::

   U(x, 0) = \sum_{n=1}^\infty B_n \sin(n \pi x) = 1.

To find the values for :math:`B_n`, 
we exploit the :index:`orthogonality` of the basis :math:`\sin(n \pi x)`,
:math:`n = 1, 2, \ldots`.

We compute :math:`B_n` in a :index:`Fourier series` approximation to

.. math::

   1 = \sum_{n=1}^\infty B_n \sin(n \pi x).

The inner product

.. math::

   \left\langle \vphantom{\frac{1}{2}}
   \sin(n \pi x), \sin(m \pi x) \right\rangle
   = \int_0^1 \sin(n \pi x) \sin(m \pi x) dx
   = \left\{
        \begin{array}{ccc}
           0 & \mbox{if} & n \not=m \\
           \frac{1}{2} & \mbox{if} & n = m
        \end{array}
     \right.

shows that sine functions running at different frequencies
are orthogonal to each other.
Applying the :index:`inner product` to

.. math::

   1 = \sum_{n=1}^\infty B_n \sin(n \pi x)

leads to

.. math::

   \left\langle \vphantom{\frac{1}{2}} 1, \sin(m \pi x) \right\rangle
   = \sum_{n=1}^\infty B_n 
   \left\langle \vphantom{\frac{1}{2}} \sin(n \pi x), \sin(m \pi x)
   \right\rangle 
   =  B_m \frac{1}{2}.

Evaluating the integral goes then as follows:

.. math::

   \frac{1}{2} B_m
   = \left\langle \vphantom{\frac{1}{2}} 1, \sin(m \pi x) \right\rangle
   = \int_0^1 \sin(m \pi x) dx
   = \frac{1}{m \pi}
   \left( \vphantom{\frac{1}{2}} 1 - \cos(m \pi) \right)

The cosine is :math:`+1` at even multiples of :math:`\pi`,
and :math:`-1` at odd multiples:

.. math::

   \cos(m \pi) = \left\{
    \begin{array}{ccl}
       -1 & \mbox{if} & m = 2 k+1, \\
       +1 & \mbox{if} & m = 2 k.
    \end{array}
   \right.

For even :math:`m`, :math:`B_m = 0`.
For odd :math:`m`, 
:math:`\displaystyle B_m = \frac{4}{m \pi}` or
:math:`\displaystyle B_{2k+1} = \frac{4}{(2k+1)\pi}`.

To conclude, the imposition of :math:`u(x, 0) = 1` on

.. math::

   U(x, t) = \sum_{n=1}^\infty B_n
   \left( e^{- n^2 \pi^2 t} \sin(n \pi x) \right)

has given us

.. math::

   B_{2k+1} = \frac{4}{(2k+1)\pi}, \quad k=0, 1, 2, \ldots

which leads to the general solution

.. math::

   U(x, t) = \sum_{k=0}^\infty \frac{4}{(2k+1)\pi}
   \left( e^{- (2k+1)^2 \pi^2 t} \sin((2k+1) \pi x) \right).

Looking at the general solution
we ask two convergence questions:

1. How many terms are needed in the sum?

2. How fast does the temperature go to zero?

These questions are answered by looking at plots.
We start with 10 terms in the sum, and then double
the number of terms till we arrive at 160.
The plots at time :math:`t = 0` are in 
:numref:`fig10terms`,
:numref:`fig20terms`,
:numref:`fig40terms`,
:numref:`fig80terms`, and
:numref:`fig160terms`.

.. _fig10terms:

.. figure:: ./fig10terms.png
   :align: center

   At :math:`t = 0`, with 10 terms in the sum.

.. _fig20terms:

.. figure:: ./fig20terms.png
   :align: center

   At :math:`t = 0`, with 20 terms in the sum.

.. _fig40terms:

.. figure:: ./fig40terms.png
   :align: center

   At :math:`t = 0`, with 40 terms in the sum.

.. _fig80terms:

.. figure:: ./fig80terms.png
   :align: center

   At :math:`t = 0`, with 80 terms in the sum.

.. _fig160terms:

.. figure:: ./fig160terms.png
   :align: center

   At :math:`t = 0`, with 80 terms in the sum.

From :numref:`fig160terms`, we see that the errors of 
the Fourier series approximation spikes at the edges,
but those spikes disappear as soon as :math:`t > 0`,
as can be observed from the evolution in time 
at :numref:`figHeatEvolution`.
We see that at time :math:`t = 0.5` the temperature
of the tube becomes very close to zero.

.. _figHeatEvolution:

.. figure:: ./figHeatEvolution.png
   :align: center

   The evolution in time with 160 terms in the sum.

Proposals for Topics of Projects
--------------------------------

1. the Fourier ring problem

   A long rod is bent and its ends welded to form a ring.
   Given an initial temperature profile around the ring,
   how will this temperature profile evolve over time?

   Study the solution of Problem 3 in our text book.

   Do the following:

   * Describe the solution using your own words.

   * Do the Fourier series computations.

   * Write a report on your computational experiments.

2. the heat diffusion problem in polar coordinates

   Consider the preparation for a camping trip.
   To not freeze to death during the night in a tent,
   a heater for tent camping is needed.
   Our tent has a circular base, with the heater at its center.
   Develop a mathematical model to prepare for a safe camping trip.

   Do the following:

   * Gather data on tent heaters.

   * Describe the heat diffusion equation in polar coordinates.

   * Interpret the computational results with the data.

Exercises
---------

1. Throughout this lecture we consider a cooling off process.

   Consider the heating of a tube filled with cold water.
   The tube is insulated everywhere, with heat provided at the ends.
   Answer the following:

   1. Formulate the boundary conditions for this model.

   2. What is the geometric interpretation?

2. The Laplacian with zero boundary conditions is negative definite:

   .. math::

      \mbox{For any } \phi \not=0: \langle \nabla^2 \phi, \phi \rangle < 0.

   Prove this statement.  
   *Hint:* See the outline of exercise 11.6 in the text book.

3. Instead of :math:`u(x, 0) = 1`, consider a bell shaped profile, with

   1. the center is at the middle :math:`1/2` 
      of the interval :math:`[0, 1] \ni x`, and

   2. the standard deviation of the Gaussian is :math:`1/2` as well.

   Do the following.

   1. Set up the model for the evolution of the temperature distribution.

   2. Compute the Fourier series of the profile.

   3. Provide plots to display the evolution.

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

1. Richard Haberman:
   *Applied Partial Differential Equations with
   Fourier Series and Boundary Value Problems*.  Pearson, 2013.

2. Charles R. MacCluer, problem 1 of section 11.3 of
   *Industrial Mathematics.  Modeling in Industry, Science, and Government*,
   Prentice Hall, 2000.