Partial Differential Equations - Executable Engineering: Interactive Numerical Methods in Python

Partial differential equations (PDEs) expand ordinary differential equations (ODEs) to include multiple independent variables (coordinates):

Additional spatial coordinates are typically extensions of the ODE case and are usually posed and treated as boundary value problems (BVPs).
The time coordinate is almost always posed as an initial value problem (IVP) and tends to introduce complexity in the form of numerical stability constraints.

Classification of PDEs¶

Introducing time as a coordinate in PDEs raises the fundamental question of how information propagates from the initial condition or from events that occur during the solution.

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Elliptic Equations¶

Second-order equations of the form (with $B^2 - 4AC < 0$ ):

A\frac{\partial^2 u}{\partial x^2} + B\frac{\partial^2 u}{\partial x \partial y} + C\frac{\partial^2 u}{\partial y^2} = f(x, y)

(1)

are termed elliptic. The absence of a time-derivative term implies that information is not being propagated over time, but rather is instantaneously available everywhere in space—all transients have decayed and the system has reached a steady state (even if coefficients such as $A(t)$ vary quasi-statically).

Laplace’s and Poisson’s equations ( $\nabla^2 u = 0$ , $\nabla^2 u = f(\mathbf{x})$ ): electrostatic potential with charge density, steady-state heat conduction with a source, or velocity potential of an incompressible, irrotational fluid.
Helmholtz equation ( $\nabla^2 u + k^2 u = 0$ ): acoustic vibrations of a membrane or cavity, and quantum-mechanical wave functions in a potential well (obtained after separating time and space variables).

Solution Approaches¶

Solving elliptic equations requires solving a coupled system of algebraic (possibly nonlinear) equations simultaneously across the entire spatial mesh:

Initial guess: Nonlinear root finders require a sufficiently close initial guess. A common mitigation strategy is continuation (parametrically ramping up a physical parameter, such as gradually increasing the heat load in natural convection).
Convergence: Discretizing a fine mesh yields a high-dimensional system where root finders may wander into local minima or leave the feasible region; trust-region methods help mitigate this.
Parallelization: Under domain decomposition, elliptic coupling implies every subdomain instantaneously influences all others, creating heavy communication overhead unless accelerated by multigrid methods and preconditioners.

Parabolic Equations¶

Equations of the form:

A \frac{\partial u}{\partial t} - C \frac{\partial^2 u}{\partial x^2} = 0

(2)

propagate information in one temporal direction (forward in time) and govern diffusive processes. Diffusion naturally smooths out sharp discontinuities and high-frequency spatial oscillations, though numerical time-stepping can still suffer from severe stability restrictions.

Solution Schemes¶

Diffusion problems are solved via a time-marching scheme in which time is discretized as an IVP while the spatial operator is discretized as a BVP at each step:

Stiffness: Fine spatial grids introduce widely separated decay time scales, making parabolic PDEs stiff.
Ill-posed inverse problems: When back-calculating model parameters or initial states from later measurements (e.g., inferring thermal conductivity from a cooling profile), the smoothing nature of forward diffusion works against us—many distinct initial states diffuse into nearly identical smooth profiles.

Hyperbolic Equations¶

Equations of the form:

A \frac{\partial^2 u}{\partial t^2} - C \frac{\partial^2 u}{\partial x^2} = 0

(3)

propagate information both forward and backward along characteristic lines at a finite wave speed. Replacing $t$ with $-t$ leaves the governing equation invariant (time-reversal symmetry). Hyperbolic PDEs describe wave propagation (second-order) as well as pure advection (first-order, e.g., $A \frac{\partial u}{\partial t} + C \frac{\partial u}{\partial x} = 0$ ).

Solution Schemes¶

Reduction of order: Second-order time derivatives can be decomposed into a coupled system of first-order equations in time.
Spectral methods: Systems with periodic boundary conditions are well-suited to Fourier spectral analysis, which avoids artificial numerical dispersion.
Reversibility: Time-reversal symmetry makes hyperbolic systems amenable to inverse reconstruction of earlier states from later wave measurements.

Solving Time-Dependent PDEs¶

The standard approach to solving time-dependent PDEs combines IVP time-marching with spatial BVP discretization:

Discretize the time derivative as an IVP relating the new state $T^{t+1}$ to the known state $T^t$ .
Discretize the spatial derivatives on a mesh of spacing $\Delta x$ using finite differences.
March forward in time, enforcing boundary conditions at each time step.

Consider the 1D parabolic heat equation with thermal diffusivity $\alpha > 0$ :

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \qquad \alpha \frac{\partial^2 T}{\partial x^2} \approx \alpha \frac{T_{i-1} - 2 T_i + T_{i+1}}{\Delta x^2}

(4)

Notation¶

Spatial grid points are indexed with a subscript $i$ ( $x_i = i \Delta x$ ).
Discrete time levels are indexed with a superscript $t$ , marching from known $T_i^t$ to unknown $T_i^{t+1}$ .

