# 15 PDEs: Solution with Time Stepping (Students)

## Heat Equation
The **heat equation** can be derived from Fourier's law and energy conservation (see the [lecture notes on the heat equation (PDF)](15_PDEs_LectureNotes_HeatEquation.pdf))

$$
\frac{\partial T(\mathbf{x}, t)}{\partial t} = \frac{K}{C\rho} \nabla^2 T(\mathbf{x}, t),
$$

## Problem: insulated metal bar (1D heat equation)
A metal bar of length $L$ is insulated along it lengths and held at 0ºC at its ends. Initially, the whole bar is at 100ºC. Calculate $T(x, t)$ for $t>0$.

### Analytic solution
Solve by separation of variables and power series: The general solution that obeys the boundary conditions $T(0, t) = T(L, t) = 0$ is

$$
T(x, t) = \sum_{n=1}^{+\infty} A_n \sin(k_n x)\, \exp\left(-\frac{k_n^2 K t}{C\rho}\right), \quad k_n = \frac{n\pi}{L}
$$

The specific solution that satisfies $T(x, t=0) = T_0 = 100^\circ\text{C}$ leads to $A_n = 4 T_0/n\pi$ for $n$ odd:

$$
T(x, t) = \sum_{n=1,3,5,\dots}^{+\infty} \frac{4 T_0}{n \pi} \sin(k_n x)\, \exp\left(-\frac{k_n^2 K t}{C\rho}\right)
$$


```python
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
```


```python
%matplotlib inline
```

Analytical solution:


```python
def T_bar(x, t, T0, L, K=237, C=900, rho=2700, nmax=1000):
    T = np.zeros_like(x)
    eta = K / (C*rho)
    for n in range(1, nmax, 2):
        kn = n*np.pi/L
        T += 4*T0/(np.pi * n) * np.sin(kn*x) * np.exp(-kn*kn * eta * t)
    return T
```


```python
T0 = 100.
L = 1.0
X = np.linspace(0, L, 100)
for t in np.linspace(0, 3000, 50):
    plt.plot(X, T_bar(X, t, T0, L))
plt.xlabel(r"$x$ (m)")
plt.ylabel(r"$T$ ($^\circ$C)");
```

### Numerical solution: Leap frog
Discretize (finite difference):

For the time domain we only have the initial values so we use a simple forward difference for the time derivative:

$$
\frac{\partial T(x,t)}{\partial t} \approx \frac{T(x, t+\Delta t) - T(x, t)}{\Delta t}
$$

For the spatial derivative we have initially all values so we can use the more accurate central difference approximation:

$$
\frac{\partial^2 T(x, t)}{\partial x^2} \approx \frac{T(x+\Delta x, t) + T(x-\Delta x, t) - 2 T(x, t)}{\Delta x^2}
$$

Thus, the heat equation can be written as the finite difference equation

$$
\frac{T(x, t+\Delta t) - T(x, t)}{\Delta t} = \frac{K}{C\rho} \frac{T(x+\Delta x, t) + T(x-\Delta x, t) - 2 T(x, t)}{\Delta x^2}
$$

which can be reordered so that the RHS contains only known terms and the LHS future terms. Index $i$ is the spatial index, and $j$ the time index: $x = x_0 + i \Delta x$, $t = t_0 + j \Delta t$.

$$
T_{i, j+1} = (1 - 2\eta) T_{i,j} + \eta(T_{i+1,j} + T_{i-1, j}), \quad \eta := \frac{K \Delta t}{C \rho \Delta x^2}
$$

Thus we can step forward in time ("leap frog"), using only known values.

### Activity: Solve the 1D heat equation numerically for an iron bar
* $K = 237$ W/mK
* $C = 900$ J/K
* $\rho = 2700$ kg/m<sup>3</sup>
* $L = 1$ m
* $T_0 = 373$ K and $T_b = 273$ K
* $T(x, 0) = T_0$ and $T(0, t) = T(L, t) = T_b$

Implement the Leapfrog time-stepping algorithm and visualize the results.


```python
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
```

For 3D rotatable output:


```python
%matplotlib widget
```

For HTML/nbviewer output, use inline:


```python
%matplotlib inline
```

Numerical solution


```python
L_rod = 1.    # m
t_max = 3000. # s

Dx = 0.02   # m
Dt = 2    # s

Nx = int(L_rod // Dx)
Nt = int(t_max // Dt)

Kappa = 237 # W/(m K)
CHeat = 900 # J/K
rho = 2700  # kg/m^3

T0 = 373 # K
Tb = 273 # K

raise NotImplementedError
# eta = 


step = 20  # plot solution every n steps

print("Nx = {0}, Nt = {1}".format(Nx, Nt))
print("eta = {0}".format(eta))

T = np.zeros(Nx)
T_new = np.zeros_like(T)
T_plot = np.zeros((Nt//step + 1, Nx))

raise NotImplementedError
# initial conditions
# ...

# boundary conditions
# ...

t_index = 0
T_plot[t_index, :] = T
for jt in range(1, Nt):
    
    raise NotImplementedError
    
    if jt % step == 0 or jt == Nt-1:
        t_index += 1
        # save the new solution for later plotting
        # T_plot[t_index, :] = 
        print("Iteration {0:5d}".format(jt), end="\r")
else:
    print("Completed {0:5d} iterations: t={1} s".format(jt, jt*Dt))
```

#### Visualization
Visualize (you can use the code as is). 

