<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Lecture 18
Step1:  We define the symbols we want to use in our linear system.  For this demonstration, I don't add the time index but I keep my subscripts consistent with the figure above.
Step2:  In this cell we create the square matrix holding the coefficients that multiply the unknown quantities.  Note the structure of the matrix.  It is a tridiagonal matrix.  The function in NumPy is very compact, in SymPy not so much.  So I apologize for the syntax in the SymPy/Python code below, but the more compact version can be difficult to read
Step3:  Here is our vector of unknown quantities.  We know the solution to the $j-1$ time step.  All of these symbols represent the value of our field (e.g. concentration, temperature, etc.) at the $j$'th time step.
Step4:  If we've got everything correct, this matrix product will reproduce the discrete diffusion equation outlined above.  You'll note that the boundary equations are not formed correctly.  For reference, here is the discrete form
Step5:  It should start to become clear that we can write this linear system (of a tridiagonal matrix and a column vector of unknowns) as a matrix equation
Step6:  So the full form of our system is therefore
Step7:  All that remains is to solve the above linear system.  Instead of using SymPy, we will use some tools in a different Python library.
Top of Page
An Implicit Numerical Solution
General setup in this section
Step8:  Simulation parameters
Step9:  A simple function to plot the initial condition
Step10:  It is worth noting that these schemes are unconditionally stable - so any choice of time step will produce a solution.  The accuracy of the solution does depend on this choice, though.
Step11:  We create two solution vectors rather than one whole array to hold all of our solution.  This is not particular to the implicit method, but it demonstrates another technique for saving memory and speeding up the calculation.  We will fill these matrices and swap them (move data from new into old and overwrite new) at each time step.
Step12:  First, some syntax
Step13:  The matrix has to be square.  It should have the same dimensions as the nubmer of points in the system.  The following code snippet was inspired by this post.
Step14:  We first need to prime the arrays by copying the initial condition into oldConcentration.  Afterwards it will be enough to swap pointers (a variable that points to a memory location).
Step15:  Top of Page
Deconstruction of the Solution Scheme
In spite of the small chunk of code a few cells below, there is a lot going on.  Let us dissect it.  In bullet points
<ASSISTANT_TASK:>
Python Code:
import sympy as sp
sp.init_session(quiet=True)
var('U_LHS U_RHS')
Explanation: Lecture 18: Numerical Solutions to the Diffusion Equation
(Implicit Methods)
Sections
Introduction
Learning Goals
On Your Own
In Class
Revisiting the Discrete Version of Fick's Law
A Linear System for Diffusion
An Implicit Numerical Solution
Deconstruction of the Solution Scheme
Homework
Summary
Looking Ahead
Reading Assignments and Practice
Possible future improvement:  Crank-Nicholson.
Introduction
This lecture introduces the implicit scheme for solving the diffusion equation. Examine the descritization for our explicit scheme that we covered in the previous lecture:
$$
\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}
$$
This expression uses a forward difference in time where we are subtracting the value of our dependent variable at time index $j$ from the value of our dependent variable at time-index $j+1$.  If, instead, we perform a backward difference (replace $j$ with $j-1$) we will be subtracting our dependent variable at $j-1$ from the value at the index $j$.  For example:
$$
\frac{u_{i,\, j} - u_{i,\, j-1}}{\Delta t} = D \frac{u_{i - 1,\, j} - 2 u_{i,\, j} + u_{i + 1,\, j}}{\Delta x^2}
$$
Attempting to repeat our previous algebraic manipulations we find that the solution to this equation is in terms of three unknown quantities at the index $j$.  These quantities depend on indices $i-1$, $i$ and $i+1$ and our solution is only known to the index $j-1$.  This seems like a complication that cannot be resolved however, examination of all the resulting equations in our grid will revel that this is a linear system that can be solved with the inclusion of the boundary conditions.
The point of this lecture is to develop your understanding for how the use of matrices and linear algebra can be used to solve this problem.  The system of equations and the method for solving the equations is known as an "implicit method".  There is a good discussion in Numerical Recipes by Teukolsky, et al. to provide a foundation for these methods.
Top of Page
Learning Goals
Re-develop the descretizaton of Fick's law such that the solution scheme is implicit
Write the method as a linear system
Incorporate boundary conditions
Develop a solution strategy using linear algebra and Numpy or SciPy as appropriate.
Top of Page
On Your Own
Suggestions for improvement of this section:
Develop numpy methods for linear algebra (e.g. creating matrices efficiently)
Matrix operations relevant to LA.
Solve a simple linear system.
Top of Page
In Class
Re-derive the discrete form of Fick's law.
Examine the structure of the resulting matrix.
Write a solver.
Revisiting the Discrete Version of Fick's Law
We start with a re-statement of Fick's second law in finite difference form that uses a FORWARD difference in time:
$$
\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}
$$
This choice of time differencing led to the EXPLICIT scheme.  This time, we choose a BACKWARD difference in time as follows:
$$
\frac{u_{i,\, j} - u_{i,\, j-1}}{\Delta t} = D \frac{u_{i - 1,\, j} - 2 u_{i,\, j} + u_{i + 1,\, j}}{\Delta x^2}
$$
This choice leads to a set of linear equations.  To illustrate how this develops we will write the equation above for a grid of eight points that represent the quantity of diffusing substance.  See the next figure.
In the above figure we represent a grid of two dimensions - this grid is identical to the explicit finite difference grid.  The main difference between the explicit and implicit method is the way we fill the grid to arrive at our solution.  In the spatial dimension we have 8 columns (the "$i$" index) and in the time dimension we show only three rows (the "$j$" index).  The sizes of the grid in the two dimensions are arbitrary.  Keep the following in mind:
The solution is known to the $j-1$ index.
The unknowns are the $j$ indcies.
Algebraiclly rearranging this differencing scheme, we can write down:
$$
u_{i,\, j-1} = \frac{\Delta t D}{\Delta x^2} \left( - u_{i - 1,\, j} + 2 u_{i,\, j} - u_{i + 1,\, j} \right) + u_{i,\, j}
$$
one additional re-arrangment (substitute $\beta$ for the factor containing the diffusion coefficient) and we get:
$$
- \beta u_{i - 1,\, j} + (1 + 2 \beta) u_{i,\, j} - \beta u_{i + 1,\, j} = u_{i,\, j-1}
$$
We include "ghost cells" in grey above to enforce the boundary conditions.  We can use fixed value (setting the ghost cells to a particular number) or fixed flux (setting the value of the ghost cell based on a pair of interior cells) to produce a linear system with an equal number of unknowns and equations.
A Linear System for Diffusion
We begin as usual by importing SymPy into the namespace and using init_session to define some helpful symbols.  We also define a pair of symbols $U_{LHS}$ and $U_{RHS}$ that will be used to define values in the ghost cells.
End of explanation
var('dt dx beta u1:7 b1:7')
Explanation: We define the symbols we want to use in our linear system.  For this demonstration, I don't add the time index but I keep my subscripts consistent with the figure above.
End of explanation
hpad = ones(0, 1); vpad = ones(1, 0)
mainDiag = 2*beta+1; offDiag = -beta
M = (sp.diag(vpad, offDiag, offDiag, offDiag, offDiag, offDiag, hpad)+ \
     sp.diag(hpad, offDiag, offDiag, offDiag, offDiag, offDiag, vpad)+ \
     sp.diag(mainDiag,mainDiag,mainDiag,mainDiag,mainDiag,mainDiag))