Explicit Euler Time-Stepping¶

Evaluating the spatial second derivative at the current time level $t$ yields the explicit (Forward Euler) scheme:

\begin{aligned} \frac{T^{t+1}_i - T^{t}_i}{\Delta t} &= \alpha \frac{T^t_{i-1} - 2 T^t_i + T^t_{i+1}}{\Delta x^2} \\ T^{t+1}_i &= T_i^t + r \left(T^t_{i-1} - 2 T^t_i + T^t_{i+1}\right), \qquad r = \frac{\alpha \Delta t}{\Delta x^2} \end{aligned}

(5)

![forward Euler stencil.svg](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0iVVRGLTgiIHN0YW5kYWxvbmU9Im5vIj8+CjwhLS0gQ3JlYXRlZCB3aXRoIElua3NjYXBlIChodHRwOi8vd3d3Lmlua3NjYXBlLm9yZy8pIC0tPgoKPHN2ZwogICB3aWR0aD0iNDUuMjA2OTgybW0iCiAgIGhlaWdodD0iMTkuNDI3MTE2bW0iCiAgIHZpZXdCb3g9IjAgMCA0NS4yMDY5ODIgMTkuNDI3MTE2IgogICB2ZXJzaW9uPSIxLjEiCiAgIGlkPSJzdmcxIgogICBpbmtzY2FwZTpleHBvcnQtZmlsZW5hbWU9ImZvcndhcmQgRXVsZXIgc3RlbmNpbC5zdmciCiAgIGlua3NjYXBlOmV4cG9ydC14ZHBpPSI5NiIKICAgaW5rc2NhcGU6ZXhwb3J0LXlkcGk9Ijk2IgogICB4bWxuczppbmtzY2FwZT0iaHR0cDovL3d3dy5pbmtzY2FwZS5vcmcvbmFtZXNwYWNlcy9pbmtzY2FwZSIKICAgeG1sbnM6c29kaXBvZGk9Imh0dHA6Ly9zb2RpcG9kaS5zb3VyY2Vmb3JnZS5uZXQvRFREL3NvZGlwb2RpLTAuZHRkIgogICB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciCiAgIHhtbG5zOnN2Zz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciPgogIDxzb2RpcG9kaTpuYW1lZHZpZXcKICAgICBpZD0ibmFtZWR2aWV3MSIKICAgICBwYWdlY29sb3I9IiNmZmZmZmYiCiAgICAgYm9yZGVyY29sb3I9IiMwMDAwMDAiCiAgICAgYm9yZGVyb3BhY2l0eT0iMC4yNSIKICAgICBpbmtzY2FwZTpzaG93cGFnZXNoYWRvdz0iMiIKICAgICBpbmtzY2FwZTpwYWdlb3BhY2l0eT0iMC4wIgogICAgIGlua3NjYXBlOnBhZ2VjaGVja2VyYm9hcmQ9IjAiCiAgICAgaW5rc2NhcGU6ZGVza2NvbG9yPSIjZDFkMWQxIgogICAgIGlua3NjYXBlOmRvY3VtZW50LXVuaXRzPSJtbSIgLz4KICA8ZGVmcwogICAgIGlkPSJkZWZzMSIgLz4KICA8ZwogICAgIGlua3NjYXBlOmxhYmVsPSJMYXllciAxIgogICAgIGlua3NjYXBlOmdyb3VwbW9kZT0ibGF5ZXIiCiAgICAgaWQ9ImxheWVyMSIKICAgICB0cmFuc2Zvcm09InRyYW5zbGF0ZSgtNDIuMzc2MzgyLC02OC45Nzg3ODEpIj4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9InBhdGgxIgogICAgICAgY3g9IjQ0LjIwMjk1IgogICAgICAgY3k9Ijg2LjU3OTMzIgogICAgICAgcj0iMS44MjY1NjgyIgogICAgICAgaW5rc2NhcGU6ZXhwb3J0LWZpbGVuYW1lPSIuXHBhdGgxLnN2ZyIKICAgICAgIGlua3NjYXBlOmV4cG9ydC14ZHBpPSI5NiIKICAgICAgIGlua3NjYXBlOmV4cG9ydC15ZHBpPSI5NiIgLz4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTEiCiAgICAgICBjeD0iNTguMjI1ODgiCiAgICAgICBjeT0iODYuNTc5MzMiCiAgICAgICByPSIxLjgyNjU2ODIiIC8+CiAgICA8Y2lyY2xlCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgIGlkPSJjaXJjbGUyIgogICAgICAgY3g9IjcyLjI0ODgxIgogICAgICAgY3k9Ijg2LjU3OTMzIgogICAgICAgcj0iMS44MjY1NjgyIiAvPgogICAgPGNpcmNsZQogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBpZD0iY2lyY2xlMyIKICAgICAgIGN4PSI1OC4yMjU4OCIKICAgICAgIGN5PSI3MS4yMzM0ODIiCiAgICAgICByPSIxLjgyNjU2ODIiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBkPSJNIDU4LjA4NDg3LDcxLjYwMTQ3NSBWIDg2LjIxNDAyIgogICAgICAgaWQ9InBhdGgzIiAvPgogICAgPHBhdGgKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgZD0iTSA0NC4yMDI5NTMsODYuNTc5MzM2IEggNzIuNjk3NDE2IgogICAgICAgaWQ9InBhdGg0IiAvPgogICAgPGNpcmNsZQogICAgICAgaWQ9InBhdGg1IgogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6bm9uZTtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBjeD0iNzIuNjk3NDE4IgogICAgICAgY3k9Ijg2LjU3OTMzOCIKICAgICAgIHI9IjAuNzkzNzQ5OTkiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODM7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgIGQ9Im0gNDQuMjAyOTUzLDg2LjU3OTMzNiBoIDI4LjEyOTE1IgogICAgICAgaWQ9InBhdGg2IiAvPgogICAgPGNpcmNsZQogICAgICAgaWQ9InBhdGg3IgogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6bm9uZTtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBjeD0iNzIuMzMyMTA4IgogICAgICAgY3k9Ijg2LjIxNDAyIgogICAgICAgcj0iMC4zOTY4NzQ5OSIgLz4KICAgIDxwYXRoCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MztzdHJva2Utb3BhY2l0eToxIgogICAgICAgZD0iTSA1OC4yMjU4OCw4Ni40MDY2MjEgViA3MS4yOTUxNTEiCiAgICAgICBpZD0icGF0aDgiIC8+CiAgICA8dGV4dAogICAgICAgeG1sOnNwYWNlPSJwcmVzZXJ2ZSIKICAgICAgIHN0eWxlPSJmb250LXNpemU6NC4yMzMzM3B4O2ZpbGw6IzAwMDAwMDtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODM7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgIHk9Ijg3LjgyNzM2MiIKICAgICAgIGlkPSJ0ZXh0OCI+PHRzcGFuCiAgICAgICAgIHNvZGlwb2RpOnJvbGU9ImxpbmUiCiAgICAgICAgIGlkPSJ0c3BhbjgiCiAgICAgICAgIHN0eWxlPSJzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgICAgeT0iODcuODI3MzYyIj50PC90c3Bhbj48L3RleHQ+CiAgICA8dGV4dAogICAgICAgeG1sOnNwYWNlPSJwcmVzZXJ2ZSIKICAgICAgIHN0eWxlPSJmb250LXNpemU6NC4yMzMzM3B4O2ZpbGw6IzAwMDAwMDtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODM7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgIHk9IjcyLjE5OTI1NyIKICAgICAgIGlkPSJ0ZXh0OSI+PHRzcGFuCiAgICAgICAgIHNvZGlwb2RpOnJvbGU9ImxpbmUiCiAgICAgICAgIGlkPSJ0c3BhbjkiCiAgICAgICAgIHN0eWxlPSJzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgICAgeT0iNzIuMTk5MjU3Ij50KzE8L3RzcGFuPjwvdGV4dD4KICA8L2c+Cjwvc3ZnPgo=)

