Skip to main content

Numerical dispersion, stability; four-point implicit (500 m)

Short answer

Numerical diffusion smears fronts; dispersion causes oscillations; explicit schemes require Courant number ≤1. Using four-point implicit (θ=0.55), Δt=1 h, Δx=500 m, ∂Q/∂t = 6.94×10^-4 m^3/s^2 and ∂Q/∂x = -0.00690 m^2/s.

Describe numerical dispersion, diffusion, and stability of finite difference schemes. The value of flow rate QQ at four points in the space-time grid is shown in the figure below. Determine the value of first-order derivatives ∂Q/∂t\partial Q/\partial t and ∂Q/∂x\partial Q/\partial x by using the four-point implicit method. Given: Δt=1\Delta t = 1 hour, Δx=500\Delta x = 500 m and θ=0.55\theta = 0.55.

[Figure: space-time grid with time levels n and n+1 and distance points i and i+1; QQ at (i, n+1) =75= 75 m3^3/s, at (i+1, n+1) =72= 72 m3^3/s, at (i, n) =73= 73 m3^3/s, at (i+1, n) =69= 69 m3^3/s]

Numerical diffusion

The leading truncation error of an odd-order scheme (e.g. first-order upwind) behaves like a physical diffusion term. For Qt+cQx=0Q_t+cQ_x=0 with backward-space, forward-time:

Qt+cQx=cΔx2(1−C) Qxx+…Q_t+cQ_x=\frac{c\Delta x}{2}(1-C)\,Q_{xx}+\dots

The solution is smeared: peaks are lowered and sharp fronts spread out, although the true equation has no diffusion.

Numerical dispersion

For second-order (central) schemes the leading error term has an odd derivative, ∝Qxxx\propto Q_{xxx}. Different wavelengths then travel at different speeds, so a sharp front produces spurious oscillations (ripples ahead of or behind the front). It is a phase error, with little change in amplitude.

Stability

A scheme is stable if errors (round-off or truncation) do not grow from step to step. For explicit schemes, the von Neumann or Courant-Friedrichs-Lewy condition requires

C=c ΔtΔx≤1C=\frac{c\,\Delta t}{\Delta x}\le1

Implicit schemes (and the four-point scheme with θ≥0.5\theta\ge0.5) are unconditionally stable, but a large Δt\Delta t still reduces accuracy.

Four-point implicit method

∂Q∂t≈(Qin+1−Qin)+(Qi+1n+1−Qi+1n)2Δt\frac{\partial Q}{\partial t}\approx\frac{(Q_i^{n+1}-Q_i^n)+(Q_{i+1}^{n+1}-Q_{i+1}^n)}{2\Delta t} ∂Q∂x≈θ Qi+1n+1−Qin+1Δx+(1−θ) Qi+1n−QinΔx\frac{\partial Q}{\partial x}\approx\theta\,\frac{Q_{i+1}^{n+1}-Q_i^{n+1}}{\Delta x}+(1-\theta)\,\frac{Q_{i+1}^n-Q_i^n}{\Delta x}

Data: Qin+1=75Q_i^{n+1}=75, Qi+1n+1=72Q_{i+1}^{n+1}=72, Qin=73Q_i^n=73, Qi+1n=69Q_{i+1}^n=69; Δt=3600\Delta t=3600 s, Δx=500\Delta x=500 m, θ=0.55\theta=0.55.

Time derivative

∂Q∂t=(75−73)+(72−69)2×3600=57200=6.9444×10−4 m3/s2\frac{\partial Q}{\partial t}=\frac{(75-73)+(72-69)}{2\times3600}=\frac{5}{7200}=6.9444\times10^{-4}\ \text{m}^3/\text{s}^2

Space derivative

∂Q∂x=0.55 72−75500+0.45 69−73500=−0.00330+(−0.00360)=−0.00690 m3/s per m\frac{\partial Q}{\partial x}=0.55\,\frac{72-75}{500}+0.45\,\frac{69-73}{500}=-0.00330+(-0.00360)=-0.00690\ \text{m}^3/\text{s per m}

Answer: ∂Q/∂t=6.9444×10−4\partial Q/\partial t=6.9444\times10^{-4} m3^3/s2^2 and ∂Q/∂x=−0.00690\partial Q/\partial x=-0.00690 m2^2/s.

✓ Answer checked by Claude Opus 5.5 (high effort)

All Computational Techniques in Civil Engineering chapter 5 questions →