Discretization Methods

The finite difference method (FDM) replaces derivatives by algebraic differences of function values at grid points, spaced $h=\Delta x$ apart. The basis is the Taylor series.

Taylor expansions about $x_i$

$$f_{i+1}=f_i+h f'_i+\frac{h^2}{2!}f''_i+\frac{h^3}{3!}f'''_i+\cdots$$ $$f_{i-1}=f_i-h f'_i+\frac{h^2}{2!}f''_i-\frac{h^3}{3!}f'''_i+\cdots$$

Forward difference

Rearrange the first series for $f'_i$:

$$f'_i=\frac{f_{i+1}-f_i}{h}-\frac h2f''_i-\frac{h^2}{6}f'''_i-\cdots$$ $$\boxed{f'_i\approx\frac{f_{i+1}-f_i}{h}},\qquad \text{error}=-\frac h2f''_i=O(h)$$

Backward difference

Rearrange the second series:

$$f'_i=\frac{f_i-f_{i-1}}{h}+\frac h2f''_i-\cdots$$ $$\boxed{f'_i\approx\frac{f_i-f_{i-1}}{h}},\qquad \text{error}=+\frac h2f''_i=O(h)$$

Central difference

Subtract the second series from the first; the even-order terms cancel:

$$f_{i+1}-f_{i-1}=2hf'_i+\frac{h^3}{3}f'''_i+\cdots$$ $$\boxed{f'_i\approx\frac{f_{i+1}-f_{i-1}}{2h}},\qquad \text{error}=-\frac{h^2}{6}f'''_i=O(h^2)$$

Summary

Halving $h$ halves the error of the first-order schemes but reduces the central-difference error to one quarter. Forward and backward differences are used at boundaries and in upwind schemes, where only one side is available.

Exact: $f'(1)=\cos 1=0.54030$.

Values needed: $\sin0.9=0.783327$, $\sin0.95=0.813416$, $\sin1.0=0.841471$, $\sin1.05=0.867423$, $\sin1.1=0.891207$.

h = 0.1

• Forward: $\dfrac{\sin1.1-\sin1.0}{0.1}=\dfrac{0.891207-0.841471}{0.1}=0.49736$
• Backward: $\dfrac{0.841471-0.783327}{0.1}=0.58144$
• Central: $\dfrac{0.891207-0.783327}{0.2}=0.53940$

h = 0.05

• Forward: $\dfrac{0.867423-0.841471}{0.05}=0.51904$
• Backward: $\dfrac{0.841471-0.813416}{0.05}=0.56111$
• Central: $\dfrac{0.867423-0.813416}{0.1}=0.54008$

Errors

Order of accuracy

Ratio of errors when $h$ is halved:

• Forward: $0.04294/0.02126=2.02\approx2^1$, so first order.
• Backward: $0.04114/0.02081=1.98\approx2^1$, so first order.
• Central: $0.00090/0.000225=4.0=2^2$, so second order.

Here error means (numerical $-$ exact). The signs agree with theory: the forward formula equals $f'+\tfrac h2f''$, and since $f''=-\sin x<0$ its error is negative; the backward error has the opposite sign; the central error is far smaller.

Answer: Central difference is best (0.5394 and 0.5401 against exact 0.5403); forward and backward are $O(h)$ and central is $O(h^2)$.

Derivation

Add the Taylor series of $f_{i+1}$ and $f_{i-1}$:

$$f_{i+1}+f_{i-1}=2f_i+h^2f''_i+\frac{h^4}{12}f''''_i+\cdots$$

Solve for $f''_i$:

$$\boxed{f''_i\approx\frac{f_{i+1}-2f_i+f_{i-1}}{h^2}},\qquad \text{error}=-\frac{h^2}{12}f''''_i=O(h^2)$$

It is second-order accurate, and it uses three points (stencil $i-1,\ i,\ i+1$). It is the basic formula for diffusion terms.

Numerical for $f=x^4$ at $x=1.2$

Exact: $f''=12x^2=12\times1.44=17.28$.

h = 0.2: $f(1.0)=1.0$, $f(1.2)=2.0736$, $f(1.4)=3.8416$

$$f''\approx\frac{1.0-2(2.0736)+3.8416}{0.2^2}=\frac{0.6944}{0.04}=17.36$$

h = 0.1: $f(1.1)=1.4641$, $f(1.3)=2.8561$

$$f''\approx\frac{1.4641-2(2.0736)+2.8561}{0.01}=\frac{0.173}{0.01}=17.30$$

The error falls by a factor of 4 when $h$ is halved, confirming second-order accuracy. It also matches theory: $\dfrac{h^2}{12}f''''=\dfrac{0.04}{12}\times24=0.08$ for $h=0.2$.

Answer: $f''(1.2)\approx17.36$ ($h=0.2$) and $17.30$ ($h=0.1$); exact 17.28.

Truncation error and order

The truncation error is the part of the Taylor series dropped when a derivative is replaced by a difference formula. If the leading dropped term is proportional to $h^p$, the scheme is of order $p$, written $O(h^p)$. For small $h$, the error behaves as $E\approx Ch^p$, so halving $h$ reduces the error by a factor $2^p$.

The order is found from two grids by

$$p=\frac{\ln(E_1/E_2)}{\ln(h_1/h_2)}$$

First order against second order

Consistency, stability, convergence

• Consistency: the truncation error tends to zero as $h,\Delta t\to0$; the difference equation then approaches the original PDE.
• Stability: errors (round-off) introduced during computation do not grow without bound. Explicit schemes have conditions such as $r=\alpha\Delta t/\Delta x^2\le0.5$ and $C=u\Delta t/\Delta x\le1$.
• Convergence: the numerical solution tends to the exact solution as the mesh is refined.

Lax equivalence theorem: for a well-posed linear problem, a consistent scheme converges if and only if it is stable.

Other numerical errors include round-off error (finite computer precision) and iteration (convergence) error. A mesh-independence test uses systematic refinement to estimate discretisation error.

Governing equation and discretisation

Steady 1D conduction with generation:

$$k\frac{d^2T}{dx^2}+\dot q=0$$

Central difference at an interior node $i$:

$$k\frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}+\dot q=0\ \Rightarrow\ -T_{i-1}+2T_i-T_{i+1}=\frac{\dot q\Delta x^2}{k}$$

 x=0     0.025    0.05    0.075    0.1
 (0)------(1)------(2)------(3)------(4)
 100 C                               200 C

Data: $\Delta x=0.1/4=0.025$ m.

$$\frac{\dot q\Delta x^2}{k}=\frac{10^6\times(0.025)^2}{20}=31.25\ \text{K}$$

with $T_0=100$ and $T_4=200$.

Nodal equations

• Node 1: $2T_1-T_2=31.25+100=131.25$
• Node 2: $-T_1+2T_2-T_3=31.25$
• Node 3: $-T_2+2T_3=31.25+200=231.25$

In matrix form:

$$\begin{bmatrix}2&-1&0\\-1&2&-1\\0&-1&2\end{bmatrix}\begin{bmatrix}T_1\\T_2\\T_3\end{bmatrix}=\begin{bmatrix}131.25\\31.25\\231.25\end{bmatrix}$$

Solution (elimination)

From node 1: $T_1=(131.25+T_2)/2$. From node 3: $T_3=(231.25+T_2)/2$. Substitute in node 2:

$$-\frac{131.25+T_2}{2}+2T_2-\frac{231.25+T_2}{2}=31.25\ \Rightarrow\ T_2-181.25=31.25$$ $$T_2=212.5\ ^\circ\text{C},\quad T_1=\frac{131.25+212.5}{2}=171.875\ ^\circ\text{C},\quad T_3=\frac{231.25+212.5}{2}=221.875\ ^\circ\text{C}$$

Exact solution

$$T(x)=T_0+(T_L-T_0)\frac xL+\frac{\dot q}{2k}x(L-x)$$

At $x=0.025$: $100+25+25000\times0.025\times0.075=171.875$. Similarly 212.5 and 221.875 at 0.05 and 0.075.

The two agree exactly because the exact solution is a quadratic, and the central difference formula for $d^2T/dx^2$ has an error proportional to $T''''$, which is zero here. The maximum temperature lies slightly beyond node 3, between 0.075 and 0.1 m.

Answer: $T_1=171.9$ °C, $T_2=212.5$ °C, $T_3=221.9$ °C (equal to the exact values).

(a) FTCS scheme

Use a forward difference in time and a central difference in space:

$$\frac{T_i^{n+1}-T_i^n}{\Delta t}=\alpha\frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}$$ $$T_i^{n+1}=T_i^n+r\,(T_{i+1}^n-2T_i^n+T_{i-1}^n),\qquad r=\frac{\alpha\Delta t}{\Delta x^2}$$