Each grid point $T_i^{t+1}$ at the new time step depends only on three known values at time level $t$ . Use the slider below to increase $\Delta t$ past the stability threshold $r = 0.5$ .

differential_equations_module.py

def forward_euler_heat(T0, alpha, dx, dt, num_steps):
    T = T0.copy()
    T_history = [T.copy()]
    r = alpha * dt / dx**2
    for _ in range(num_steps):
        T[1:-1] = T[1:-1] + r * (T[:-2] - 2.0 * T[1:-1] + T[2:])
        T_history.append(T.copy())
    return np.array(T_history)

exe.display_source(exe.forward_euler_heat)
fig = exe.plot_forward_euler_heat()
fig.show()

Von Neumann Stability Analysis¶

To determine when spatial oscillations grow or decay, von Neumann stability analysis expands the spatial profile in Fourier modes ${}^k T_i^t = e^{j k i \Delta x}$ (with wave number $k$ and $j = \sqrt{-1}$ ) and defines the amplification factor $G$ via ${}^k T_i^{t+1} = G \, {}^k T_i^t$ . Stability requires $|G| \le 1$ for all wave numbers $k$ .

Substituting a Fourier mode into the Forward Euler update gives:

\begin{aligned} {}^k T^{t+1}_i &= e^{j k i \Delta x} + r \left(e^{j k (i-1) \Delta x} - 2 e^{j k i \Delta x} + e^{j k (i+1) \Delta x}\right) \\ &= {}^k T^{t}_i \left[1 + r \left(e^{j k \Delta x} + e^{-j k \Delta x} - 2\right)\right] \\ \implies G &= 1 + 2 r \left[\cos(k \Delta x) - 1\right] \end{aligned}

(6)