Note how we are making the plot use proper units by mutiplying with `Dt * step` and `Dx`.


```python
X, Y = np.meshgrid(range(T_plot.shape[0]), range(T_plot.shape[1]))
Z = T_plot[X, Y]
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_wireframe(X*Dt*step, Y*Dx, Z)
ax.set_xlabel(r"time $t$ (s)")
ax.set_ylabel(r"position $x$ (m)")
ax.set_zlabel(r"temperature $T$ (K)")
fig.tight_layout()
```

## Stability of the solution

### Empirical investigation of the stability
Investigate the solution for different values of `Dt` and `Dx`. Can you discern patters for stable/unstable solutions?

Report `Dt`, `Dx`, and `eta`
* for 3 stable solutions 
* for 3 unstable solutions


Wrap your heat diffusion solver in a function so that it becomes easier to run it:


```python
def calculate_T(L_rod=1, t_max=3000, Dx=0.02, Dt=2, T0=373, Tb=273,
               step=20):
    Nx = int(L_rod // Dx)
    Nt = int(t_max // Dt)

    Kappa = 237 # W/(m K)
    CHeat = 900 # J/K
    rho = 2700  # kg/m^3

    raise NotImplementedError
    
    return T_plot

def plot_T(T_plot, Dx, Dt, step):
    X, Y = np.meshgrid(range(T_plot.shape[0]), range(T_plot.shape[1]))
    Z = T_plot[X, Y]
    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    ax.plot_wireframe(X*Dt*step, Y*Dx, Z)
    ax.set_xlabel(r"time $t$ (s)")
    ax.set_ylabel(r"position $x$ (m)")
    ax.set_zlabel(r"temperature $T$ (K)")
    fig.tight_layout()
    return ax
```


```python
T_plot = calculate_T(Dx=0.02, Dt=2, step=20)
plot_T(T_plot, 0.02, 2, 20)
```

For which values of $\Delta t$ and $\Delta x$ does the solution become unstable?

### Von Neumann stability analysis 

If the difference equation solution diverges then we *know* that we have a bad approximation to the original PDE. 

Von Neumann stability analysis starts from the assumption that *eigenmodes* of the difference equation can be written as

$$
T_{m,j} = \xi(k)^j e^{ikm\Delta x}, \quad t=j\Delta t,\ x=m\Delta x 
$$

with the unknown wave vectors $k=2\pi/\lambda$ and unknown complex functions $\xi(k)$.

Solutions of the difference equation can be written as linear superpositions of these basis functions. But they are only stable if the eigenmodes are stable, i.e., will not grow in time (with $j$). This is the case when 
$$
|\xi(k)| < 1
$$
for all $k$.

Insert the eigenmodes into the finite difference equation

$$
T_{m, j+1} = (1 - 2\eta) T_{m,j} + \eta(T_{m+1,j} + T_{m-1, j})
$$

to obtain 

\begin{align}
\xi(k)^{j+1} e^{ikm\Delta x} &= (1 - 2\eta) \xi(k)^{j} e^{ikm\Delta x} 
    + \eta(\xi(k)^{j} e^{ik(m+1)\Delta x} + \xi(k)^{j} e^{ik(m-1)\Delta x})\\
\xi(k) &= (1 - 2\eta) + \eta(e^{ik\Delta x} + e^{-ik\Delta x})\\
\xi(k) &= 1 - 2\eta + 2\eta \cos k\Delta x\\
\xi(k) &= 1 + 2\eta\big(\cos k\Delta x - 1\big)
\end{align}

For $|\xi(k)| < 1$ (and all possible $k$):

\begin{align}
|\xi(k)| < 1 \quad &\Leftrightarrow \quad \xi^2(k) < 1\\
(1 + 2y)^2 = 1 + 4y + 4y^2 &< 1 \quad \text{with}\ \  y = \eta(\cos k\Delta x - 1)\\
y(1 + y) &< 0 \quad \Leftrightarrow \quad -1 < y < 0\\
\eta(\cos k\Delta x - 1) &\leq 0 \quad \forall k \quad (\eta > 0, -1 \leq \cos x \leq 1)\\
\eta(\cos k\Delta x - 1) &> -1\\
\eta &< \frac{1}{1 - \cos k\Delta x}\\
\eta &= \frac{K \Delta t}{C \rho \Delta x^2} < \frac{1}{2}
\end{align}

Thus, solutions are only stable for $\eta < 1/2$. In particular, decreasing $\Delta t$ will always improve stability, But decreasing $\Delta x$ requires an quadratic *increase* in $\Delta t$!

Note
* Perform von Neumann stability analysis when possible (depends on PDE and the specific discretization).
* Test different combinations of $\Delta t$ and $\Delta x$.
* Not guarantee that decreasing both will lead to more stable solutions!


```python

```

09d631bb150c8681122b52e3574f5d599bc08336
ipynb
Jupyter Notebook
15_PDEs/15_PDEs-Students.ipynb
Py4Phy/PHY432-resources
c26d95eaf5c28e25da682a61190e12ad6758a938
CC-BY-4.0
15_PDEs/15_PDEs-Students.ipynb
Py4Phy/PHY432-resources
c26d95eaf5c28e25da682a61190e12ad6758a938
CC-BY-4.0
2022-03-03T21:47:56.000Z
2022-03-03T21:47:56.000Z
15_PDEs/15_PDEs-Students.ipynb
Py4Phy/PHY432-resources
c26d95eaf5c28e25da682a61190e12ad6758a938
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES