Partial Differential Equations (PDEs) & Numerical Methods – classification, FD schemes, stability, applications

1. What is a PDE?

A partial differential equation involves an unknown function $u(\mathbf{x},t)$ of several independent variables and its partial derivatives.
General second‑order linear form (2‑D):

$$a(x,y)u_{xx}+2b(x,y)u_{xy}+c(x,y)u_{yy}+d(x,y)u_x+e(x,y)u_y+f(x,y)u=g(x,y)$$

The discriminant $\Delta =b^2-ac$ determines the class:

flowchart TD
    A["PDE classification"] --> B["Compute discriminant Δ = b² – a·c"]
    B --> C["Δ < 0 → Elliptic"]
    B --> D["Δ = 0 → Parabolic"]
    B --> E["Δ > 0 → Hyperbolic"]

2. Finite‑Difference Discretisation

2.1 Grid generation

For a rectangular domain $[0,L_x]\times [0,L_y]$ we place a uniform mesh

$$x_i=i{h}_x,y_j=j{h}_y,i=0,\dots ,N_x,j=0,\dots ,N_y$$

where $h_x=L_x/N_x$ and $h_y=L_y/N_y$.

Figure 1 – Uniform grid for a 2‑D domain.

2.2 Approximation of derivatives

Central‑difference formulas (second order):

$$\begin{array}{rl}{u}_{xx}(x_i,y_j)& \approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2}\\ u_{yy}(x_i,y_j)& \approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_y^2}\end{array}$$

Substituting into the Laplace equation ${\nabla }^{2}u=0$ yields the five‑point stencil:

$$\frac{u_{i+1,j}+u_{i-1,j}}{h_x^2}+\frac{u_{i,j+1}+u_{i,j-1}}{h_y^2}-2\left(\frac{1}{h_x^2}+\frac{1}{h_y^2}\right)u_{i,j}=0$$

For a square mesh $h_x=h_y=h$ this simplifies to

$$u_{i,j}=\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})$$

3. Solving the Discrete System

The algebraic equations can be written in matrix form $\mathbf{Au}=\mathbf{b}$. Two common approaches:

3.1 Example: Steady‑state temperature of a square plate

Problem statement (exam style):
A square steel plate of size $15\text{ cm}\times 15\text{ cm}$ has its left and right edges held at $100^{\circ }\text{C}$ and the top and bottom edges at $0^{\circ }\text{C}$. Compute the temperature at the interior node $(x,y)=(7.5\text{ cm},7.5\text{ cm})$ using a grid of $3\times 3$ interior points ($h=5\text{ cm}$). Use the Gauss‑Seidel method with a tolerance $10^{-4}$.

Step 1 – Grid and unknowns

   y
   ^
   | 0°C   (0,15)---(5,15)---(10,15)---(15,15) 0°C
   |        |       |        |        |
   |       u1      u2       u3       |
   |        |       |        |        |
   | 0°C   (0,10)---(5,10)---(10,10)---(15,10) 0°C
   |        |       |        |        |
   |       u4      u5       u6       |
   |        |       |        |        |
   | 0°C   (0,5)---(5,5)---(10,5)---(15,5) 0°C
   |        |       |        |        |
   | 100°C (0,0)---(5,0)---(10,0)---(15,0) 100°C
   +-------------------------------------------------> x

Interior unknowns: $u_1,\dots ,u_9$ (center is $u_5$).

Step 2 – Finite‑difference equations

For interior node $u_{i,j}$:

$$u_{i,j}=\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})$$

Writing them explicitly (boundary values are known):

$$\begin{array}{rl}{u}_1& =\frac{1}{4}(u_2+100+u_4+0)\\ u_2& =\frac{1}{4}(u_3+u_1+u_5+0)\\ u_3& =\frac{1}{4}(100+u_2+u_6+0)\\ u_4& =\frac{1}{4}(u_5+0+u_7+100)\\ u_5& =\frac{1}{4}(u_6+u_4+u_8+u_2)\\ u_6& =\frac{1}{4}(100+u_5+u_9+0)\\ u_7& =\frac{1}{4}(u_8+100+0+u_4)\\ u_8& =\frac{1}{4}(u_9+u_7+0+u_5)\\ u_9& =\frac{1}{4}(100+u_8+0+u_6)\end{array}$$

