###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth.

# Spreading out

Welcome to the second lesson of the course module: "_Spreading out: parabolic PDEs."_ We're studying the heat equation in one spatial dimension:

$$
\begin{equation}
\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}
\end{equation}
$$

where $\alpha$ is the thermal diffusivity and $T$ is the temperature.

In the previous lesson, we reviewed the numerical solution of the 1D diffusion equation with a forward-time, centered-space scheme: an _explicit_ scheme. What does that mean?  

The solution for $T$ at time step $t^{n+1}$ was calculated using different combinations of $T$ values from the *previous* time step $t^n$.  We have complete knowledge of the parts that feed into the solution update at each spatial point.  

*Implicit* methods work differently: we will use more data from the "future" in the update, including several values of $T$ at $t^{n+1}$. This will make the scheme more difficult to apply, but there are several reasons why it may be worth the effort.

In lesson 1, we discussed two disadvantages of explicit methods: (1) boundary effects drag behind by one time step; (2) stability requirements constrain the time-step to very small values. Both of these issues are resolved by implicit schemes.

## Implicit schemes

Let's move things around a bit and try combining the Euler time step with an evaluation of the spatial derivative on the *updated* solution at $t^{n+1}$. The discretized form of the equation is now as follows:

$$
\begin{equation}
\frac{T_{i}^{n+1}-T_{i}^{n}}{\Delta t}=\alpha\frac{T_{i+1}^{n+1}-2T_{i}^{n+1}+T_{i-1}^{n+1}}{\Delta x^2}
\end{equation}
$$

The stencil for this discretization doesn't look anything like the other stencils we've used until now. Check it out.


.
#### Figure 1. Stencil of the implicit central scheme.

From the previous time-step, we only know $T_i^{n}$, but what about $T_i^{n+1}$, $T_{i-1}^{n+1}$ and $T_{i+1}^{n+1}$?  What can we do?

No need to panic! Let's start by putting what we *do know* on the right-hand side of the equation and what we *don't know* on the left. We get:

$$
\begin{equation}
-T_{i-1}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{i}^{n+1} - T_{i+1}^{n+1} = T_{i}^{n}\frac{\Delta x^2}{\alpha\Delta t}
\end{equation}
$$

It looks like there are a lot of unknowns and just one equation! 

What does it look like with $i=1$?

$$
\begin{equation}
-T_{0}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{1}^{n+1} - T_{2}^{n+1} = T_{1}^{n}\frac{\Delta x^2}{\alpha\Delta t}
\end{equation}
$$

and $i=2$?

$$
\begin{equation}
-T_{1}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{2}^{n+1} - T_{3}^{n+1} = T_{2}^{n}\frac{\Delta x^2}{\alpha\Delta t}
\end{equation}
$$

What about $i=3$?

$$
\begin{equation}
-T_{2}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{3}^{n+1} - T_{4}^{n+1} = T_{3}^{n}\frac{\Delta x^2}{\alpha\Delta t}
\end{equation}
$$

Can you see the common element across equations? Here's a little help:



$T_{i}^{n+1}$ also appears in the equation for $T_{i-1}^{n+1}$ and $T_{i+1}^{n+1}$. We might have enough equations if we apply this for all $i$-values at the same time, don't you think? In fact, this is a linear system of equations for the unknown values $T_{i}^{n+1}$ on the spatial grid.

### What about boundary conditions? 

Let's look at the boundary points of the example from the previous lesson with a Dirichlet BC at $x=0$ and a Neumann BC at $x=1$, discretizing with $N$ mesh points.

The value $T_0^{n+1}$ is known at every time-step from the BC, so putting all unknown terms on the left-hand side of the equation and the known values on the right side yields the following for the $i=1$ equation:

$$
\begin{equation}
-T_{2}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{1}^{n+1} = T_{1}^{n}\frac{\Delta x^2}{\alpha\Delta t} + T_{0}^{n+1}
\end{equation}
$$