M
Explanation: In this cell we create the square matrix holding the coefficients that multiply the unknown quantities.  Note the structure of the matrix.  It is a tridiagonal matrix.  The function in NumPy is very compact, in SymPy not so much.  So I apologize for the syntax in the SymPy/Python code below, but the more compact version can be difficult to read:
End of explanation
xmatrix = sp.Matrix([u1,u2,u3,u4,u5,u6])
xmatrix
Explanation: Here is our vector of unknown quantities.  We know the solution to the $j-1$ time step.  All of these symbols represent the value of our field (e.g. concentration, temperature, etc.) at the $j$'th time step.
End of explanation
M*xmatrix
Explanation: If we've got everything correct, this matrix product will reproduce the discrete diffusion equation outlined above.  You'll note that the boundary equations are not formed correctly.  For reference, here is the discrete form:
$$
- \beta u_{i - 1,\, j} + (1 + 2 \beta) u_{i,\, j} - \beta u_{i + 1,\, j} = u_{i,\, j-1}
$$
End of explanation
bmatrix = sp.Matrix([(b1+beta*U_LHS),b2,b3,b4,b5,(b6+beta*U_RHS)]) 
bmatrix
Explanation: It should start to become clear that we can write this linear system (of a tridiagonal matrix and a column vector of unknowns) as a matrix equation:
$$
M \cdot \overline{x} = \overline{b}
$$
Where M is the square matrix, x is the vector of unknown quantities and b is the last known value of the system variables (the $u_{i,j}$ are the unknowns, the $j-1$ are the last known values).  There is still some work to be done before we can use linear algebra to get the solution.  We need to implement the boundary conditions.
Fixed Value Boundary Conditions
Start with the form at the interior of the grid:
$$
- \beta u_{i - 1,\, j} + (1 + 2 \beta) u_{i,\, j} - \beta u_{i + 1,\, j} = u_{i,\, j-1}
$$
To get the form correct at the top and bottom of this solution vector we need to imagine adding "ghost cells" to the boundaries of our domain at $i=0$ and $i=7$.  Using the above expression, let $i = 1$:
$$
- \beta u_{0,\, j} + (1 + 2 \beta) u_{1,\, j} - \beta u_{2,\, j} = u_{1,\, j-1}
$$
If we have fixed value boundary conditions, we then know the value of $u_0$.  This is the boundary condition of our simulation.  We will call this value $U_{LHS}$, substitute $U_{LHS} = u_0$ and move the known quantities to the RHS of the equation:
$$
(1 + 2 \beta) u_{1,\, j} - \beta u_{2,\, j} = u_{1,\, j-1} + \beta U_{LHS}
$$
Fixed Flux Boundary Conditions
If we have fixed flux boundary conditions we can write the flux as a central difference on the cell $u_1$ that uses the "ghost" point at $u_0$:
$$
\frac{u_{2,\, j} - u_{0,\, j}}{2 \Delta x} = F
$$
Proceeding as before with $i=1$:
$$
- \beta u_{0,\, j} + (1 + 2 \beta) u_{1,\, j} - \beta u_{2,\, j} = u_{1,\, j-1}
$$
This time we know the relationship of $u_0$ to the other unknowns due to the specification of the defined flux boundary condition.  Solving for $u_0$ we get:
$$
u_{0,\, j} = u_{2,\, j} - {2 \Delta x} F
$$
Substituting this into our expression that includes the ghost cell gives us:
$$
- \beta (u_{2,\, j} - {2 \Delta x} F) + (1 + 2 \beta) u_{1,\, j} - \beta u_{2,\, j} = u_{1,\, j-1}
$$
Simplifying:
$$
(1 + 2 \beta) u_{1,\, j} - 2 \beta u_{2,\, j} = u_{1,\, j-1} - \beta 2 \Delta x F 
$$
So in this case we have to modify the matrix $M$ entries AND the solution vector $b$ recalling that the $j-1$ index is the known solution.
We have now recovered the form of the equation in the dot product $M \cdot x$ and the form of this equation is telling us that we need to modify the solution vector $b$ with information about the boundary conditions before we find the inverse of the matrix and compute the new solution vector.
Modifying the $b$ matrix with the known ghost cell values for the fixed value boundary conditions we get:
End of explanation
sp.Eq(M*xmatrix,bmatrix)
Explanation: So the full form of our system is therefore:
$$
\left[\begin{matrix}2 \beta + 1 & - \beta & 0 & 0 & 0 & 0\- \beta & 2 \beta + 1 & - \beta & 0 & 0 & 0\0 & - \beta & 2 \beta + 1 & - \beta & 0 & 0\0 & 0 & - \beta & 2 \beta + 1 & - \beta & 0\0 & 0 & 0 & - \beta & 2 \beta + 1 & - \beta\0 & 0 & 0 & 0 & - \beta & 2 \beta + 1\end{matrix}\right] \cdot \left[\begin{matrix}u_{1}\u_{2}\u_{3}\u_{4}\u_{5}\u_{6}\end{matrix}\right] = \left[\begin{matrix}U_{LHS} \beta + b_{1}\b_{2}\b_{3}\b_{4}\b_{5}\U_{RHS} \beta + b_{6}\end{matrix}\right]
$$
SymPy can evaluate the LHS for us.
End of explanation
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Explanation: All that remains is to solve the above linear system.  Instead of using SymPy, we will use some tools in a different Python library.
Top of Page
An Implicit Numerical Solution
General setup in this section:
End of explanation
numberOfPoints = 100
lengthOfDomain = 1.0
dx = lengthOfDomain/numberOfPoints
xPoints = np.linspace(0.0, lengthOfDomain, numberOfPoints)
initialCondition = np.sin(xPoints*np.pi/lengthOfDomain)
Explanation: Simulation parameters:
End of explanation
def plotIC():
    fig = plt.figure()
    axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
    axes.plot(xPoints, initialCondition, 'ro')
    axes.set_xlabel('Distance $x$')
    axes.set_ylabel('Concentration of Stuff $c(x,t)$')
    axes.set_title('Initial Conditions');
