Interactive Lecture
Solving PDEs Using Finite Difference Methods
by Dr Arnaldo Delli Carri
Partial Differential Equations (PDEs) describe phenomena involving multiple independent variables and their rates of change. They govern heat transfer, wave propagation, and fluid dynamics.
PDEs govern continuous processes in space and time. Analytical solutions exist only for simple cases, so numerical methods — finite differences, finite elements, finite volumes — are how real engineering problems get solved.
PDEs are classified into three main types based on their mathematical properties:
Parabolic: Heat/Diffusion Equation
\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\]
Examples: heat conduction, mass diffusion
Hyperbolic: Wave Equation
\[\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}\]
Examples: mechanical vibrations, acoustic waves, electromagnetic waves
Elliptic: Laplace/Poisson Equation
\[\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)\]
Examples: steady-state heat conduction, electrostatics, potential flow
The finite difference method approximates derivatives using discrete grid points. For the heat equation, we discretise space and time:
Spatial Discretisation: Central difference for second derivative
\[\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\]
Temporal Discretisation: Forward difference for time derivative
\[\frac{\partial u}{\partial t} \approx \frac{u_i^{n+1} - u_i^n}{\Delta t}\]
Two main approaches exist for solving the discretised equations:
Explicit Method (FTCS - Forward Time, Central Space):
\[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)\]
Simple to implement, but requires stability condition:
\[r = \frac{\alpha \Delta t}{\Delta x^2} \leq \frac{1}{2}\]
Implicit Method (Crank-Nicolson):
\[u_i^{n+1} - \frac{r}{2}(u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}) = u_i^n + \frac{r}{2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)\]
Unconditionally stable but requires solving a system of equations at each time step.
The Courant-Friedrichs-Lewy (CFL) condition ensures numerical stability:
For Heat Equation (Parabolic):
\[r = \frac{\alpha \Delta t}{\Delta x^2} \leq \frac{1}{2} \quad \text{(explicit method)}\]
For Wave Equation (Hyperbolic):
\[CFL = \frac{c \Delta t}{\Delta x} \leq 1\]
These conditions ensure that the numerical domain of dependence includes the physical domain of dependence.
Truncation Error: The error introduced by approximating derivatives with finite differences. For central difference in space and forward in time: O(Δt) + O(Δx²)
Convergence: As Δx → 0 and Δt → 0, the numerical solution approaches the exact solution, provided the stability condition is satisfied.
This simulation solves the 1D heat equation using the explicit finite difference method. Watch how heat diffuses from the initial temperature distribution over time.
Numerical PDE methods show up in these engineering fields:
| Concept | Equation | Description |
|---|---|---|
| Heat Equation (1D) | \(\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\) | Parabolic PDE for diffusion |
| Wave Equation | \(\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}\) | Hyperbolic PDE for propagation |
| Laplace Equation | \(\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0\) | Elliptic PDE for steady-state |
| Central Difference (2nd) | \(\frac{\partial^2 u}{\partial x^2} \approx \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\) | Spatial discretisation O(Δx²) |
| Explicit Method Update | \(u_i^{n+1} = u_i^n + r(u_{i+1}^n - 2u_i^n + u_{i-1}^n)\) | Forward Euler time integration |
| Stability Parameter | \(r = \frac{\alpha \Delta t}{\Delta x^2} \leq \frac{1}{2}\) | CFL condition for stability |
For a 1D heat equation solver with α = 0.001 m²/s and spatial step Δx = 0.02 m, what is the maximum time step Δt that ensures stability using the explicit method? If you want to simulate 10 seconds, how many time steps are required?
Show Solution
Step 1: Apply stability condition
\[r = \frac{\alpha \Delta t}{\Delta x^2} \leq \frac{1}{2}\]
\[\Delta t \leq \frac{\Delta x^2}{2\alpha} = \frac{(0.02)^2}{2 \times 0.001} = \frac{0.0004}{0.002} = 0.2 \text{ s}\]
Step 2: Calculate number of time steps
For safety, use Δt = 0.15 s (r = 0.375 < 0.5)
\[N = \frac{t_{final}}{\Delta t} = \frac{10}{0.15} = 66.7 \approx 67 \text{ time steps}\]
Answer: Δt_max = 0.2 s; using a safety margin of Δt = 0.15 s (r = 0.375), N ≈ 67 time steps are needed to simulate 10 s.
Note: Using a smaller time step (e.g., Δt = 0.1 s) would require more steps (100) but provide better accuracy and stability margin.
The central difference approximation for the second derivative has a truncation error. Using Taylor series expansion, show that the error is O(Δx²). What does this mean for solution accuracy?
Show Solution
Taylor Series Expansion:
\[u(x + \Delta x) = u(x) + u'(x)\Delta x + \frac{u''(x)}{2}\Delta x^2 + \frac{u'''(x)}{6}\Delta x^3 + O(\Delta x^4)\]
\[u(x - \Delta x) = u(x) - u'(x)\Delta x + \frac{u''(x)}{2}\Delta x^2 - \frac{u'''(x)}{6}\Delta x^3 + O(\Delta x^4)\]
Add the two equations:
\[u(x + \Delta x) + u(x - \Delta x) = 2u(x) + u''(x)\Delta x^2 + O(\Delta x^4)\]
Solve for u''(x):
\[u''(x) = \frac{u(x + \Delta x) - 2u(x) + u(x - \Delta x)}{\Delta x^2} + O(\Delta x^2)\]
Answer: The truncation error of the central difference approximation is O(Δx²). The error decreases quadratically with grid refinement — halving Δx reduces the error by a factor of 4.
Compare the computational cost and stability of explicit vs implicit methods. If an explicit method requires Δt = 0.001 s for stability, but an implicit method is stable for any Δt, what is the breakeven point where implicit becomes more efficient?
Show Solution
Explicit Method:
• Cost per time step: O(N) operations (N grid points)
• Δt limited by stability: Δt = 0.001 s
• For 1 second simulation: 1000 time steps
• Total cost: 1000 × N = 1000N operations
Implicit Method:
• Cost per time step: O(N) for tridiagonal solver
• No stability limit: can use Δt = 0.1 s (100× larger)
• For 1 second simulation: 10 time steps
• Total cost: 10 × (5N) ≈ 50N operations (factor of 5 for solver)
Step 3: Breakeven time step
Total cost is equal when (N/Δt_explicit) = (5N/Δt_implicit), i.e. when:
\[\Delta t_{implicit} = 5 \times \Delta t_{explicit} = 5 \times 0.001 = 0.005 \text{ s}\]
Once the implicit method can take a time step larger than 0.005 s, its five-fold per-step cost is repaid by taking fewer steps.
Answer: The breakeven time step is Δt = 0.005 s (5× the explicit stability limit). At Δt_implicit = 0.1 s, the implicit method costs ≈50N operations against ≈1000N for explicit — about 20× more efficient for this case.