The scheme is first order in time and second order in space. By von Neumann analysis the amplification factor is $G=1-4r\sin^2(\theta/2)$; $|G|\le1$ requires

$$r\le\frac12$$

(equivalently the coefficient of $T_i^n$, $1-2r$, must not be negative.)

(b) Maximum time step

$$\Delta t_{max}=\frac{0.5\,\Delta x^2}{\alpha}=\frac{0.5\times(0.02)^2}{9.7\times10^{-5}}=2.06\ \text{s}$$

(c) Two steps with $\Delta t=1$ s

$$r=\frac{9.7\times10^{-5}\times1}{(0.02)^2}=0.2425\quad(<0.5,\ \text{stable})$$

Step 1 (from 100, 20, 20, 20, 20, 0):

• $T_1=20+0.2425(20-40+100)=20+19.40=39.40$
• $T_2=20+0.2425(20-40+20)=20.00$
• $T_3=20+0.2425(20-40+20)=20.00$
• $T_4=20+0.2425(0-40+20)=20-4.85=15.15$

Step 2 (from 100, 39.40, 20, 20, 15.15, 0):

• $T_1=39.40+0.2425(20-78.8+100)=39.40+9.99=49.39$
• $T_2=20+0.2425(20-40+39.40)=20+4.70=24.70$
• $T_3=20+0.2425(15.15-40+20)=20-1.18=18.82$
• $T_4=15.15+0.2425(0-30.3+20)=15.15-2.50=12.65$