plotIC()
Explanation: A simple function to plot the initial condition:
End of explanation
diffusionCoefficient = 10.0
dt = dx**2/(diffusionCoefficient)
numberOfIterations = 1000
Explanation: It is worth noting that these schemes are unconditionally stable - so any choice of time step will produce a solution.  The accuracy of the solution does depend on this choice, though.
End of explanation
newConcentration = np.zeros((numberOfPoints), dtype='float32')
oldConcentration = np.zeros((numberOfPoints), dtype='float32')
Explanation: We create two solution vectors rather than one whole array to hold all of our solution.  This is not particular to the implicit method, but it demonstrates another technique for saving memory and speeding up the calculation.  We will fill these matrices and swap them (move data from new into old and overwrite new) at each time step.
End of explanation
['h','h','h']*3
Explanation: First, some syntax:
End of explanation
def tridiag(a, b, c, k1=-1, k2=0, k3=1):
    # Here we use Numpy addition to make the job easier.
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)
a = [-dt*diffusionCoefficient/dx/dx]*(numberOfPoints-1)
b = [2*dt*diffusionCoefficient/dx/dx+1]*(numberOfPoints)
c = [-dt*diffusionCoefficient/dx/dx]*(numberOfPoints-1)
A = tridiag(a, b, c)
A
Explanation: The matrix has to be square.  It should have the same dimensions as the nubmer of points in the system.  The following code snippet was inspired by this post.
End of explanation
np.copyto(oldConcentration,initialCondition)
Explanation: We first need to prime the arrays by copying the initial condition into oldConcentration.  Afterwards it will be enough to swap pointers (a variable that points to a memory location).
End of explanation
uLHS = 0.0
uRHS = 0.0
numIterations = 200
for i in range(numIterations):
    # enforce boundary conditions
    oldConcentration[0] = oldConcentration[0] + uLHS*dt*diffusionCoefficient/dx/dx
    oldConcentration[-1] = oldConcentration[-1] + uRHS*dt*diffusionCoefficient/dx/dx
    # solve the system
    np.copyto(newConcentration,np.linalg.solve(A,oldConcentration))
    # swap pointers
    oldConcentration, newConcentration = newConcentration, oldConcentration
