# Lecture 15:  Solution to the Diffusion Equation (Separation of Variables)

## Reading and Reference

* Essential Mathematical Methods for Physicists, H. Weber and G. Arfken, Academic Press, 2003
* Advanced engineering Mathematics, E. Kreyszig, John wiley and Sons, 2010
* S. Farlow, Partial Differential Equations for Scientists and Engineers, Dover, 1993

## What to learn?

* See the solution to a partial differential equation (PDE) by reducing it to two ordinary differential equations (ODE).
* Determine constants by applying the boundary conditions.
* Develop an initial condition that is consistent with the physical problem.
* Visualize the results.

## What to do?

* Solve the diffusion PDE using separation of variables.
* Practice using `dsolve` and applying boundary conditions using substitutions.
* Develop code to help you visualize solutions.

### Introduction
----

Separation of variables is a technique that reduces a partial differential equation to a pair of ordinary differential equations.  Solutions of the resulting ODEs are more straightforward than for the PDE.  This lecture and computational notebook will help you analyze one such solution to the diffusion equation.  To start we will select a problem that is best solved with a Fourier `sin` series.  Making this selection at the start is arbitrary, however the consequences will define the physics of the problem.  This is most easily seen in the restrictions on the boundary conditions.  For this problem our boundary conditions will be chosen to be compatible with a Fourier `sin` series.


```python
%matplotlib notebook
import sympy as sp
# You may have to comment out the init_sesison() line if you have not used my session.py file.
sp.init_session(quiet=True)
sp.init_printing()
dum = sp.symbols('dum')
```


```python
def initial_condition(x):
    return x

def b_m_amplitudes(m, funToProject, center=0.5, width=1.0):
    return (2/width)*sp.integrate(funToProject(dum)*sp.sin(2*m*sp.pi*dum/width), (dum,center-width/2,center+width/2))

def b_m_vectorspace_element(m, var, width=1.0):
    return sp.sin(2*m*sp.pi*var/width)

terms = 10
width = 2
center = 0

sinSeries = sum([b_m_amplitudes(m, initial_condition, center=center, width=width)*b_m_vectorspace_element(m,x,width=width) for m in range(terms)])
```


```python
sp.plot(sinSeries,(x,0,1))
```

The physical problem of interest is one where an extensive quantity (energy, mass, etc.) evolves diffusivly in some medium.  Using mass diffusion as an example, we will identify the concentration of species in our medium at any point and tiime, $c(x,t)$.  We must choose a domain (the geometry of the medium) in which this diffusion process will take place.  In the solution developed below we assume that we have a long bar of material with a constant cross section and this bar is insulated so that none of our extensive quantity can escape through the ends or through the cylindrical surface.  These assumptions permit simplification of the problem to one dimension.


We will assume that our bar of material has a fixed length $L = 1.0$ units long.  The medium will contain an initial amount of the diffusing substance whose distribution will be given by a function that represents the initial condition.  We will represent our initial condition using the symbol $\phi(x,0)$ where the zero indicates time $t=0$.  For this exercise we will choose:

$$
\phi(x,0) = x
$$

The solution to the problem depends in part on our ability to represent the initial condition as a Fourier series.  The choice of $\phi(x)=x$ permits us demonstrate something about the series solution at the boundary.  Let us visualize the solution to get a feel for what we're dealing with first and then we will develop the solution ourselves later in this notebook.  In the `interact` function we can change the time, the number of terms, and the diffusion constant ($\alpha$).  

* We define a function called `solutionAtX` that takes a pair of values ($x$, $t$), a diffusion coefficient $\alpha$ and the number of terms to use in the Fourier series.
* We define another function to produce a plot (`make_plot`) that we pass to `interact` so that we have an interactive plot to examine the solution and a set of parameters.


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact

def solutionAtX(x, t, alpha, terms):
    # This is an already developed solution to the diffusion equation.
    # If you can follow along in the lecture - your solution will look
    # identical.
    return sum([
            (2/(np.pi*n)**2)*
            (np.sin(n*np.pi*x))*
            (-np.pi*n*np.cos(np.pi*n)+np.sin(np.pi*n))*
            np.exp(-t*(alpha*n*np.pi)**2)
         for n in range(1,terms+1,1)
        ])

def make_plot(t=0.01, alpha=0.2, terms=20):
    plt.figure(figsize=(7,5))
    x = np.arange(0, 1, 0.001)
    y = solutionAtX(x, t, alpha, terms)
    plt.plot(x, y, x, x, 'orange')
    plt.fill_between(x, 0, y, alpha=0.2)
    plt.ylim(0.0, 1.0)
    plt.xlabel('Distance')
    plt.ylabel('Concentration')
    plt.show()