Since $\cos(k \Delta x) \in [-1, 1]$ , the amplification factor lies in $[1 - 4r, 1]$ . Requiring $1 - 4r \ge -1$ yields the 1D stability limit (and its 2D/3D extensions on a uniform grid of spacing $h = \Delta x$ ):

\text{1D: } \Delta t \le \frac{\Delta x^2}{2 \alpha}, \qquad \text{2D: } \Delta t \le \frac{h^2}{4 \alpha}, \qquad \text{3D: } \Delta t \le \frac{h^2}{6 \alpha}

(7)

Because $\Delta t$ scales with $\Delta x^2$ , halving the grid spacing forces a fourfold reduction in time step—a severe stiffness penalty shared by explicit diffusion schemes.

Implicit Euler Time-Stepping¶

Evaluating the spatial finite difference at the new time level $t+1$ yields the Backward (Implicit) Euler scheme:

-r T^{t+1}_{i-1} + (1 + 2r) T^{t+1}_i - r T^{t+1}_{i+1} = T_i^{t}

(8)

![backward Euler stencil.svg](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0iVVRGLTgiIHN0YW5kYWxvbmU9Im5vIj8+CjwhLS0gQ3JlYXRlZCB3aXRoIElua3NjYXBlIChodHRwOi8vd3d3Lmlua3NjYXBlLm9yZy8pIC0tPgoKPHN2ZwogICB3aWR0aD0iNDUuMjA2OTgybW0iCiAgIGhlaWdodD0iMTkuNDI3MTQxbW0iCiAgIHZpZXdCb3g9IjAgMCA0NS4yMDY5ODIgMTkuNDI3MTQxIgogICB2ZXJzaW9uPSIxLjEiCiAgIGlkPSJzdmcxIgogICBpbmtzY2FwZTpleHBvcnQtZmlsZW5hbWU9ImZvcndhcmQgRXVsZXIgc3RlbmNpbC5zdmciCiAgIGlua3NjYXBlOmV4cG9ydC14ZHBpPSI5NiIKICAgaW5rc2NhcGU6ZXhwb3J0LXlkcGk9Ijk2IgogICB4bWxuczppbmtzY2FwZT0iaHR0cDovL3d3dy5pbmtzY2FwZS5vcmcvbmFtZXNwYWNlcy9pbmtzY2FwZSIKICAgeG1sbnM6c29kaXBvZGk9Imh0dHA6Ly9zb2RpcG9kaS5zb3VyY2Vmb3JnZS5uZXQvRFREL3NvZGlwb2RpLTAuZHRkIgogICB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciCiAgIHhtbG5zOnN2Zz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciPgogIDxzb2RpcG9kaTpuYW1lZHZpZXcKICAgICBpZD0ibmFtZWR2aWV3MSIKICAgICBwYWdlY29sb3I9IiNmZmZmZmYiCiAgICAgYm9yZGVyY29sb3I9IiMwMDAwMDAiCiAgICAgYm9yZGVyb3BhY2l0eT0iMC4yNSIKICAgICBpbmtzY2FwZTpzaG93cGFnZXNoYWRvdz0iMiIKICAgICBpbmtzY2FwZTpwYWdlb3BhY2l0eT0iMC4wIgogICAgIGlua3NjYXBlOnBhZ2VjaGVja2VyYm9hcmQ9IjAiCiAgICAgaW5rc2NhcGU6ZGVza2NvbG9yPSIjZDFkMWQxIgogICAgIGlua3NjYXBlOmRvY3VtZW50LXVuaXRzPSJtbSIgLz4KICA8ZGVmcwogICAgIGlkPSJkZWZzMSIgLz4KICA8ZwogICAgIGlua3NjYXBlOmxhYmVsPSJMYXllciAxIgogICAgIGlua3NjYXBlOmdyb3VwbW9kZT0ibGF5ZXIiCiAgICAgaWQ9ImxheWVyMSIKICAgICB0cmFuc2Zvcm09InRyYW5zbGF0ZSgtNDIuMzc2MzgyLC0xMDEuNzg3MTEpIj4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTkiCiAgICAgICBjeD0iNDQuMjAyOTUiCiAgICAgICBjeT0iLTEwNC4wNDE4NSIKICAgICAgIHI9IjEuODI2NTY4MiIKICAgICAgIGlua3NjYXBlOmV4cG9ydC1maWxlbmFtZT0iYmFja3dhcmQgRXVsZXIgc3RlbmNpbC5zdmciCiAgICAgICBpbmtzY2FwZTpleHBvcnQteGRwaT0iOTYiCiAgICAgICBpbmtzY2FwZTpleHBvcnQteWRwaT0iOTYiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPGNpcmNsZQogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBpZD0iY2lyY2xlMTAiCiAgICAgICBjeD0iNTguMjI1ODgiCiAgICAgICBjeT0iLTEwNC4wNDE4NSIKICAgICAgIHI9IjEuODI2NTY4MiIKICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMSwtMSkiIC8+CiAgICA8Y2lyY2xlCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgIGlkPSJjaXJjbGUxMSIKICAgICAgIGN4PSI3Mi4yNDg4MSIKICAgICAgIGN5PSItMTA0LjA0MTg1IgogICAgICAgcj0iMS44MjY1NjgyIgogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxLC0xKSIgLz4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTEyIgogICAgICAgY3g9IjU4LjIyNTg4IgogICAgICAgY3k9Ii0xMTkuMzg3NjgiCiAgICAgICByPSIxLjgyNjU2ODIiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPHBhdGgKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgZD0iTSA1OC4wODQ4NywxMTkuMDE5NjggViAxMDQuNDA3MDYiCiAgICAgICBpZD0icGF0aDEyIiAvPgogICAgPHBhdGgKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgZD0iTSA0NC4yMDI5NTMsMTA0LjA0MTc0IEggNzIuNjk3NDE2IgogICAgICAgaWQ9InBhdGgxMyIgLz4KICAgIDxjaXJjbGUKICAgICAgIGlkPSJjaXJjbGUxMyIKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlOm5vbmU7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgY3g9IjcyLjY5NzQxOCIKICAgICAgIGN5PSItMTA0LjA0MTg1IgogICAgICAgcj0iMC43OTM3NDk5OSIKICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMSwtMSkiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODM7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgIGQ9Im0gNDQuMjAyOTUzLDEwNC4wNDE3NCBoIDI4LjEyOTE1IgogICAgICAgaWQ9InBhdGgxNCIgLz4KICAgIDxjaXJjbGUKICAgICAgIGlkPSJjaXJjbGUxNCIKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlOm5vbmU7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgY3g9IjcyLjMzMjEwOCIKICAgICAgIGN5PSItMTA0LjQwNzE1IgogICAgICAgcj0iMC4zOTY4NzQ5OSIKICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMSwtMSkiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODM7c3Ryb2tlLW9wYWNpdHk6MSIKICAgICAgIGQ9Im0gNTguMjI1ODgsMTA0LjIxNDQ2IHYgMTUuMTExNTUiCiAgICAgICBpZD0icGF0aDE1IiAvPgogICAgPHRleHQKICAgICAgIHhtbDpzcGFjZT0icHJlc2VydmUiCiAgICAgICBzdHlsZT0iZm9udC1zaXplOjQuMjMzMzNweDtmaWxsOiMwMDAwMDA7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICB4PSI4MC4wNzc4ODgiCiAgICAgICB5PSIxMjAuNjM1NyIKICAgICAgIGlkPSJ0ZXh0MTUiPjx0c3BhbgogICAgICAgICBzb2RpcG9kaTpyb2xlPSJsaW5lIgogICAgICAgICBpZD0idHNwYW4xNSIKICAgICAgICAgc3R5bGU9InN0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgICAgeD0iODAuMDc3ODg4IgogICAgICAgICB5PSIxMjAuNjM1NyI+dDwvdHNwYW4+PC90ZXh0PgogICAgPHRleHQKICAgICAgIHhtbDpzcGFjZT0icHJlc2VydmUiCiAgICAgICBzdHlsZT0iZm9udC1zaXplOjQuMjMzMzNweDtmaWxsOiMwMDAwMDA7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICB4PSI4MC4wNzc4ODgiCiAgICAgICB5PSIxMDUuMDA3NTgiCiAgICAgICBpZD0idGV4dDE2Ij48dHNwYW4KICAgICAgICAgc29kaXBvZGk6cm9sZT0ibGluZSIKICAgICAgICAgaWQ9InRzcGFuMTYiCiAgICAgICAgIHN0eWxlPSJzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgICAgeT0iMTA1LjAwNzU4Ij50KzE8L3RzcGFuPjwvdGV4dD4KICA8L2c+Cjwvc3ZnPgo=)