# plot the results
fig2 = plt.figure()
axes = fig2.add_axes([0.1, 0.1, 0.8, 0.8])
axes.plot(xPoints, newConcentration, 'ro')
axes.set_ylim(0,1)
axes.set_xlabel('Distance $x$')
axes.set_ylabel('Concentration of Stuff $c(x,t)$')
axes.set_title('Solution');
Explanation: Top of Page
Deconstruction of the Solution Scheme
In spite of the small chunk of code a few cells below, there is a lot going on.  Let us dissect it.  In bullet points:
Before the first solution step we enforce the boundary conditions.  Our choice of matrix means that we are using "fixed value" boundary conditions.  So we need to modify the b vector accordingly.  The indexing notation of Numpy that permits us to find the first ([0]) and last cell ([-1]) of an array is very helpful here.
```python
    oldConcentration[0] = oldConcentration[0] + uLHSdtdiffusionCoefficient/dx/dx
    oldConcentration[-1] = oldConcentration[-1] + uRHSdtdiffusionCoefficient/dx/dx
```
Recall:
$$
\left[\begin{matrix}2 \beta + 1 & - \beta & 0 & 0 & 0 & 0\- \beta & 2 \beta + 1 & - \beta & 0 & 0 & 0\0 & - \beta & 2 \beta + 1 & - \beta & 0 & 0\0 & 0 & - \beta & 2 \beta + 1 & - \beta & 0\0 & 0 & 0 & - \beta & 2 \beta + 1 & - \beta\0 & 0 & 0 & 0 & - \beta & 2 \beta + 1\end{matrix}\right] \cdot \left[\begin{matrix}u_{1}\u_{2}\u_{3}\u_{4}\u_{5}\u_{6}\end{matrix}\right] = \left[\begin{matrix}U_{LHS} \beta + b_{1}\b_{2}\b_{3}\b_{4}\b_{5}\U_{RHS} \beta + b_{6}\end{matrix}\right]
$$
Solving the system involves using the built in NumPy functions to invert the matrix.  What is returned is the solution vector.  Please note that I'm using an internal Numpy (an optimized function!) function to COPY the results of the linear algebra solution into the newConcentration vector. 
python
    np.copyto(newConcentration,np.linalg.solve(A,oldConcentration))
Rather than storing ALL the data, we instead store just the current and the old concentrations.  There are efficiencies in doing this, but if we want the older values, we need to store them on disk or in memory.
Tuple unpacking in Python leads to the A,B=B,A syntax below.  This switches the references to the arrays.  This is important for efficiency - you don't want to move any data if you don't have to.  If you are running big calculations then moving that data around is a waste of time/resources.  Better to just swap references.
python
    oldConcentration, newConcentration = newConcentration, oldConcentration
Repeat the process and after a specified number of iterations, plot the results.
End of explanation
<END_TASK>