Heat diffuses in from the hot end, and node 4 cools towards the 0 °C end.

Answer: $r\le0.5$, $\Delta t_{max}=2.06$ s; after 2 s the interior nodes are 49.39, 24.70, 18.82, 12.65 °C.

For $\partial T/\partial t=\alpha\,\partial^2T/\partial x^2$, let $r=\alpha\Delta t/\Delta x^2$ and $\delta^2T_i=T_{i+1}-2T_i+T_{i-1}$.

Discretised forms

• Explicit (FTCS): $T_i^{n+1}=T_i^n+r\,\delta^2T_i^{\,n}$. The right side uses only old values, so each new value is found directly.
• Fully implicit (BTCS): $T_i^{n+1}=T_i^n+r\,\delta^2T_i^{\,n+1}$, or $-rT_{i-1}^{n+1}+(1+2r)T_i^{n+1}-rT_{i+1}^{n+1}=T_i^n$. This needs a tridiagonal system solved each step (Thomas algorithm).
• Crank-Nicolson: average of the two, $T_i^{n+1}=T_i^n+\dfrac r2\left(\delta^2T_i^{\,n}+\delta^2T_i^{\,n+1}\right)$. Also tridiagonal.

Comparison

Remarks

• With the explicit scheme, halving $\Delta x$ forces $\Delta t$ to be one quarter, so cost rises about 8 times in 1D.
• Implicit methods are used in most commercial CFD solvers because stable large time steps are possible, especially for steady-state solutions and stiff problems.
• Unconditional stability does not mean accuracy: too large a $\Delta t$ still gives a poor transient.

Courant number

$$C=\frac{u\Delta t}{\Delta x}=\frac{2\times0.04}{0.1}=0.8$$

The physical shift in one step is $u\Delta t=0.08$ m, i.e. 0.8 of a cell, so the exact solution moves the profile to the right by 0.8 cell.

(a) First-order upwind (flow in +x, use the upstream node)

$$\phi_i^{n+1}=\phi_i^n-C(\phi_i^n-\phi_{i-1}^n)$$

