4.3 Discretization in space: stability and convergence

After dealing with temporal discretization, we now turn to another fundamental question of numerical analysis of initial value problems, which is to find conditions under which the discrete (spatial) approximation u (x,t) N converges to the right solution u(x, t) of the PDE (116View Equation) as N → ∞ and t ∈ [0, T]. The time derivative term is treated formally, as one would treat a source term on the right-hand side, that we do not consider here, for better clarity.

A given spatial scheme of solving the PDE is said to be convergent if any numerical approximation uN (x,t), obtained through this scheme, to the solution u(x,t)

∥PN u − uN∥L2w → 0 as N → ∞. (141 )
Two more concepts are helpful in the convergence analysis of numerical schemes:

4.3.1 Lax–Richtmyer theorem

The direct proof of convergence of a given scheme is usually very difficult to obtain. Therefore, a natural approach is to use the Lax–Richtmyer equivalence theorem: “a consistent approximation to a well-posed linear problem is stable if and only if it is convergent”. Thus, the study of convergence of discrete approximations can be reduced to the study of their stability, assuming they are consistent. Hereafter, we sketch out the proof of this equivalence theorem.

The time-evolution PDE (116View Equation) is approximated by

∂uN-- ∂t = LN uN. (144 )
To show that stability implies convergence, we subtract it from the exact solution (116View Equation)
∂-(u-−-uN-) ∂t = LN (u − uN ) + Lu − LN u,
and obtain, after integration, (the dependence on the space coordinate x is skipped)
∫ t u (t) − u (t) = eLN t[u (0) − u (0)] + eLN(t−s)[Lu (s) − L u(s)]ds. (145 ) N N 0 N
Using the stability property (143View Equation), the norm (2 Lw) of this equation implies
∫ t ∥u (t) − u (t)∥ ≤ C(t)∥u(0) − u (0)∥ + C (t − s)∥Lu (s) − L u(s)∥ds. (146 ) N N 0 N
Since the spatial approximation scheme is consistent and C(t) is a bounded function independent of N, for a given t ∈ [0, T] the left-hand side goes to zero as N → ∞, which proves the convergence.

Conversely, to show that convergence implies stability, we use the triangle inequality to get

|∥ ∥ ∥ ∥| ∥ ∥ 0 ≤ |∥eLN tu∥ − ∥eLtu ∥| ≤ ∥eLNtu − eLtu∥ → 0.
From the well-posedness Lt ∥e u∥ is bounded and therefore L t ∥e N u∥ is bounded as well, independent of N.

The simplest stability criterion is the von Neumann stability condition: if we define the adjoint L∗ of the operator L, using the inner product, with weight w of the Hilbert space

∀(u,v ) ∈ ℋ, (u, Lv)w = (L ∗u,v)w ,
then the matrix representation ∗ LN of ∗ L is also the adjoint of the matrix representation of LN. The operator LN is said to be normal if it commutes with its adjoint ∗ L N. The von Neumann stability condition states that, for normal operators, if there exists a constant K independent of N, such that
∀i,1 ≤ i ≤ N, Re (λi) < K, (147 )
with (λi) being the eigenvalues of the matrix LN, then the scheme is stable. This condition provides an operational technique for checking the stability of normal approximations. Unfortunately, spectral approximations using orthogonal polynomials have, in general, strongly non-normal matrices LN and therefore, the von Neumann condition cannot be applied. Some exceptions include Fourier-based spectral approximations for periodic problems.

4.3.2 Energy estimates for stability

The most straightforward technique for establishing the stability of spectral schemes is the energy method: it is based on choosing the approximate solution itself as a test function in the evaluation of residual (60View Equation). However, this technique only provides a sufficient condition and, in particular, crude energy estimates indicating that a spectral scheme might be unstable can be very misleading for non-normal evolution operators (see the example in Section 8 of Gottlieb and Orszag [94]).

Some sufficient conditions on the spatial operator L and its approximation LN are used in the literature to obtain energy estimates and stability criteria, including:

As an illustration, we now consider a Galerkin method applied to the solution of Equation (116View Equation), in which the operator L is semibounded, following the definition (148View Equation). The discrete solution uN is such that the residual (60View Equation) estimated on the approximate solution uN itself verifies

( ) ∂uN -∂t--− LuN ,uN = 0. (153 ) w
Separating the time derivative and the spatial operator:
1-d- 2 1- ∗ 2dt ∥uN (t)∥w = 2 ((L + L )uN (t),uN (t))w ,
which shows that the “energy”
2 γt 2 ∥uN (t)∥ ≤ e ∥uN (0)∥ (154 )
grows at most exponentially with time. Since uN (t) = eLNtuN (0) for any uN (0), we obtain
∥∥eLN t∥∥ ≤ e12γt, (155 )
which gives stability and therefore convergence, provided that the approximation is consistent (thanks to the Lax–Richtmyer theorem).

4.3.3 Examples: heat equation and advection equation Heat equation
We first study the linear heat equation

∂u ∂2u --- − ---2 = 0, with − 1 < x < 1,t > 0, (156 ) ∂t ∂x
with homogeneous Dirichlet boundary conditions
∀t ≥ 0, u(− 1,t) = u(1,t) = 0, (157 )
and initial condition
∀ − 1 ≤ x ≤ 1, u(x,0) = u0(x ). (158 )
In the semidiscrete approach, the Chebyshev collocation method for this problem (see Section 2.5.3) can de devised as follows: the spectral solution uN (t > 0 ) is a polynomial of degree N on the interval [− 1,1], vanishing at the endpoints. On the other Chebyshev–Gauss–Lobatto collocation points {xk} k=1...N− 1 (see Section 2.4.3), uN (t) is defined through the collocation equations
∂u ∂2u ∀k = 1...N − 1, --(xk, t) − ---2(xk,t) = 0, (159 ) ∂t ∂x
which are time ODEs (discussed in Section 4.1) with the initial conditions
0 ∀k = 0...N, uN (xk,0) = u (xk). (160 )

