User:Compsonheir/Numerical diffusion

In numerical analysis, computer simulations of continua (such as fluids, solids, or electromagnetic waves) can exhibit a higher diffusivity than the true medium; this is called numerical diffusion. In some instances, numerical diffusion is a nuisance because the object of study is not diffusive at all, such as an ideal fluid. However, numerical diffusion is sometimes introduced on purpose for its stabilizing and smoothing effects.

In addition to numerical diffusion, some discretizations of differential equations introduce numerical dispersion instead. Oftentimes there is a tradeoff between the two effects.

The modified differential equation

One can ascertain whether a numerical PDE solution exhibits numerical diffusion or dispersion by finding a different differential equation to which it is a more accurate solution. If that modified equation is diffusive, then the approximation scheme exhibits numerical diffusion.

An approximation with artificial diffusion

Consider the one-dimensional advection equation

${\frac {\partial u}{\partial t}}+c{\frac {\partial u}{\partial x}}=0$

for the unknown function u. This problem has the analytical solution $u=u_{0}(x-ct)$, where u₀ is the initial condition; this makes it very easy to use as a benchmark. Let U be an approximation to u computed using the finite difference method, taking the spatial mesh width to be Δx and the timestep to be Δt, with j indexing the spatial coordinate and n the time coordinate.

First, consider the upwind method:

${\frac {U_{j}^{n+1}-U_{j}^{n}}{\Delta t}}=-c{\frac {U_{j}^{n}-U_{j-1}^{n}}{\Delta x}}$

It is known that, provided Δx/Δt > c, the upwind method is accurate to order O(Δx). However, by appealing to Taylor's theorem, one can show that U is a more accurate solution of another partial differential equation, the modified equation. Specifically, U is an O(Δx²)-accurate approximation to the solution of

${\frac {\partial u}{\partial t}}+c{\frac {\partial u}{\partial x}}={\frac {1}{2}}c\Delta x\left(1-{\frac {c\Delta t}{\Delta x}}\right){\frac {\partial ^{2}u}{\partial x^{2}}},$

which is now an advection-diffusion equation with diffusion coefficient

$\kappa ={\frac {1}{2}}c\Delta x\left(1-{\frac {c\Delta t}{\Delta x}}\right).$

The character of the modified equation tells us a lot about the nature of our numerical approximation. The presence of the diffusive terms means that U will appear smoother and more spread out over time than the true solution u, which is especially pronounced in the propagation of shock waves.

An approximation with artificial dispersion

Now consider a different discretization, the Lax-Wendroff method:

${\frac {U_{j}^{n+1}-U_{j}^{n}}{\Delta t}}=-{\frac {c}{2}}{\frac {U_{j+1}^{n}-U_{j-1}^{n}}{2\Delta x}}+c^{2}\Delta t{\frac {U_{j-1}^{n}-2U_{j}^{n}+U_{j+1}^{n}}{2\Delta x^{2}}}.$

Using Taylor's theorem again, the same analysis as for the upwind method shows that the Lax-Wendroff scheme is an O(Δx³)-accurate solution of the equation

${\frac {\partial u}{\partial t}}+c{\frac {\partial u}{\partial x}}=-{\frac {c\Delta x^{2}}{6}}\left(1-\left({\frac {c\Delta t}{\Delta x}}\right)^{2}\right){\frac {\partial ^{3}u}{\partial x^{3}}}.$

This equation exhibits dispersion rather than diffusion, so that different wavelengths have difference group velocities. In order to ascertain the behavior of the numerical method, we can compute the dispersion relation:

$\omega =ck-{\frac {c\Delta x^{2}}{6}}\left(1-\left({\frac {c\Delta t}{\Delta x}}\right)^{2}\right)k^{3}.$

For small k, the frequency ω is approximately linear in k, as it should be. However, for large k, the group velocity moves in the opposite direction as the wave packet should be.

Numerical dispersion relation

An alternative approach to finding the modified differential equation is to find the approximation's numerical dispersion relation. In this method, one treats the numerical solution as a superposition of waves in a manner akin to Fourier analysis. The same approach is used in von Neumann stability analysis.

A diffusive approximation

Consider again the upwind discretization

${\frac {U_{j}^{n+1}-U_{j}^{n}}{\Delta t}}=-c{\frac {U_{j}^{n}-U_{j-1}^{n}}{\Delta x}}$

and suppose we had some function of the form $\phi _{j}^{n}=e^{i(kj\Delta x-\omega n\Delta t)}$ which solved the difference equation. Given a specific wave number k, what must ω be?