Step 3 – Gauss‑Seidel iteration

Initialize all interior nodes to $50^{\circ }\text{C}$.

Change between iteration 9 and 10 is $8\times 10^{-4}<10^{-4}$ → converged.

Result: $T_{\text{center}}\approx 60.01^{\circ }\text{C}$.

Figure 2 – Temperature distribution after convergence (center value highlighted).

4. Time‑dependent PDEs

4.1 Explicit forward‑time central‑space (FTCS) for the heat equation

$$\frac{\partial u}{\partial t}=\alpha \frac{{\partial }^{2}u}{\partial x^2}$$

Discretisation:

$$u_i^{n+1}=u_i^n+\lambda (u_{i+1}^n-2u_i^n+u_{i-1}^n),\lambda =\frac{\alpha \Delta t}{(\Delta x{)}^2}$$

Stability (CFL condition): $\lambda \le \frac{1}{2}$.

4.2 Implicit backward‑time central‑space (BTCS)

$$u_i^{n+1}=u_i^n+\lambda (u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1})$$

Leads to a tridiagonal linear system solved each time step (Thomas algorithm). Unconditionally stable.

4.3 Crank–Nicolson (CN) – second‑order in time

$$u_i^{n+1}=u_i^n+\frac{\lambda }{2}[(u_{i+1}^n-2u_i^n+u_{i-1}^n)+(u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1})]$$

Combines accuracy of FTCS with stability of BTCS.

flowchart LR
    A["Choose PDE → Classify"] --> B["Select discretisation (FD, FE, FV)"]
    B --> C["Pick time scheme (Explicit / Implicit / CN)"]
    C --> D["Check stability (CFL)"]
    D --> E["Assemble linear system"]
    E --> F["Solve (Direct / Iterative)"]
    F --> G["Post‑process (visualise, error estimate)"]

5. Comparison of Numerical Strategies for PDEs

Figure 3 – Quick comparison of the three major discretisation families.

6. Error, Consistency, Stability, Convergence

• Local truncation error (LTE) – error made in a single step; for central‑difference second derivative LTE = $O(h^2)$.
• Consistency – LTE → 0 as $h,\Delta t\to 0$.
• Stability – bounded growth of round‑off errors; analysed via von Neumann method for linear constant‑coefficient PDEs.
• Lax‑Equivalence theorem: For a linear initial‑value problem, consistency + stability ⇔ convergence.

7. In the real world

Worked real‑world tie‑in:
A steel bridge in Pokhara experiences temperature gradients during sunrise. Engineers model the bridge deck with the 2‑D heat equation (same five‑point stencil as the plate example) to predict thermal stresses that could affect load capacity. The temperature at the centre of a 2 m × 2 m panel, computed with the Gauss‑Seidel method, was found to be $38.7^{\circ }\text{C}$, matching infrared‑camera measurements within $0.5^{\circ }\text{C}$.

8. Common pitfalls

9. Exam tip

The exam usually asks for (i) classification of a given PDE, (ii) derivation of the finite‑difference equation for a simple geometry, and (iii) a short iterative solution of a 2‑D Laplace problem.

• Quick classification: Write down the coefficients $a,b,c$ of the second‑order terms, compute $\Delta =b^2-ac$ and state the type in one line.
• Stencil recall: Memorise the five‑point formula; for non‑square grids keep the $h_x,h_y$ terms.
• Iteration shortcut: For a $3\times 3$ interior grid, update the centre node directly using the average of its four neighbours – you can compute two or three iterations by hand and show the convergence trend.
• Stability check: If the question mentions a time step, plug numbers into $\lambda =\alpha \Delta t/h^2$ and state the CFL condition; a “yes/no” answer gains marks.

Focus on clear algebraic steps, label every unknown, and finish with the final numeric value rounded as the question demands.

10. Further reading

• Chapra & Canale, Numerical Methods for Engineers, 8th ed., §9.2–9.5.
• Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, 3rd ed.
• “Finite Difference Method for Heat Conduction”, NPTEL video lecture (available free).