Part II

Partial Differential Equations

In this part, we discuss the standard numerical techniques used to integrate partial
differential equations (PDEs). Here we focus on finite difference method.
Chapter 4
Introduction

4.1 Definition, Notation and Classification

A differential equation involving more than one independent variable and its partial
derivatives with respect to those variables is called a partial differential equation
(PDE).
Consider a simple PDE of the form:

∂
u(x, y) = 0.
∂x
This equation implies that the function u(x, y) is independent of x. Hence the general
solution of this equation is u(x, y) = f (y), where f is an arbitrary function of y. The
analogous ordinary differential equation is

du
= 0,
dx
its general solution is u(x) = c, where c is a constant. This example illustrates that
general solutions of ODEs involve arbitrary constants, whereas solutions of PDEs
involve arbitrary functions.
In general, one can classify PDEs with respect to different criterion, e.g.:
• Order;
• Dimension;
• Linearity;
• Initial/Boundary value problem, etc.
By order of PDE we will understand the order of the highest derivative that occurs.
A PDE is said to be linear if it is linear in unknown functions and their derivatives,
with coefficients depending on the independent variables. The independent variables
typically include one or more space dimensions and sometimes time dimension as
well.
For example, the wave equation

49
 ∂ 2 u(x,t) ∂ 2 u(x,t)
= a2
∂t 2 ∂ x2
is a one-dimensional, second-order linear PDE. In contrast, the Fisher Equation of
the form
∂ u(r,t)
= △u(r,t) + u(r,t) − u(r,t)2,
∂t
where r = (x, y) is a two-dimensional, second-order nonlinear PDE.

Linear Second-Order PDEs

For linear PDEs in two dimensions there is a simple classification in terms of the
general equation

auxx + buxy + cuyy + dux + euy + f u + g = 0 , u = u(x, y) ,

where the coefficients a, b, c, d, e, f and g are real and in general can also be func-
tions of x and y. The PDE’s of this type are classified by the value of discriminant
Dλ = b2 − 4a c of the eigenvalue problem for the matrix
 
a b/2
A= ,
b/2 c

build from the coefficients of the highest derivatives. A simple classification is

shown on the following table [10, 6]:
Dλ Typ Eigenvalues
Dλ < 0 elliptic the same sign
Dλ > 0 hyperbolic different signs
Dλ = 0 parabolic zero is an eigenvalue
For instance, the Laplace equation for the electrostatic potential ϕ in the space
without a charge
∂ 2ϕ ∂ 2ϕ
+ 2 =0
∂ x2 ∂y
is elliptic, as a = c = 1, b = 0, Dλ = −4 < 0. In general, elliptic PDEs describe
processes that have already reached steady state, and hence are time-independent.
One-dimensional wave equation for some amplitude A(x,t)

∂ 2A 2∂ A
2
− v =0
∂ t2 ∂ x2
with the positive dispersion velocity v is hyperbolic (a = 1, b = 0, c = −v2 , Dλ =
4v2 > 0). Hyperbolic PDEs describe time-dependent, conservative processes, such
as convection, that are not evolving toward steady state.
The next example is a diffusion equation for the patricle’s density ρ (x,t)
 ∂ρ ∂ 2ρ
=D 2,
∂t ∂x
where D > 0 is a diffusion coefficient. This equation is called to be parabolic
(a = −D, b = c = 0,Dλ = 0). Parabolic PDEs describe time-dependent, dissipative
processes, such as diffusion, that are evolving toward steady state.
Each of these classes should be investigated separately as different methods are
required for each class. The next point to emphasize is that as all the coefficients of
the PDE can depend on x and y, this classification concept is local.

Initial and Boundary-Value Problems

As it was mentioned above the solution of PDEs involve arbitrary functions. That
is, in order to solve the system in question completely, additional conditions are
needed. These conditions can be given in the form of initial and boundary condi-
tions. Initial conditions define the values of the dependent variables at the initial
stage (e.g. at t = 0), whereas the boundary conditions give the information about
the value of the dependent valiable or its derivative on the boundary of the area of
interest. Typically, one distinguishes
- Dirichlet conditions specify the values of the dependent variables of the boundary
points.
- Neumann conditions specify the values of the normal gradients of the boundary.
- Robin conditions defines a linear combination of the Drichlet and Neumann con-
ditions.
Moreover, it is useful to classify the PDE in question in terms of initial value
problem (IVP) and boundary value problem (BVP).
- Initial value problem: PDE in question describes time evolution, i.e., the initial
space-distribution is given; the goal is to find how the dependent variable propa-
gates in time ( e.g., the diffusion equation).
- Boundary value problem: A static solution of the problem should be found in
some region-and the dependent variable is specified on its boundary ( e.g., the
Laplace equation).

4.2 Finite difference method

Let us consider a one-dimensional PDE for the unknown function u(x,t). One way to
numerically solve the PDE is to approximate all the derivatives by finite differences.
We divide the domain in space using a mesh x0 , x1 , . . . , xN and in time using a mesh
t0 ,t1 , . . . ,tT . Fisrt we assume a uniform partition both in space and in time, so that
the difference between two consecutive space points will be △x and between two
consecutive time points will be △t, i.e.,
 xi = x0 + i△x, i = 0, 1, . . . , M;