• $\phi_1=300-0.8(300-300)=300.0$
• $\phi_2=320-0.8(320-300)=304.0$
• $\phi_3=360-0.8(360-320)=328.0$
• $\phi_4=400-0.8(400-360)=368.0$

(b) Central difference (FTCS)

$$\phi_i^{n+1}=\phi_i^n-\frac C2(\phi_{i+1}^n-\phi_{i-1}^n)$$

• $\phi_1=300-0.4(320-300)=292.0$
• $\phi_2=320-0.4(360-300)=296.0$
• $\phi_3=360-0.4(400-320)=328.0$
• $\phi_4=400-0.4(400-360)=384.0$

Comments

• Upwind results stay between the old minimum (300) and maximum (400): the scheme is bounded and stable for $C\le1$.
• The central scheme gives 292 at node 1, below the lowest old value 300 (an unphysical undershoot). Forward Euler with central space differencing is unconditionally unstable for pure convection, so these oscillations grow with time.
• Upwind is first order. Its truncation error contains $\dfrac{u\Delta x}{2}(1-C)\dfrac{\partial^2\phi}{\partial x^2}$, which acts like an extra diffusion (numerical or false diffusion) with coefficient

$$\Gamma_{num}=\frac{u\Delta x}{2}(1-C)=\frac{2\times0.1}{2}(1-0.8)=0.02\ \text{m}^2/\text{s}$$

that smears sharp gradients. For $C=1$ it vanishes and the scheme is exact in 1D.

• Remedies: finer grids, higher-order bounded schemes (QUICK, TVD, second-order upwind) with limiters.

Answer: $C=0.8$; upwind 300, 304, 328, 368; central 292, 296, 328, 384; numerical diffusion 0.02 m²/s.

Why FVM is preferred

1. It is derived from the integral form, so mass, momentum and energy are conserved locally in each control volume, even on a coarse mesh.
2. It works on unstructured meshes, so complex industrial geometry can be meshed automatically.
3. The physical meaning of each term (flux through a face) is clear, which makes boundary conditions and upwind schemes natural.
4. It is efficient in memory and suited to segregated pressure-based solvers such as SIMPLE.

FDM remains useful for teaching, simple domains and high-order schemes, because its formulae follow directly from the Taylor series.

Derivation

Expand $f_1$ and $f_2$ about $x_0$:

$$f_1=f_0+hf'_0+\frac{h^2}{2}f''_0+\frac{h^3}{6}f'''_0$$ $$f_2=f_0+2hf'_0+2h^2f''_0+\frac{4h^3}{3}f'''_0$$

Eliminate $f''_0$: compute $4f_1-f_2$:

$$4f_1-f_2=3f_0+2hf'_0-\frac{2h^3}{3}f'''_0$$ $$\boxed{f'_0\approx\frac{-3f_0+4f_1-f_2}{2h}},\qquad f'_0=\text{approx}+\frac{h^2}{3}f'''_0+\cdots\ \Rightarrow\ O(h^2)$$

The first-order formula is $f'_0\approx(f_1-f_0)/h$. The second-order formula is important for evaluating wall shear stress and wall heat flux without needing a node outside the domain.

Numerical

Data: $f_0=50,\ f_1=46,\ f_2=40$ °C, $h=0.5\times10^{-3}$ m.

First-order:

$$\frac{dT}{dx}\approx\frac{46-50}{0.5\times10^{-3}}=-8000\ \text{K/m}$$

Second-order:

$$\frac{dT}{dx}\approx\frac{-3(50)+4(46)-40}{2\times0.5\times10^{-3}}=\frac{-150+184-40}{10^{-3}}=-6000\ \text{K/m}$$

The successive differences are -4 and -6, so the temperature curve is getting steeper away from the wall (curvature is present). The first-order formula effectively gives the slope at mid-interval (-8000 K/m), so it overestimates the magnitude of the wall slope; the second-order formula extrapolates to the wall itself.

Wall heat flux, $q''=-k\,dT/dx$:

$$q''=-0.6\times(-6000)=3600\ \text{W/m}^2$$

(first-order would give 4800 W/m², about 33% higher.)

Answer: gradient $-8000$ K/m (1st order), $-6000$ K/m (2nd order); $q''=3.6$ kW/m² (flux directed in the +x direction).