Von Neumann analysis gives the amplification factor $G = \frac{1}{1 + 2r[1 - \cos(k \Delta x)]} \in \left[\frac{1}{1 + 4r}, 1\right]$ , which satisfies $|G| \le 1$ for all $\Delta t > 0$ . Backward Euler is unconditionally stable, at the cost of solving a tridiagonal linear system $\mathbf{A} \mathbf{T}^{t+1} = \mathbf{T}^t$ at each time step.

differential_equations_module.py

def backward_euler_heat(T0, alpha, dx, dt, num_steps):
    T = T0.copy()
    T_history = [T.copy()]
    n = len(T)
    r = alpha * dt / dx**2
    A = (
        np.diag((1.0 + 2.0 * r) * np.ones(n), 0)
        + np.diag(-r * np.ones(n - 1), -1)
        + np.diag(-r * np.ones(n - 1), 1)
    )
    A[0, :] = 0.0
    A[0, 0] = 1.0
    A[-1, :] = 0.0
    A[-1, -1] = 1.0
    for _ in range(num_steps):
        b = T.copy()
        T = np.linalg.solve(A, b)
        T_history.append(T.copy())
    return np.array(T_history)

exe.display_source(exe.backward_euler_heat)
fig = exe.plot_backward_euler_heat()
fig.show()