t j = t0 + j△t, j = 0, 1, . . . , T ;

The Taylor series method

Consider a Taylor expansion of an analytical function u(x).

∞
△xn ∂ n u ∂ u △x2 ∂ 2 u △x3 ∂ 3 u
u(x + △x) = u(x) + ∑ = u(x) + △x + + + ...
n=1 n! ∂ x ∂x 2! ∂ x2 3! ∂ x3
n

(4.1)
Then for the first derivative one obtains:
∂ u u(x + △x) − u(x) △x ∂ 2 u △x2 ∂ 3 u
= − − − ... (4.2)
∂x △x 2! ∂ x2 3! ∂ x3
If we break the right hand side of the last equation after the first term, for △x ≪ 1
the last equation becomes

∂ u u(x + △x) − u(x) △i u

= + O(△x) = + O(△x) , (4.3)
∂x △x △x

where
△i u = u(x + △x) − u(x) := ui+1 − ui
is called a forward difference. The backward expansion of the function u can be
written as △x ≪ 1 the last equation reads

∂ u △x2 ∂ 2 u △x3 ∂ 3 u
u(x + (−△x)) = u(x) − △x + − + . . ., (4.4)
∂x 2! ∂ x2 3! ∂ x3
so for the first derivative one obtains

∂ u u(x) − u(x − △x) ∇i u

= + O(△x) = + O(△x) , (4.5)
∂x △x △x

where
∇i u = u(x) − u(x − △x) := ui − ui−1
is called a backward difference. One can see that both forward and backward dif-
ferences are of the order O(△x). We can combine these two approaches and derive
a central difference, which yields a more accurate approximation. If we substract
Eq. (4.5) from Eq. (4.3) one obtains

∂u △x3 ∂ 3 u
u(x + △x) − u(x − △x) = 2△x +2 + ..., (4.6)
∂x 3! ∂ x3
what is equivalent to

∂ u u(x + △x) − u(x − △x)

= + O(△x2 ) (4.7)
∂x 2△x

Note that the central difference (4.7) is of the order of O(△2 x).
The second derivative can be found in the same way using the linear combination of
different Taylor expansions. For instance, consider

∂ u (2△x)2 ∂ 2 u (2△x)3 ∂ 3 u
u(x + 2△x) = u(x) + 2△x + + + ... (4.8)
∂x 2! ∂ x2 3! ∂ x3
Substracting from the last equation Eq. (4.1), multiplied by two, one gets the fol-
lowing equation

∂ 2u 3∂ u
3
u(x + 2△x) − 2u(x + △x) = −u(x) + △x2 + △x + ... (4.9)
∂ x2 ∂ x3
Hence one can approximate the second derivative as

∂ 2 u u(x + 2△x) − 2u(x + △x) + u(x)

= + O(△x). (4.10)
∂ x2 △x2

Similarly one can obtain the expression for the second derivative in terms of back-
ward expansion, i.e.,

∂ 2 u u(x − 2△x) − 2u(x − △x) + u(x)

= + O(△x). (4.11)
∂ x2 △x2

Finally, if we add Eqn. (4.3) and (4.5) an expression for the cental second derivative
reads
∂ 2 u u(x + △x) − 2u(x) + u(x − △x)
= + O(△x2 ). (4.12)
∂ x2 △x2
One can see that approximation (4.12) provides more accurate approximation as
(4.10) and (4.11).
In an analogous way one can obtain finite difference approximations to higher or-
der derivatives and differential operators. The coefficients for first three deriva-
tives for the case of forward, backward and central differences are given in Ta-
bles 4.1, 4.2, 4.3.

Mixed derivatives

A finite difference approximations for the mixed partial derivatives can be calculated
in the same way. For example, let us find the central approximation for the derivative
 ui ui+1 ui+2 ui+3 ui+4
△x ∂∂ ux -1 1
△x2 ∂∂ x2u
2
1 -2 1
3 ∂3u
△x ∂ x3 -1 3 -3 1
△x4 ∂∂ x4u
4
1 -4 6 -4 1
Table 4.1 Forward difference quotient, O(△x)

ui−4 ui−3 ui−2 ui−1 ui

△x ∂∂ ux -1 1
△x2 ∂∂ x2u
2
1 -2 1
∂ 3
△x3 ∂ x3u -1 3 -3 1
△x4 ∂∂ x4u
4
1 -4 6 -4 1
Table 4.2 Backward difference quotient, O(△x)

ui−2 ui−1 ui ui+1 ui+2

2△x ∂∂ ux -1 0 1
△x2 ∂∂ x2u
2
1 -2 1
2△x3 ∂∂ x3u
3
-1 2 0 -2 1
△x4 ∂∂ x4u
4
1 -4 6 -4 1
Table 4.3 Central difference quotient, O(△x2 )

   
∂ 2u ∂ ∂u ∂ u(x,y+△y)−u(x,y−△y)
∂ x∂ y = ∂x ∂y = ∂x 2△y + O(△y2 ) =
u(x+△x,y+△y)−u(x−△x,y+△y)−u(x+△x,y−△y)+u(x−△x,y−△y)
= 4△x△y + O(△x2 △y2 ) .