We will now discuss the stability of such a scheme, with the computation of the energy bound to the solution. Multiplying the k-th equation of the system (159View Equation) by uN (xk,t)wk, where {wk }k=0...N are the discrete weights of the Chebyshev–Gauss–Lobatto quadrature (Section 2.4.3), and summing over k, one gets:

N N 1-d-∑ 2 ∑ ∂2uN-- 2 dt (uN (xk,t)) wk − ∂x2 (xk,t)uN (xk,t)wk = 0. (161 ) k=0 k=0
Boundary points (k = 0, N) have been included in the sum since uN is zero there due to boundary conditions. The product uN × ∂2uN ∕∂x2 is a polynomial of degree 2N − 2, so the quadrature formula is exact
∑N 2 ∫ 1 2 ∂--uN-(x ,t)u (x ,t)w = ∂-uN-(x ,t)u (x ,t)w (x)dx, (162 ) ∂x2 k N k k −1 ∂x2 k N k k=0
and integrating by parts twice, one gets the relation
∫ 1 2 ∫ 1( )2 ∫ 1 2 ∂-uN-(x ,t)u (x ,t)w (x)dx = ∂uN-- w (x)dx − 1- u2 ∂-w-dx. (163 ) −1 ∂x2 k N k −1 ∂x 2 −1 N∂x2
By the properties of the Chebyshev weight
2 ( )2 2 ( ) ∂-w-− 2- ∂w- = 0 and ∂-w- = 1 + 2x2 w5, (164 ) ∂x2 w ∂x ∂x2
it is possible to show that
∫ 1 2 ∫ 1 ∫ 1 2 u2 ∂-w-dx ≤ 3 u2 w5dx ≤ 6 ∂-uN--(xk, t)uN (xk,t)w (x )dx, (165 ) −1 N ∂x2 −1 N − 1 ∂x2
and thus that
∫ ∫ ( ) 1 ∂2uN-- 1- 1 ∂uN-- 2 ∂x2 (xk,t)uN(xk, t)w (x)dx ≥ 4 ∂x w(x)dx ≥ 0. (166 ) −1 −1
Therefore, integrating relation (161View Equation) over the time interval [0,t], one obtains
∑N ∑N (uN (xk,t))2wk ≤ (u0(xk))2wk ≤ 2 max |u0(x )|2. (167 ) x∈[0,1] k=0 k=0
The left-hand side represents the discrete norm of uN (t)2, but since this is a polynomial of degree 2N, one cannot apply the Gauss–Lobatto rule. Nevertheless, it has been shown (see, e.g., Section 5.3 of Canuto et al. [57Jump To The Next Citation Point]) that discrete and L2 w-norms are uniformly equivalent, therefore:
∫ 1 (uN (x,t))2w (x) ≤ 2 max |u0(x )|2, (168 ) −1 x∈[0,1]
which proves the stability of the Chebyshev collocation method for the heat equation. Convergence can again be deduced from the Lax–Richtmyer theorem, but a detailed analysis cf. Section 6.5.1 of Canuto et al. [57Jump To The Next Citation Point]) shows that the numerical solution obtained by the method described here converges to the true solution and one can obtain the convergence rate. If the solution u(x, t) is m-times differentiable with respect to the spatial coordinate x (see Section 2.4.4) the energy norm of the error decays as 1−m N. In particular, if the solution is ∞ 𝒞, the error decays faster than any power of N. Advection equation
We now study the Legendre-tau approximation to the simple advection equation

∂u- ∂u- ∂t + ∂x = 0, with − 1 < x < 1, t > 0, (169 )
with homogeneous Dirichlet boundary condition
∀t ≥ 0, u (− 1,t) = 0, (170 )
and initial condition
∀ − 1 ≤ x ≤ 1, u(x,0) = u0(x ). (171 )
If we seek the solution as the truncated Legendre series:
∑N uN (x, t) = ai(t)Pi (x ) i=0
by the tau method, then u N satisfies the equation:
∂uN ∂uN -∂t--+ -∂x--= τN(t)PN (x). (172 )
Equating coefficients of PN on both sides of (172View Equation), we get
daN τN = ----. dt
Applying the L2w scalar product with uN to both sides of Equation (172View Equation), we obtain
∫ 1 ∂ ( 2 2 ) 1 ∂uN 1 2 2∂t- ∥uN ∥ − aN = − u-∂x-dx = − 2u N(1) ≤ 0, −1
which implies the following inequality:
d N∑ −1 --- a2i ≤ 0. (173 ) dt i=0
Finally, aN (t) is bounded because it is determined in terms of {ai}i=0...N− 1 from the boundary condition (170View Equation), and thus, stability is proved. In the same way as before, for the heat equation, it is possible to derive a bound for the error ∥u(x,t) − uN (x,t)∥, if the solution u(x,t) is m-times differentiable with respect to the spatial coordinate x; the energy norm of the error decays like N 1− m (see also Section 6.5.2 of Canuto et al. [57]). In particular, if the solution is ∞ 𝒞, the error decays faster than any power of N.
  Go to previous page Go up Go to next page