Crank-Nicolson Method¶

Backward Euler is only first-order accurate in time ( $\mathcal{O}(\Delta t)$ ), whereas central differencing in space is second-order ( $\mathcal{O}(\Delta x^2)$ ). The Crank-Nicolson method achieves second-order accuracy in both time and space by averaging the explicit and implicit spatial discretizations at levels $t$ and $t+1$ :

\begin{aligned} \frac{T^{t+1}_i - T^{t}_i}{\Delta t} &= \frac{\alpha}{2 \Delta x^2} \left[\left(T^{t+1}_{i-1} - 2 T^{t+1}_i + T^{t+1}_{i+1}\right) + \left(T^{t}_{i-1} - 2 T^{t}_i + T^{t}_{i+1}\right)\right] \\ -r T_{i-1}^{t+1} + (1 + 2r) T_i^{t+1} - r T_{i+1}^{t+1} &= r T_{i-1}^{t} + (1 - 2r) T_i^{t} + r T_{i+1}^{t}, \qquad r = \frac{\alpha \Delta t}{2 \Delta x^2} \end{aligned}

(9)

![Crank Nicholson stencil.svg](data:image/svg+xml;base64,PD94bWwgdmVyc2lvbj0iMS4wIiBlbmNvZGluZz0iVVRGLTgiIHN0YW5kYWxvbmU9Im5vIj8+CjwhLS0gQ3JlYXRlZCB3aXRoIElua3NjYXBlIChodHRwOi8vd3d3Lmlua3NjYXBlLm9yZy8pIC0tPgoKPHN2ZwogICB3aWR0aD0iNDUuMjA2OTgybW0iCiAgIGhlaWdodD0iMTkuMDI2MTczbW0iCiAgIHZpZXdCb3g9IjAgMCA0NS4yMDY5ODIgMTkuMDI2MTczIgogICB2ZXJzaW9uPSIxLjEiCiAgIGlkPSJzdmcxIgogICBpbmtzY2FwZTpleHBvcnQtZmlsZW5hbWU9ImZvcndhcmQgRXVsZXIgc3RlbmNpbC5zdmciCiAgIGlua3NjYXBlOmV4cG9ydC14ZHBpPSI5NiIKICAgaW5rc2NhcGU6ZXhwb3J0LXlkcGk9Ijk2IgogICB4bWxuczppbmtzY2FwZT0iaHR0cDovL3d3dy5pbmtzY2FwZS5vcmcvbmFtZXNwYWNlcy9pbmtzY2FwZSIKICAgeG1sbnM6c29kaXBvZGk9Imh0dHA6Ly9zb2RpcG9kaS5zb3VyY2Vmb3JnZS5uZXQvRFREL3NvZGlwb2RpLTAuZHRkIgogICB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciCiAgIHhtbG5zOnN2Zz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciPgogIDxzb2RpcG9kaTpuYW1lZHZpZXcKICAgICBpZD0ibmFtZWR2aWV3MSIKICAgICBwYWdlY29sb3I9IiNmZmZmZmYiCiAgICAgYm9yZGVyY29sb3I9IiMwMDAwMDAiCiAgICAgYm9yZGVyb3BhY2l0eT0iMC4yNSIKICAgICBpbmtzY2FwZTpzaG93cGFnZXNoYWRvdz0iMiIKICAgICBpbmtzY2FwZTpwYWdlb3BhY2l0eT0iMC4wIgogICAgIGlua3NjYXBlOnBhZ2VjaGVja2VyYm9hcmQ9IjAiCiAgICAgaW5rc2NhcGU6ZGVza2NvbG9yPSIjZDFkMWQxIgogICAgIGlua3NjYXBlOmRvY3VtZW50LXVuaXRzPSJtbSIgLz4KICA8ZGVmcwogICAgIGlkPSJkZWZzMSIgLz4KICA8ZwogICAgIGlua3NjYXBlOmxhYmVsPSJMYXllciAxIgogICAgIGlua3NjYXBlOmdyb3VwbW9kZT0ibGF5ZXIiCiAgICAgaWQ9ImxheWVyMSIKICAgICB0cmFuc2Zvcm09InRyYW5zbGF0ZSgtNDIuMzc2MzgyLC0xMjYuMTI4NzUpIj4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTE2IgogICAgICAgY3g9IjQ0LjIwMjk1IgogICAgICAgY3k9Ii0xMjguMzgzNTEiCiAgICAgICByPSIxLjgyNjU2ODIiCiAgICAgICBpbmtzY2FwZTpleHBvcnQtZmlsZW5hbWU9ImJhY2t3YXJkIEV1bGVyIHN0ZW5jaWwuc3ZnIgogICAgICAgaW5rc2NhcGU6ZXhwb3J0LXhkcGk9Ijk2IgogICAgICAgaW5rc2NhcGU6ZXhwb3J0LXlkcGk9Ijk2IgogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxLC0xKSIgLz4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTE3IgogICAgICAgY3g9IjU4LjIyNTg4IgogICAgICAgY3k9Ii0xMjguMzgzNTEiCiAgICAgICByPSIxLjgyNjU2ODIiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPGNpcmNsZQogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBpZD0iY2lyY2xlMTgiCiAgICAgICBjeD0iNzIuMjQ4ODEiCiAgICAgICBjeT0iLTEyOC4zODM1MSIKICAgICAgIHI9IjEuODI2NTY4MiIKICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMSwtMSkiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBkPSJNIDU4LjA4NDg3LDE0My4zNjE1IFYgMTI4Ljc0ODg4IgogICAgICAgaWQ9InBhdGgxOSIgLz4KICAgIDxwYXRoCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgIGQ9Ik0gNDQuMjAyOTUzLDEyOC4zODM1NiBIIDcyLjY5NzQxNiIKICAgICAgIGlkPSJwYXRoMjAiIC8+CiAgICA8Y2lyY2xlCiAgICAgICBpZD0iY2lyY2xlMjAiCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZTpub25lO3N0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgIGN4PSI3Mi42OTc0MTgiCiAgICAgICBjeT0iLTEyOC4zODM1MSIKICAgICAgIHI9IjAuNzkzNzQ5OTkiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPHBhdGgKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICBkPSJtIDQ0LjIwMjk1MywxMjguMzgzNTYgaCAyOC4xMjkxNSIKICAgICAgIGlkPSJwYXRoMjEiIC8+CiAgICA8Y2lyY2xlCiAgICAgICBpZD0iY2lyY2xlMjEiCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZTpub25lO3N0cm9rZS13aWR0aDowLjI2NDU4MyIKICAgICAgIGN4PSI3Mi4zMzIxMDgiCiAgICAgICBjeT0iLTEyOC43NDg3OSIKICAgICAgIHI9IjAuMzk2ODc0OTkiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPHBhdGgKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICBkPSJtIDU4LjIyNTg4LDEyOC41NTYyOCB2IDE1LjExMTU1IgogICAgICAgaWQ9InBhdGgyMiIgLz4KICAgIDx0ZXh0CiAgICAgICB4bWw6c3BhY2U9InByZXNlcnZlIgogICAgICAgc3R5bGU9ImZvbnQtc2l6ZTo0LjIzMzMzcHg7ZmlsbDojMDAwMDAwO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MztzdHJva2Utb3BhY2l0eToxIgogICAgICAgeD0iODAuMDc3ODg4IgogICAgICAgeT0iMTQ0Ljk3NzE0IgogICAgICAgaWQ9InRleHQyMiI+PHRzcGFuCiAgICAgICAgIHNvZGlwb2RpOnJvbGU9ImxpbmUiCiAgICAgICAgIGlkPSJ0c3BhbjIyIgogICAgICAgICBzdHlsZT0ic3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgICB4PSI4MC4wNzc4ODgiCiAgICAgICAgIHk9IjE0NC45NzcxNCI+dDwvdHNwYW4+PC90ZXh0PgogICAgPHRleHQKICAgICAgIHhtbDpzcGFjZT0icHJlc2VydmUiCiAgICAgICBzdHlsZT0iZm9udC1zaXplOjQuMjMzMzNweDtmaWxsOiMwMDAwMDA7c3Ryb2tlOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzO3N0cm9rZS1vcGFjaXR5OjEiCiAgICAgICB4PSI4MC4wNzc4ODgiCiAgICAgICB5PSIxMjkuMzQ5MjEiCiAgICAgICBpZD0idGV4dDIzIj48dHNwYW4KICAgICAgICAgc29kaXBvZGk6cm9sZT0ibGluZSIKICAgICAgICAgaWQ9InRzcGFuMjMiCiAgICAgICAgIHN0eWxlPSJzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICAgIHg9IjgwLjA3Nzg4OCIKICAgICAgICAgeT0iMTI5LjM0OTIxIj50KzE8L3RzcGFuPjwvdGV4dD4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTIzIgogICAgICAgY3g9IjQ0LjIwMjk1IgogICAgICAgY3k9Ii0xNDMuMTk5OTgiCiAgICAgICByPSIxLjgyNjU2ODIiCiAgICAgICBpbmtzY2FwZTpleHBvcnQtZmlsZW5hbWU9ImJhY2t3YXJkIEV1bGVyIHN0ZW5jaWwuc3ZnIgogICAgICAgaW5rc2NhcGU6ZXhwb3J0LXhkcGk9Ijk2IgogICAgICAgaW5rc2NhcGU6ZXhwb3J0LXlkcGk9Ijk2IgogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxLC0xKSIgLz4KICAgIDxjaXJjbGUKICAgICAgIHN0eWxlPSJmaWxsOiMwMDAwMDA7c3Ryb2tlLXdpZHRoOjAuMjY0NTgzIgogICAgICAgaWQ9ImNpcmNsZTI0IgogICAgICAgY3g9IjU4LjIyNTg4IgogICAgICAgY3k9Ii0xNDMuMTk5OTgiCiAgICAgICByPSIxLjgyNjU2ODIiCiAgICAgICB0cmFuc2Zvcm09InNjYWxlKDEsLTEpIiAvPgogICAgPGNpcmNsZQogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBpZD0iY2lyY2xlMjUiCiAgICAgICBjeD0iNzIuMjQ4ODEiCiAgICAgICBjeT0iLTE0My4xOTk5OCIKICAgICAgIHI9IjEuODI2NTY4MiIKICAgICAgIHRyYW5zZm9ybT0ic2NhbGUoMSwtMSkiIC8+CiAgICA8cGF0aAogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBkPSJNIDQ0LjIwMjk1MywxNDMuMjAwMzIgSCA3Mi42OTc0MTYiCiAgICAgICBpZD0icGF0aDI1IiAvPgogICAgPGNpcmNsZQogICAgICAgaWQ9ImNpcmNsZTI2IgogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6bm9uZTtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBjeD0iNzIuNjk3NDE4IgogICAgICAgY3k9Ii0xNDMuMTk5OTgiCiAgICAgICByPSIwLjc5Mzc0OTk5IgogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxLC0xKSIgLz4KICAgIDxwYXRoCiAgICAgICBzdHlsZT0iZmlsbDojMDAwMDAwO3N0cm9rZTojMDAwMDAwO3N0cm9rZS13aWR0aDowLjI2NDU4MztzdHJva2Utb3BhY2l0eToxIgogICAgICAgZD0ibSA0NC4yMDI5NTMsMTQzLjIwMDMyIGggMjguMTI5MTUiCiAgICAgICBpZD0icGF0aDI2IiAvPgogICAgPGNpcmNsZQogICAgICAgaWQ9ImNpcmNsZTI3IgogICAgICAgc3R5bGU9ImZpbGw6IzAwMDAwMDtzdHJva2U6bm9uZTtzdHJva2Utd2lkdGg6MC4yNjQ1ODMiCiAgICAgICBjeD0iNzIuMzMyMTA4IgogICAgICAgY3k9Ii0xNDMuNTY1MjYiCiAgICAgICByPSIwLjM5Njg3NDk5IgogICAgICAgdHJhbnNmb3JtPSJzY2FsZSgxLC0xKSIgLz4KICA8L2c+Cjwvc3ZnPgo=)

