In each of the next three lectures a different numerical technique will be studied and implemented to solve the diffusion equation. Each technique is built from mathematics that we've studied in previous lectures. These or similar methods are used in numerical software and an in-depth understanding of the basic presented will give the student a foundation for more advanced methods that may be encountered in commercial codes.
The first technique we will study is the Explicit Finite Difference method. This is one of three common finite difference methods that are easy to program and built from the definition of Taylor's polynomial.
Taylor's approximation of a first and second derivatives are defined as central differences and are second order accurate:
$$ - 2 h \left. \frac{d}{d \xi_{1}} f{\left (\xi_{1} \right )} \right|_{\substack{ \xi_{1}=c }} - f{\left (c - h \right )} + f{\left (c + h \right )} = 0 $$The diffusion equation is a partial differential equation in two independent variables, space and time. In order that we may use our Taylor's series approximations for the time and space derivatives we need to descretize the domain of the problem. One easy way to visualize time and space descretization is to use a grid construction. The figure below shows one way that the time and space variables can be represented.
Each cell in the grid holds the value of the dependent parameter at a particular value of time and space. We indicate this using the following symbols:
$$ u_{i,j} $$The $i$ and $j$ are the indices of the grid and reference a particular location where the value of the dependent parameter is stored. How much memory is required to store these values depends on the type of data and the size of the grid. In addition, the grid spacing must be specified for each of the independent variables, in this case we need both a $\delta x$ and a $\delta t$. For example, a difference in time might be represented as:
$$ u(i,j) - u(i,j+1) = c(x,t) - c(x,t+\delta t) $$Typically, the grid is a uniform grid in each independent variable - meaning that the distance in the independent variable between grid points is the same in any one variable. There are cases where non-uniform grids may be desirable.
The finite difference representation of the diffusion equation with $u_{i,j}$ as the stored value that represents $c(x,t)$ is:
$$ \frac{u_{i,\, j+1} - u_{i,\, j}}{\Delta t} = D \frac{u_{i - 1,\, j} - 2 u_{i,\, j} + u_{i + 1,\, j}}{\Delta x^2} $$Note:
The general steps required to solve a finite difference equation are as follows:
The following parameters are needed to solve the finite difference equation above:
numberOfPoints - the number of grid points within our computational domainnumberOfIterations - the number of timesteps within our computational domainlengthOfDomain - the physical length of our computational domaindx - the distance between successive grid points in our domainxPoints - a linear space divided into numberOfPoints pointsinitialCondition - our starting distribution of solute (i.e. $c(x,0)$)
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
In addition to the simulation parameters, we start with an initial seed of concentration data. Unlike our other analytical strategies there are no coefficients to compute, no functions to fit. The data could be generated from a function or from measurements. Here we choose a sin function as our initial condition.
In [ ]:
numberOfPoints = 100
numberOfIterations = 1000
lengthOfDomain = 1.0
dx = lengthOfDomain/numberOfPoints
xPoints = np.linspace(0.0, lengthOfDomain, numberOfPoints)
initialCondition = np.sin(xPoints*np.pi/lengthOfDomain)
In [ ]:
def plotIC():
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(xPoints, initialCondition, 'ro')
ax.set_xlabel(r'Distance $x$')
ax.set_ylabel(r'Concentration $c(x,t)$')
ax.set_title(r'Initial Conditions')
plt.show()
return
In [ ]:
plotIC()
We now set:
diffusionCoefficient - the diffusion coefficientdt - the discrete time step (the formulaic choice of dt is needed to satisfy stability constraints)
In [ ]:
diffusionCoefficient = 10.0
dt = dx**2/(4*diffusionCoefficient)
Our choice of timestep is restricted:
$$ dt \leq \frac{\Delta x^2}{2 \, D} $$There are potentially three strategies for storing the results of numerical computations. One is to store ALL the data, another is to store SOME of the data, and the last is to store NONE of the data except the very last computation. Each strategy has advantages and disadvantages although all strategies may seem equally difficult to implement. In this lecture we will design design a data structure that stores all the data.
Let us create a numpy array that has one dimension equal to the number of points in our grid and another dimension that is equal to the number of iterations we wish to compute:
In [ ]:
arrayWithAllTheData = np.zeros((numberOfPoints,numberOfIterations), dtype='float32')
You can think of the 2D array as having one space axis (the first index, we will use i for this one) and one time axis (the second index, we will use j for this one).
We will set our initial conditions by assigning the initialCondition array to the first row of the arrayWithAllTheData. Note the slicing in the first index so that we can copy the contents of the initialCondition into the whole first row with a single assignment statement.
In [ ]:
arrayWithAllTheData[:,0] = initialCondition
With our initial conditions in place we need to develop the computational steps to advance our solution in time. The PDE we are solving (with a constant diffusion coefficient) is:
$$ \frac{\partial c(x,t)}{\partial t} = D \frac{\partial^2 c(x,t)}{\partial x^2} $$we transform this into:
$$ \frac{u_{i,\, j+1} - u_{i,\, j}}{\Delta t} = D \frac{u_{i - 1,\, j} - 2 u_{i,\, j} + u_{i + 1,\, j}}{\Delta x^2} $$so that we can algebraically solve for $u_{i+1,\, j}$:
$$ u_{i,\, j+1} = \frac{D \Delta t}{\Delta x^2} \left( u_{i - 1,\, j} - 2 u_{i,\, j} + u_{i + 1,\, j} \right) + u_{i,\, j} $$From the expression above you can see that all the terms on the RHS of the expression are at the index $j$ (the last iteration) and all the terms on the LHS are for the $j+1$ index (the next iteration). This scheme defines a simple method (with a restrictive timestep) for integrating a PDE. Re-examine the figure below in comparison to the finite difference scheme:
To make all of this work we proceed as follows:
In [ ]:
# Note the counting for j in this loop. You may wish to print out
# the values of i,j to help build this operation.
for j in range(1,numberOfIterations):
for i in range(1,numberOfPoints-1):
arrayWithAllTheData[i,j] = 0 # What should you put here?
Doing this will help you visualize the operations and it will increase your ability to make modifications in the future and devise new more compact ways to integrate this PDE.
If you've sketched the algorithm as advised above then you see that in our development of this solution we implicitly set the boundary conditions. We initialize arrayWithAllTheData with np.zeros and then compute on all the interior rows/columns. This creates a condition where all the boundary cells are set to zero and their values remain untouched throughout the computation.
In [ ]:
%matplotlib inline
import numpy as np
from ipywidgets import interact, fixed
import matplotlib.pyplot as plt
def plotArray(xPoints, dataArray, rowID=0):
"""
This function in conjunction with interact() permits
inspection of the contents of an array, row by row. This
is useful for some small tasks such as examining the results
of a PDE solution.
"""
x = xPoints
y = dataArray[:,rowID]
fig = plt.figure(figsize=(7,4))
axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
axes.set_ylim(0,1)
axes.plot(x, y, 'ro', label=r"$c(x,t)$")
axes.legend()
axes.grid(False)
plt.show()
return
In [ ]:
interact(plotArray,
xPoints=fixed(xPoints),
dataArray=fixed(arrayWithAllTheData),
rowID=(0,numberOfIterations-1,1), );
In [ ]:
# Your solver code goes here.
In [ ]:
# Your solver code goes here.
Problem 1. Solve the following problem:
$$ \frac{\partial c(x,t)}{\partial t} = D \frac{\partial^2 c(x,t)}{\partial x^2} $$with the initial condition
$$c(x,t=0) = \sin \pi x,$$over the domain
$$( x \, \vert \, 0 \le x \le 1.0 ).$$with zero flux boundary conditions.
Problem 2. Solve the following problem:
$$ \frac{\partial c(x,t)}{\partial t} = \frac{\partial}{\partial x} D(c) \frac{\partial c(x,t)}{\partial x} $$with the initial condition:
$$ c(x \leq 0.5 , t=0) = 1.0\\ c(x \gt 0.5 , t=0) = 0.0 $$with D having the dependence:
$$ D(c) = D_0\cdot c \cdot(1-c) $$over the domain:
$$ ( x \, \vert \, 0 \le x \le 1.0 ) $$with zero flux boundary conditions.