A nonequidistant mesh

In the section above we have considered different numerical approximations for the
derivatives using the equidistant mesh. However, in many applications it is con-
vinient to use a nonequidistant mesh, where the spatial steps fulfill the following
rule:
△xi = α △xi−1 .
If α = 1 the mesh is said to be equidistant. Let us now calculate the first derivative
of the function u(x) of the second-order accurance:

∂ u (α △x)2 ∂ 2 u (α △x)3 ∂ 3 u
u(x + α △x) = u(x) + α △x + + + ... (4.13)
∂x 2! ∂ x2 3! ∂ x3
Adding the last equation with Eq. (4.4) multiplied by α one obtains the expression
for the second derivative
 ∂ 2 u u(x + α △x) − (1 + α )u(x) + α u(x − △x)
= + O(△x) (4.14)
∂ x2 2 α (α + 1)△x
1 2

Substitution of the last equation into Eq. (4.4) yields

∂ u u(x + α △x) − (1 − α 2)u(x) − α 2 u(x − △x)

= + O(△x2 ) . (4.15)
∂x α (α + 1)△x

4.3 von Neumann stability analysis

One of the central questions arising by numerical treatment of the problem in ques-
tion is stability of the numerical scheme we are interested in [7]. An algorithm for
solving an evolutionary partial differential equation is said to be stable, if the nu-
merical solution at a fixed time remains bounded as the step size goes to zero, so
the perturbations in form of, e.g., rounding error does not increase in time. Unfortu-
nately, there are no general methods to verify the numerical stability for the partial
differential equations in general form, so one restrict oneself to the case of linear
PDE’s. The standard method for linear PDE’s with periodic boundary conditions
was proposed by John von Neumann [4, 2] and is based on the representation of the
rounding error in form of the Fouirer series.
In order to illustrate the procedure, let us introduce the following notation:

u j+1 = T [u j ]. (4.16)

Here T is a nonlinear operator, depending on the numerical scheme in question.

The successive application of T results in a sequence of values

u0 , u1 , u2 , . . . ,

that approximate the exact solution of the problem. However, at each time step we
add a small error ε j , i.e., the sequence above reads

u0 + ε 0 , u1 + ε 1 , u2 + ε 2 , . . . ,

where ε j is a cumulative rounding error at time t j . Thus we obtain

u j+1 + ε j+1 = T (u j + ε j ) . (4.17)

After linearization of the last equation (we suppose that Taylor expansion of T is
possible) the linear equation for the perturbation takes the form:

∂ T (u j ) j
ε j+1 = ε := Gε j . (4.18)
∂uj
This equation is called the error propagation law, whereas the linearization matrix G
is said to be an amplification matrix [5]. Now, the stability of the numerical scheme
in question depends on the eigenvalues gµ of the matrix G. In other words, the
scheme is stable if and only if

|gµ | ≤ 1 ∀µ

Now the question is how this information can be used in practice. The first point to
emphasize is that in general one deals with the u(xi ,t j ) := uij , so one can write

εij+1 = ∑ Gii′ εij′ , (4.19)

i′

where
∂ T (u j )i
Gii′ = .
∂ uij′

Futhermore, the spatial variation of εij (rounding error at the time step t j at the point
xi ) can be expanded in a finite Fourier series in the intreval [0, L]:

εij = ∑ eikxi ε̃ j (k), (4.20)

k

where k is the wavenumber and ε̃ j (k) are the Fourier coefficients. Since the rounding
error tends to grow or decay exponentially with time, it is reasonable to assume that
ε̃ j (k) varies exponentially with time, i.e.,

εij = ∑ eω t j eikxi ,
k

where ω is a constant. The next point to emphasize is that the functions eikxi are
eigenfunctions of the matrix G, so the last expansion can be interpreted as the ex-
pansion in eigenfunctions of G. In addition, the equation for the error is linear, so it
is enough to examine the grows of the error of a typical term of the sum. Thus, from
the practical point of view one take the error εi just as
j

εij = eω t j eikxi .

The substitution of this expression into Eq. (4.20) results in the following relation

εij+1 = g(k)εij . (4.21)

That is, one can interpert eikxi as an eigenvector corresponding to the eigenvalue
g(k). The value g(k) is often called an amplification factor. Finally, the stability
criterium is then given by
|g(k)| ≤ 1 ∀k . (4.22)
This criterium is called von Neumann stablity criterium.
Notice that presented stability analysis can be applied only in certain cases. Namely,
the linear PDE in question schould be with constant coefficients and satisfies peri-
odic boundary conditions. In addition, the corresponding difference scheme should
possesses no more than two time levels [9]. However, it is often used in more com-
plicated situations as a good estimation for the step sizes used in the approximation.