The left-hand matrix $\mathbf{A}$ has the exact same tridiagonal sparsity as Backward Euler, delivering $\mathcal{O}(\Delta t^2 + \Delta x^2)$ accuracy and unconditional stability for virtually identical computational cost per step.

differential_equations_module.py

def crank_nicolson_heat(T0, alpha, dx, dt, num_steps):
    T = T0.copy()
    T_history = [T.copy()]
    n = len(T)
    r = alpha * dt / (2.0 * dx**2)
    A = (
        np.diag((1.0 + 2.0 * r) * np.ones(n), 0)
        + np.diag(-r * np.ones(n - 1), -1)
        + np.diag(-r * np.ones(n - 1), 1)
    )
    B = (
        np.diag((1.0 - 2.0 * r) * np.ones(n), 0)
        + np.diag(r * np.ones(n - 1), -1)
        + np.diag(r * np.ones(n - 1), 1)
    )
    A[0, :] = 0.0
    A[0, 0] = 1.0
    A[-1, :] = 0.0
    A[-1, -1] = 1.0
    for _ in range(num_steps):
        b = B @ T
        b[0] = T[0]
        b[-1] = T[-1]
        T = np.linalg.solve(A, b)
        T_history.append(T.copy())
    return np.array(T_history)

exe.display_source(exe.crank_nicolson_heat)
fig = exe.plot_pde_method_comparison()
fig.show()

Summary¶

Implicit time-steppers are invaluable for stiff parabolic PDEs, but each time step requires solving a coupled system of equations. For large 2D/3D meshes, direct factorization ( $\mathbf{L}\mathbf{U}$ ) becomes prohibitively expensive, necessitating iterative linear solvers.

Explicit time-steppers are conditionally stable—which restricts $\Delta t$ —but require no simultaneous equation solving or nonlinear root finding:

No dependence on initial guesses for root finders
Guaranteed per-step update without convergence failures
Trivial parallelization across high-performance computing (HPC) nodes
Minimal memory footprint (no system matrix assembly)

When does the cost of implicit linear/nonlinear iterations outweigh the small time-step restriction of an explicit method? It depends on the problem:

General multiphysics solvers (e.g., COMSOL, ANSYS) typically default to implicit time-stepping to handle coupled stiff physics over long time scales.
Specialized short-duration or massive-mesh solvers often prefer explicit methods that advance millions of degrees of freedom without matrix assembly—such as LS-DYNA (high-velocity impact, crash, and drop tests) and explicit CFD packages (e.g., FLOW-3D).