That was easy!

On the other hand, for $i=N-2$, the equation reads

$$
\begin{equation}
-T_{N-3}^{n+1} + \left( 2 + \frac{\Delta x^2}{\alpha\Delta t}\right) T_{N-2}^{n+1} - T_{N-1}^{n+1} = T_{N-2}^{n}\frac{\Delta x^2}{\alpha\Delta t}
\end{equation}
$$

The discretized Neumann boundary condition on the right side of the rod is

$$
\begin{equation}
\frac{T^{n}_{N-1} - T^{n}_{N-2}}{\Delta x} = q
\end{equation}
$$

But we can just as easily write that at time step $n+1$ (the boundary conditions apply at every time-step):

$$
\begin{equation}
\frac{T^{n+1}_{N-1} - T^{n+1}_{N-2}}{\Delta x} = q
\end{equation}
$$

Inserting the Neumann boundary condition in the equation for $i=N-2$ yields

$$
\begin{equation}
    -T_{N-3}^{n+1} + \left( 1 + \frac{\Delta x^2}{\alpha\Delta t} \right) T_{N-2}^{n+1} = T_{N-2}^{n} \frac{\Delta x^2}{\alpha\Delta t} + \Delta x q
\end{equation}
$$

Make sure you work this out with pen and paper: it's important to recognize where these terms come from!

Now we can write the linear system of equations in matrix form as follows:

$$
[A][x] = [b]+[b]_{b.c.}
$$

where the matrix of coefficients $[A]$ is a sparse matrix—most of the matrix elements are zero—with three non-zero diagonals. We write below the system expanded out, so you can see the structure of the matrix, with $\sigma=\frac{\alpha\Delta t}{\Delta x^2}$:

$$
\begin{align}
    \left[
        \begin{array}{cccccc}
            \left( 2 + \frac{1}{\sigma} \right) & -1 & 0 & \cdots & & 0 \\
            -1 & \left( 2 + \frac{1}{\sigma} \right) & -1 & 0 & \cdots & 0 \\
            0 & & \ddots & & & \vdots \\
            \vdots & & & & \left( 2 + \frac{1}{\sigma} \right) & \\
            0 & \cdots & & & -1 & \left( 1 + \frac{1}{\sigma} \right)
        \end{array}
    \right] \cdot 
    \left[
        \begin{array}{c} 
            T_1^{n+1} \\
            T_2^{n+1} \\
            \vdots \\ \\
            T_{N-2}^{n+1}
        \end{array}
    \right] =
    \left[
        \begin{array}{c} 
            T_1^n \frac{1}{\sigma} \\
            T_2^{n} \frac{1}{\sigma} \\
            \vdots \\ \\
            T_{N-2}^{n} \frac{1}{\sigma}
        \end{array}
    \right] +
    \begin{bmatrix}
        T_0^{n+1} \\
        0 \\
        \vdots \\
        0 \\
        q\Delta x
    \end{bmatrix}
\end{align}
$$
 
Notice that the Dirichlet boundary condition adds only a term to the right-hand side of the system.  The Neumann boundary condition both adds a term to the right-hand side and modifies the matrix $[A]$.

### Problem set up

We'll re-use the problem from lesson 1: we have a graphite rod, with [thermal diffusivity](http://en.wikipedia.org/wiki/Thermal_diffusivity) $\alpha=1.22\times10^{-3} {\rm m}^2/{\rm s}$, length $L=1{\rm m}$, and temperature held at $T=100{\rm C}$ on the left end, $x=0$, and $0{\rm C}$ everywhere else initially. We'll compute the evolution of temperature on the length of the rod.

Let's start like we did in the previous lesson: import your libraries and set up the discretization.


```python
import numpy
from matplotlib import pyplot
%matplotlib inline
```


```python
# Set the font family and size to use for Matplotlib figures.
pyplot.rcParams['font.family'] = 'serif'
pyplot.rcParams['font.size'] = 16
```


