Diffusion Equation
The dispersionrelationpreserving (DRP) scheme is well suited to compute numerical solutions of diffusion equations. To illustrate this, consider computing the solution of the initial value problem governed by the heat equation as follows:
дU f д2U д2u
dt д X2 + д y2
at t = 0, u = f (x, y).
Before implementing the DRP scheme, it is advantageous to cast (D1) into a system of firstorder equations. This may be done by introducing two auxiliary variables defined by
(D1) becomes
д U / д V д W дt дX^~ ду2
Now, the system of equations (D3) to (D5) is equivalent to (D1).
Upon applying the DRP scheme to (D3) to (D5), the finite difference system is
1 3 v(n) = — a U(n) Vl, m = Дх ajUl+j, m 
(D6) 
j=3 

13 W(n) = 1 a U(n) Wl, m = Ду ajUl, m+j y j=3 
(D7) 
E(n) = K h a 1 v(n) i 1 w(n) lm =3 aj AxVl+j, m + AyWlm+j 
(D8) 
3 u’,t=um+it c bjE’m j 
(D9) 
j=0 
The time step At that can be used in (D9) is constrained by stability consideration. To determine the largest time step without triggering numerical instability, the first step is to find the dispersion relation of the system of finite difference equations (D6) to (D9). This is done by generalizing the discretized system to a set of finite difference equations with continuous variables and then take their FourierLaplace transforms. This leads to (a Ф will be used to denote the FourierLaplace transform of Ф ) the following:
v — iaU 
(D10) 
w — ів U 
(D11) 
Ё — K(iav + івії) 
(D12) 
— f(a^) p і(&Н———— — Ё, 2n 
(D13) 
where f (a, в) is the Fourier transform of the initial condition f(x, y). a and в are the Fourier transform variables in the x and y directions, and ы is the Laplace transform variable. The solution of (D10) to (D13) is
і f (а, в)
U — 2 ♦
2п ы + ік (a2 + в )
The dispersion relation is given by setting the denominator of (D14) to zero, i. e.,
ы — – ік(а2 + в2). (D15)
Note that ы is purely imaginary for real a and в. According to the stability diagram of the DRP scheme (Figure 3.2) the stability limit along the negative imaginary axis of the rnAt plane is &At  < 0.29 . Thus, by dispersion relation (D15), in order to maintain numerical stability, it is required that
к(а2 + в2 )At < 0.29. (D16)
(D16) may be rewritten as
For the 7point stencil DRP scheme, a Ax and в Ay are less than 1.75. Thus, by replacing a Ax and в Ay by 1.75 in (D17), the largest permissible time step, AtMax, is given by
It is a fact that diffusion current always flows outward from the source. Mathematically speaking, this means that the solution at a point, away from the source depends only on the solution at points closer to the source. Therefore, at the boundary of a finite computation domain, the appropriate boundary treatment is to use backward difference stencils in the boundary region. There is no need to impose special boundary conditions.
As a simple numerical example, consider the case of an initial Gaussian pulse,
i. e.,
where b is the halfwidth of the pulse. The exact solution is
Figure D1 shows the result of a computation using the 7point stencil DRP scheme for the case b = 3. The length scale is Ax = Ay. The time scale is x^ . A 100 x 100 mesh was used in the computation. The numerical solutions at t = 20, 40, 80, 160, and 300 are displayed in the figure. The exact solution is shown by a dotted line. However, it cannot be easily seen because it differs from the numerical solution by less than the thickness of the lines of the figure. The numerical solution is extremely accurate.
Leave a reply