<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 17
Step1:  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.
Step2:  Set up other simulation parameters
We now set
Step3:  Our choice of timestep is restricted
Step4:  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.
Step5:  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
Step7:  DIY
Step8:  DIY
Step9:  DIY
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
Explanation: Lecture 17:  Numerical Solutions to the Diffusion Equation (Explicit Methods)
Reading and Reference
Numerical Recipes, W. Press, Cambridge University Press, 1986
Numerical Methods, R. Hornbeck, Quantum Publishers, 1975
B. Gustafsson, Fundamentals of Scientific Computing, Springer, 2011
S. Farlow, Partial Differential Equations for Scientists and Engineers, Dover, 1993 
What to Learn?
How to set up the finite difference equations
How to set up a data structure to organize the solution
How to apply the initial and boundary conditions to permit computation of the solution
What to do?
Write the solver code for three example problems
Visualize the results
Introduction
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 Series and Derivatives
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
$$
$$
h^{2} \left. \frac{d^{2}}{d \xi_{1}^{2}}  f{\left (\xi_{1} \right )} \right|{\substack{ \xi{1}=c }} + 2 f{\left (c \right )} - f{\left (c - h \right )} - f{\left (c + h \right )} = 0
$$
Dividing Space and Time
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.  
Finite Difference Form for the Diffusion Equation
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:
* A forward difference is used on the LHS
* A central difference is used on the RHS
The general steps required to solve a finite difference equation are as follows:
Identify the geometry of the system
Define the initial conditions and boundary conditions
Write the difference equation in terms of the unknowns
Compute the unknowns subject to any stability requirments and boundary conditions
Store and update the results in the grid
Write the results to disk
Visualize the results
The following parameters are needed to solve the finite difference equation above:
numberOfPoints - the number of grid points within our computational domain
numberOfIterations - the number of timesteps within our computational domain
lengthOfDomain - the physical length of our computational domain
dx - the distance between successive grid points in our domain
xPoints - a linear space divided into numberOfPoints points
initialCondition - our starting distribution of solute (i.e. $c(x,0)$)
End of explanation
numberOfPoints = 100
numberOfIterations = 1000
lengthOfDomain = 1.0
dx = lengthOfDomain/numberOfPoints
xPoints = np.linspace(0.0, lengthOfDomain, numberOfPoints)
initialCondition = np.sin(xPoints*np.pi/lengthOfDomain)
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
plotIC()
Explanation: 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.
End of explanation
diffusionCoefficient = 10.0
dt = dx**2/(4*diffusionCoefficient)
Explanation: Set up other simulation parameters
We now set:
diffusionCoefficient - the diffusion coefficient
dt - the discrete time step (the formulaic choice of dt is needed to satisfy stability constraints)
End of explanation
arrayWithAllTheData = np.zeros((numberOfPoints,numberOfIterations), dtype='float32')
Explanation: Our choice of timestep is restricted:
$$
dt \leq \frac{\Delta x^2}{2 \, D}
$$
Set up the data structure
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:
End of explanation
arrayWithAllTheData[:,0] = initialCondition
Explanation: 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.
End of explanation
# 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?
Explanation: 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:
Write the Solver and Compute the Results
Add your code below.  As a reminder, the algebraic equation we want to solve is:
$$
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}
$$
To make all of this work we proceed as follows:
Compute the prefactor $D \Delta t/ \Delta x^2$.
Apply the boundary conditions in the $j$th row of the array.
Using the $j$ row of the array (plus the boundary conditions), fill in row $j+1$ of the array with values corresponding to the new time $t + \Delta t$ according to the equation above.
Advance the index and repeat until all rows are filled.
Visualize the results.
End of explanation
%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
interact(plotArray, 
         xPoints=fixed(xPoints), 
         dataArray=fixed(arrayWithAllTheData),
         rowID=(0,numberOfIterations-1,1), );
Explanation: DIY:  Sketch the algorithm up to this point and for the cell below.
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.
Plot the results
End of explanation
# Your solver code goes here.
Explanation: DIY:  Compute a solution where you change the boundary conditions on the LHS to be $c(x=L,t) = 1.0$.
End of explanation
# Your solver code goes here.
Explanation: DIY:  Vectorize the above solution method.
End of explanation
<END_TASK>