Partial differential equations: The (parabolic) heat equation

Partial differential equations: The (parabolic) heat equation#

The heat equation is the typical example of a parabolic partial differential equation. Here we have the time and we are looking for an actual time evolution. It can be written as

(56)#\[\begin{equation} \frac{\partial u}{\partial t} = \alpha \nabla^2 u, \end{equation}\]

where \(u\) represents the temperature, or the concentration (diffusion problem) and so on. The coefficient \(\alpha\) takes into account the thermal conductivity, specific heat, and density of the material.

In the following, we will use as main references the book from Chapra and also

  • https://github.com/numerical-mooc/numerical-mooc

  • https://aquaulb.github.io/book_solving_pde_mooc/solving_pde_mooc/notebooks/01_Introduction/01_00_Preface.html

  • Wikipedia

The explicit method: Forward time centered space#

As is usual in finite differences, we discretize both the spatial (up to second order) and the time (linear order) derivatives to get

(57)#\[\begin{equation} \frac{u^{n+1}_i - u^n_i}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}, \end{equation}\]

where the superscript corresponds to a step in time, while the subscript corresponds to step in space (we are here solving a one-dimensional equation). This is the so called FTCS - Forward Time Centered Space method, which is an explicit method which allows to compute the next (in time) value in terms of the previous one.

Reorganizing terms

(58)#\[\begin{equation} u_{i}^{n+1} = u_{i}^{n}+\frac{\alpha\Delta t}{\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}) = u_{i}^{n}+\lambda(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}), \end{equation}\]

which can be visualized as

By using the von Neumann stability analysis it is possible to show that the solution is stable (errors are not amplified) only if

(59)#\[\begin{equation} \lambda = \frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2} \end{equation}\]

or, in 2D with different spacing

(60)#\[\begin{equation} \lambda = \frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2(1 + \frac{\Delta x^2}{\Delta y^2})}. \end{equation}\]

Example: Finding the temperature distribution on a thind rod#

Now let’s solve an example (Chapra 30.1): A thin rod of 10 cm, with \(\alpha = 0.835\) cm\(^2\)/s, \(\Delta t = 0.1 \)s, \(\Delta x = 2\) cm, and \(\lambda = \alpha \Delta t/\Delta x^2 = 0.020875\). The boundary conditions are \(T(0) = 100\) C, and \(T(10) = 50\) C.

# Setup imports
import numpy as np
import matplotlib.pyplot as plt

# Function to evolve
def ftcs(T, deltat, deltax, lambdac):
    Tnew = T.copy() # avoids overwriting T0
    # YOUR CODE HERE
    raise NotImplementedError()
    return Tnew

# Setup problem
L = 10 
DELTAT = 0.10
DELTAX = 2.0
NX = int(L/DELTAX)
ALPHA = 0.835
LAMBDA = ALPHA*DELTAT/(DELTAX*DELTAX)
print(f"{LAMBDA=}")

# Set the initial temperature along the rod.
T0 = np.zeros(NX)
#T0 = np.random.uniform(-100.0, 1000.0, NX)
T0[0] = 100.0
T0[-1] = 50

LAMBDA=0.020875

# Solve the problem several for different niter and plot it
NDISPLAY = 50
fig, ax = plt.subplots()
x = np.linspace(0.0, L, num=NX)
T = T0.copy()
for niter in range(800):
    T = ftcs(T, DELTAT, DELTAX, LAMBDA)
    if niter%NDISPLAY == 0:
        ax.plot(x, T, label=f"{niter=}")
ax.legend()

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[4], line 7
      5 T = T0.copy()
      6 for niter in range(800):
----> 7     T = ftcs(T, DELTAT, DELTAX, LAMBDA)
      8     if niter%NDISPLAY == 0:
      9         ax.plot(x, T, label=f"{niter=}")

Cell In[2], line 5, in ftcs(T, deltat, deltax, lambdac)
      3 Tnew = T.copy() # avoids overwriting T0
      4 # YOUR CODE HERE
----> 5 raise NotImplementedError()
      6 return Tnew

NotImplementedError:
../_images/9e13561b06f5d2c575e2e9d0ab28f3ca3efd7ae4d5738c0cd619a55149090489.png
EXERCISE

  • [ ] Play with \(\lambda\). Is there really any instability?

  • [ ] Play with the initial conditions: is there any dependence?

  • [ ] Create an animation of the temperature time evolution.

A note on derivative boundary conditions#

In case you need to fix the derivative a the boundary, let’s say an isolated wall on the left, then you could introduce that as

(61)#\[\begin{equation} u_0^{n+1} = u_0^n + \lambda(u_1^n - 2u_0^n + u_{-1}^n), \end{equation}\]

where there is a new imaginary point at \(i = -1\). But from the derivative condition you can get its value, since \((u_0 - u_{-1})/\Delta x \simeq \partial u/\partial x\).

Implicit methods#

Implicit methods improve the solution stability and, generally, result in in matrix systems that need to be solve. Examples can be found at https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/04_02_Heat_Equation_1D_Implicit.ipynb . One of the most used implicit method is the Crank-Nicholson method is an implicit method which is commonly used given that it is unconditionally stable. It is second order in time, and its evolution is represented as

For more info check https://github.com/numerical-mooc/numerical-mooc/tree/master/lessons/04_spreadout

Parabolic equations in two dimensions#

In this case you could be solving problems like

(62)#\[\begin{equation} \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}\right). \end{equation}\]

In this case, stability analysis implies that

(63)#\[\begin{equation} \Delta t \le \frac{1}{8} \frac{\Delta x^2 + \Delta y^2}{\alpha}. \end{equation}\]

As with the previous 1D case, implicit schemes improve the stability , but the matrix systems are no longer tri-diagonal and therefore the matrix storage and computational times increase a lot.

There are alternative approaches , like the ADI-Alternating direction implicit scheme which can be use in very complex PDE, even non-linear.

Exercises#

  1. Solve the 1D heat equation with a Newmann condition at the left border , where the derivative is null.

  2. Solve the heat equation in 2D using the ftcs method

Contents