```python
# Set parameters.
L = 1.0  # length of the rod
nx = 51  # number of locations on the rod
dx = L / (nx - 1)  # distance between two consecutive locations
alpha = 1.22e-3  #  thermal diffusivity of the rod
q = 0.0  # temperature gradient on the right side of the rod

# Define the locations along the rod.
x = numpy.linspace(0.0, L, num=nx)

# Set the initial temperature along the rod.
T0 = numpy.zeros(nx)
T0[0] = 100.0
```

### Solving a linear system

We need to solve the linear system of equations written above to advance the solution in time. Luckily, we can rely on our friends from SciPy who have developed some nice linear solvers, so we don't need to write our own.

From `scipy.linalg`, let's import `solve`: a function to solve linear systems. Make sure to explore the documentation of [`scipy.linalg`](https://docs.scipy.org/doc/scipy/reference/linalg.html). We'll need to define our own custom functions to generate the coefficient matrix and the right-hand side of the linear system. You should carefully study the code below.


```python
from scipy import linalg
```


```python
def lhs_operator(N, sigma):
    """
    Computes and returns the implicit operator
    of the system for the 1D diffusion equation.
    We use backward Euler method, Dirichlet condition
    on the left side of the domain and zero-gradient
    Neumann condition on the right side.
    
    Parameters
    ----------
    N : integer
        Number of interior points.
    sigma : float
        Value of alpha * dt / dx**2.
    
    Returns
    -------
    A : numpy.ndarray
        The implicit operator as a 2D array of floats
        of size N by N.
    """
    # Setup the diagonal of the operator.
    D = numpy.diag((2.0 + 1.0 / sigma) * numpy.ones(N))
    # Setup the Neumann condition for the last element.
    D[-1, -1] = 1.0 + 1.0 / sigma
    # Setup the upper diagonal of the operator.
    U = numpy.diag(-1.0 * numpy.ones(N - 1), k=1)
    # Setup the lower diagonal of the operator.
    L = numpy.diag(-1.0 * numpy.ones(N - 1), k=-1)
    # Assemble the operator.
    A = D + U + L
    return A
```


```python
def rhs_vector(T, sigma, qdx):
    """
    Computes and returns the right-hand side of the system
    for the 1D diffusion equation, using a Dirichlet condition
    on the left side and a Neumann condition on the right side.
    
    Parameters
    ----------
    T : numpy.ndarray
        The temperature distribution as a 1D array of floats.
    sigma : float
        Value of alpha * dt / dx**2.
    qdx : float
        Value of the temperature flux at the right side.
    
    Returns
    -------
    b : numpy.ndarray
        The right-hand side of the system as a 1D array of floats.
    """
    b = T[1:-1] / sigma
    # Set Dirichlet condition.
    b[0] += T[0]
    # Set Neumann condition.
    b[-1] += qdx
    return b
```

Next, we'll define a function that steps in time using the implicit central-space scheme.  Remember that for an implicit method, a step in time is performed by solving the entire linear system. This is a fundamental difference between implicit and explicit methods, and implies a considerable computational cost.


```python
def btcs_implicit(T0, nt, dt, dx, alpha, q):
    """
    Computes and returns the temperature along the rod
    after a given number of time steps.
    
    The function uses Euler implicit in time,
    central differencing in space, a Dirichlet condition
    on the left side, and a Neumann condition on the
    right side.
    
    Parameters
    ----------
    T0 : numpy.ndarray
        The initial temperature distribution as a 1D array of floats.
    nt : integer
        Number of time steps to compute.
    dt : float
        Time-step size.
    dx : float
        Distance between two consecutive locations.
    alpha : float
        Thermal diffusivity of the rod.
    q : float
        Value of the temperature gradient on the right side.
    
    Returns
    -------
    T : numpy.ndarray
        The temperature distribution as a 1D array of floats.
    """
    sigma = alpha * dt / dx**2
    # Create the implicit operator of the system.
    A = lhs_operator(len(T0) - 2, sigma)
    # Integrate in time.
    T = T0.copy()
    for n in range(nt):
        # Generate the right-hand side of the system.
        b = rhs_vector(T, sigma, q * dx)
        # Solve the system with scipy.linalg.solve.
        T[1:-1] = linalg.solve(A, b)
        # Apply the Neumann boundary condition.
        T[-1] = T[-2] + q * dx
    return T
```

We solve the linear system for every time step, but the $A$ matrix does not change. Thus, you can generate it only once and then use it for all time steps. Let's try this out!


```python
# Set the time-step size based on CFL limit.
sigma = 0.5
dt = sigma * dx**2 / alpha  # time-step size
nt = 1000  # number of time steps to compute

# Compute the temperature along the rod.
T = btcs_implicit(T0, nt, dt, dx, alpha, q)
```

Now plot the solution!


```python
# Plot the temperature along the rod.
pyplot.figure(figsize=(6.0, 4.0))
pyplot.xlabel('Distance [m]')
pyplot.ylabel('Temperature [C]')
pyplot.grid()
pyplot.plot(x, T, color='C0', linestyle='-', linewidth=2)
pyplot.xlim(0.0, L)
pyplot.ylim(0.0, 100.0);
```

Not too impressive, this looks just like the result from *explicit* forward in time, centered in space for $\alpha\frac{\Delta t}{\Delta x^2} = \frac{1}{2}$. 

But try $\alpha\frac{\Delta t}{\Delta x^2} = 5$, which violates the stability condition of the *explicit* scheme:


```python
# Increase the CFL number.
sigma = 5.0
dt = sigma * dx**2 / alpha  # time-step size
nt = 100  # number of time steps to compute

# Compute the temperature along the rod.
T = btcs_implicit(T0, nt, dt, dx, alpha, q)

# Plot the temperature along the rod.
pyplot.figure(figsize=(6.0, 4.0))
pyplot.xlabel('Distance [m]')
pyplot.ylabel('Temperature [C]')
pyplot.grid()
pyplot.plot(x, T, color='C0', linestyle='-', linewidth=2)
pyplot.xlim(0.0, L)
pyplot.ylim(0.0, 100.0);
```

**It didn't blow up!**

We were not able to use such a large time step with the explicit scheme. You can try out other values of `sigma` and you'll get a stable solution. In fact, this is an *unconditionally stable* scheme—the most valuable feature of implicit methods is that they give stable solutions without a constraint on the choice of time step. 

Using the implicit scheme, we can always advance in time using larger time steps. But each time step requires the solution of a linear system, which is computationally expensive. This is the trade-off between explicit and implicit methods. 
To experiment further, set different values of the Neumann boundary flux and see if the solution behaves as you expect.

##### A word of warning

Implicit methods allow you to use significantly larger time steps, because they are not subject to stability constraints. But that doesn't mean you can use just _any_ large time step! Remember that Euler's method is a first-order method, so the _accuracy_ gets worse as you increase the time step, in direct proportion. In fact, you can lose the ability to capture the correct physics if your time step is too large. Numerical stability does not imply accuracy!

##### Dig deeper

You see how matrix `A` is mostly full of zeros?  We call such a matrix *sparse*, and there are many ways to make more efficient calculations taking advantage of their particular structure. First of all, you can optimize the memory usage. Check out SciPy's [sparse-matrix storage formats](https://docs.scipy.org/doc/scipy/reference/sparse.html): you don't need too store $(N-2)^2$ elements! For example, a `coo_matrix` format stores only $3*N_\text{nonzero}$, where $N_\text{nonzero}$ is the number of non-zero elements in `A`. Make sure to explore this topic a little more. It's an important topic in numerical PDEs.

---
###### The cell below loads the style of the notebook


```python
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, 'r').read())
```