```


```python
interact(make_plot, t=(0.001,2,0.01), alpha=(0.2,0.5,0.01), terms=(1,40,1));
```

Using the interactive widgets above, explore different combinations of parameters and connect your physical intuition for diffusion with the solution.  The orange line is the line $\phi(x,0)=x$ and the blue line is the solution to the diffusion equation.  As you advance time try and connect your physical intuition of diffusion with the evolution of the area under the blue curve in the solution.  At the boundaries the diffusing species is allowed to escape the medium.  Over time, we should see a gradual decrease in the amount of substance (shaded area under the blue curve) in the solution to the equation.

### 1D Diffusion Problem - Separating Variables
----

There are many ways to solve a PDE.  Separating variables is one such method to reduce a PDE to a pair (or more) of ODEs.  We start by assuming that the solution has the form of a time dependent piece and a spatially dependent piece.  So we could write:

$$ c(x,t) = X(x) T(t) $$

We substitute the assumed form of the solution into Fick's law:


$$
\frac{\partial c(x,t)}{\partial t} = D \frac{\partial^2 c(x,t)}{\partial x^2}
$$

this gives us:

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


We will differentate the assumed solution of $c(x,t)$ above according to the derivatives in the differential equation.  It is usual to let D be $\alpha^2$ to guarantee that D is always positive, we can do this in `sympy` with a keyword in the `symbols` function call.


```python
import sympy as sp
sp.init_session(quiet=True)
sp.init_printing()
```


```python
var('X,T, c',cls=Function)
```


```python
var('alpha gamma', positive=True)
```

With our substitutions in place:

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

After taking the derivatives we have:

$$ X(x)T'(t) = \alpha^2 X''(x)T(t) $$

We can arrange these expressions so that all the functions of "t" are on one side and the functions of "x" are on the other side.  Divide both the LHS and RHS of the above expression by $\alpha^2 X(x) T(t)$ to seperate variables on either sides of the equals sign.  Note also that this means that the product $\alpha^2 X(x) T(t)$ should never be zero.  Add a seperation constant $\gamma$ and integrate.

The rationalization for the seperation constant is that since the LHS and RHS each depend on one variable then the only way the LHS and RHS can be equal is if they are both equal to a constant.  Our ODEs are then written as:


```python
timeODE = sp.Eq(T(t).diff(t,1)/(alpha**2*T(t)),-gamma**2)
timeODE
```


```python
spaceODE = sp.Eq(X(x).diff(x,2)/X(x),-gamma**2)
spaceODE
```

Invoke `dsolve()`


```python
sp.dsolve(timeODE,T(t))
```

This is not the form you are likely to find in a textbook.  It is the correct solution, but it may not be the best solution for us.  We can examine all the possible forms of the solution so that we can inspect the list and pick the one we like.  We create a list of hints and then use a list comprehension to solve the equation for each hint.  The notebook displays all the solutions so that we can find the one that suits our needs. 


```python
list_of_hints = sp.classify_ode(timeODE)
list_of_solutions = [sp.dsolve(timeODE,T(t), hint=hint) for hint in list_of_hints]
list_of_solutions
```

 A bit of matching up with the list of hints and we find that the one we are interested in is `1st_linear`:


```python
timePiece = sp.dsolve(timeODE,T(t), hint='1st_linear')
timePiece
```


```python
spacePiece = sp.dsolve(spaceODE,X(x))
spacePiece
```

The solution (with constants) is the product of the two solutions above:

$$ 
c(x,t) = X(x)T(t) = \Big( C_{1} \sin{\left (\gamma x \right )} + C_{2} \cos{\left (\gamma x \right )} \Big)\Big ( C_{3} e^{- \alpha^{2} \gamma^{2} t} \Big ) 
$$

We will need to "absorb" the time solution constant ($C_3$) into the space solution constants.  There are an infinity of solutions depending on our choices of $C_1$ and $C_2$.  These will ultimately be determined by the boundary conditions.

I would also note here the following:  the CAS (SymPy in this case) did all the work for us above.  The code below is really to just keep ourselves organized.  I don't intend to do much with `generalSolution` except to make substitutions to find the boundary conditions.  Working with a CAS is equal parts setup, organization, and deciphering the results.

Clean up the solution for use below.  Eliminate one of the constants as the product of C1 and C2 is also a constant.  Sympy does not automatically create the symbols corresponding to the constants for us.  So we need to create those before proceeding:


```python
var('C1 C2')
```


```python
generalSolution = Eq(c(x,t),(timePiece.rhs).subs([(C1,1)])*spacePiece.rhs)
generalSolution
```

### Boundary Conditions
----

All that is left to do is impose the boundary conditions to find the constants $C_1$ and $C_2$ and get the particular solution that satisfies our initial condition.

The following conditions apply to this problem:

$$ c(0,t) = 0 $$
$$ c(1,t) = 0 $$
$$ c(x,0) = x $$

These conditions apply for the left-hand boundary, the right-hand boundary and the initial condition (the time boundary).  The domain of the experiment is $\{\,x\, \mid 0 < x < 1\,\}$


```python
generalSolution
```

Start by substituting the boundary condition $c(0,t) = 0$:


```python
generalSolution.subs({x:0})
```

The only way this can be true is if $C_2$ is zero.  Make that substition and insert the second boundary condition where $c(1,t)$ must also be zero:


```python
generalSolution.subs({x:1, C2:0})
```

_So - wait, what?  Can can we let $C1=0$ too?_

That IS a solution.  It just is not the one we're looking for, so we need to examine how the third factor (the sin term) can help us.  

What about the exponential?  At the boundary there is no dependence of $x$ in this factor however we know that the parameters in the exponential will not cause that factor to go to zero.  It starts at one and just gets smaller - so looking at the second factor will not help us much.  

This means that we have to let $\sin(\gamma)=0$.  We know this can happen with the right choice of $\gamma$.  What choices can we use?

What if:

$$\gamma = \pm \pi, \pm 2\pi, \pm 3\pi, ... $$

and in general:

$$\gamma = \pm n\pi$$

Let us go back an examine the initial condition where $t=0$.  


```python
generalSolution.subs({t:0, C2:0, gamma:n*sp.pi})
```

Now this statement tells us that we need to find a constant $C_1$ such that $c(x,0)=x$.  You may be wondering how we can do that - but you already know how to do that.

We identify the diffusion PDE (Fick's Law) as a homogeneous PDE. Kreyszig defines homogeneous as an equation in $u$ where "each of its terms contains either $u$ or one of its partial derivatives".  Because the equation is a homogeneous PDE we can invoke the principle of superposition that combinations of solutions to the PDE are also solutions to the PDE.

So it is therefore possible to write (without proof) that the initial condition:

$$c(x,0) = \sum_{n=1}^N C_n \sin(n\pi x) = \phi(x) = x$$

and the general solution to the equation is:


$$c(x,t) = \sum_{n=1}^N C_n \sin(n\pi x) e^{- \alpha^2 n^2 \pi^2 t}$$

The summation is the Fourier sine series - so as long as you can compute the coefficients $C_n$ then you can write down a summation that will be an approximate solution to the problem.  All that remains is to find the Fourier coefficient for the function $\phi(x)$ and substitute those into the general solution.  Go back and look at the plotted solution and see if you can write your own function derived from the general solution and a Fourier series.

Your homework assignment is to set up a different initial condition and solve the problem again.

### Homework
----

From Farlow's book entitled "Partial Differential Equations for Scientists and Engineers" he poses the following two problems:

Find the Fourier sine expansion of $\phi(x)=1$ on the domain $0 \le x \le 1$.  Structure your solution/code so that you can investigate the series as a function of the number of terms it contains.  It will be helpful to pay attention to the values of the wavenumber and coefficients.  (Maybe symbolic computation will help?)  By choosing hte `sin` series you will also be requiring the values of the function to be zero at the boundaries.  Hint:  The domain of the series may be larger than the domain of the physical problem.

Using the results of the previous problem, what is the solution to the following initial boundary value problem (PDE):

$$
\begin{array}{lll}
\mathrm{PDE} & u_t = u_{xx} & 0 \le x \le 1 \\
\mathrm{BCs} & u(0,t) = 0 & 0 \le t \le \infty \\
 & u(1,t) = 0 & \\
\mathrm{IC}  & u(x,0) = 1 & 0 \le x \le 1
\end{array}
$$

Note that there is a physical impossibility in this problem.  Can you identify it and discuss as part of your solution?


```python

```

00f8cee798949222790a797a8bc9ae3fc5e66244
ipynb
Jupyter Notebook
Lecture-15-Separation-of-Variables.ipynb
mathinmse/mathinmse.github.io
837e508bfeeb7d108019fb9bc499066b2b653551
MIT
2017-07-19T04:04:38.000Z
2022-02-18T19:33:43.000Z
Lecture-15-Separation-of-Variables.ipynb
mathinmse/mathinmse.github.io
837e508bfeeb7d108019fb9bc499066b2b653551
MIT
2019-04-08T15:21:45.000Z
2020-03-03T20:19:00.000Z
Lecture-15-Separation-of-Variables.ipynb
mathinmse/mathinmse.github.io
837e508bfeeb7d108019fb9bc499066b2b653551
MIT
2017-07-27T02:27:49.000Z
2022-01-27T08:16:40.000Z
Qwen/Qwen-72B
1. YES
2. YES