<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>
<link href='https://fonts.googleapis.com/css?family=Source+Code+Pro' rel='stylesheet' type='text/css'>
<style>

@font-face {
    font-family: "Computer Modern";
    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}

#notebook_panel { /* main background */
    background: rgb(245,245,245);
}

div.cell { /* set cell width */
    width: 750px;
}

div #notebook { /* centre the content */
    background: #fff; /* white background for content */
    width: 1000px;
    margin: auto;
    padding-left: 0em;
}

#notebook li { /* More space between bullet points */
    margin-top:0.8em;
}

/* draw border around running cells */
div.cell.border-box-sizing.code_cell.running { 
    border: 1px solid #111;
}

/* Put a solid color box around each cell and its output, visually linking them*/
div.cell.code_cell {
    background-color: rgb(256,256,256); 
    border-radius: 0px; 
    padding: 0.5em;
    margin-left:1em;
    margin-top: 1em;
}

div.text_cell_render{
    font-family: 'Alegreya Sans' sans-serif;
    line-height: 140%;
    font-size: 125%;
    font-weight: 400;
    width:600px;
    margin-left:auto;
    margin-right:auto;
}


/* Formatting for header cells */
.text_cell_render h1 {
    font-family: 'Nixie One', serif;
    font-style:regular;
    font-weight: 400;    
    font-size: 45pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.5em;
    margin-top: 0.5em;
    display: block;
}

.text_cell_render h2 {
    font-family: 'Nixie One', serif;
    font-weight: 400;
    font-size: 30pt;
    line-height: 100%;
    color: rgb(0,51,102);
    margin-bottom: 0.1em;
    margin-top: 0.3em;
    display: block;
}	

.text_cell_render h3 {
    font-family: 'Nixie One', serif;
    margin-top:16px;
    font-size: 22pt;
    font-weight: 600;
    margin-bottom: 3px;
    font-style: regular;
    color: rgb(102,102,0);
}

.text_cell_render h4 {    /*Use this for captions*/
    font-family: 'Nixie One', serif;
    font-size: 14pt;
    text-align: center;
    margin-top: 0em;
    margin-bottom: 2em;
    font-style: regular;
}

.text_cell_render h5 {  /*Use this for small titles*/
    font-family: 'Nixie One', sans-serif;
    font-weight: 400;
    font-size: 16pt;
    color: rgb(163,0,0);
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.8em;
    display: block;
}

.text_cell_render h6 { /*use this for copyright note*/
    font-family: 'PT Mono', sans-serif;
    font-weight: 300;
    font-size: 9pt;
    line-height: 100%;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

.CodeMirror{
    font-family: "Source Code Pro";
    font-size: 90%;
}

.alert-box {
    padding:10px 10px 10px 36px;
    margin:5px;
}

.success {
    color:#666600;
    background:rgb(240,242,229);
}
</style>






2596c36d2365fcadef6c17a623fe6f0006e92d19
ipynb
Jupyter Notebook
lessons/04_spreadout/04_02_Heat_Equation_1D_Implicit.ipynb
mcarpe/numerical-mooc
62b3c14c2c56d85d65c6075f2d7eb44266b49c17
CC-BY-3.0
2015-01-04T22:50:56.000Z
2022-03-30T20:42:16.000Z
lessons/04_spreadout/04_02_Heat_Equation_1D_Implicit.ipynb
mcarpe/numerical-mooc
62b3c14c2c56d85d65c6075f2d7eb44266b49c17
CC-BY-3.0
2015-02-02T01:06:07.000Z
2020-11-09T12:27:41.000Z
lessons/04_spreadout/04_02_Heat_Equation_1D_Implicit.ipynb
mcarpe/numerical-mooc
62b3c14c2c56d85d65c6075f2d7eb44266b49c17
CC-BY-3.0
2015-01-02T19:19:52.000Z
2022-02-27T01:02:44.000Z
Qwen/Qwen-72B
1. YES
2. YES