Discretization in space: stability and convergence

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 (116

) 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)$

Two more concepts are helpful in the convergence analysis of numerical schemes:

consistency: an approximation to the PDE (116 ) is consistent if $∀v ∈ ℋ$ both
stability: with the formal notations of Equation (117 ), an approximation to the PDE (116 ) is stable if

where $C (t)$ is independent of $N$ and bounded for $t ∈ [0,T ]$ .

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 (116 ) is approximated by

To show that stability implies convergence, we subtract it from the exact solution (116

)

∂-(u-−-uN-) ∂t = LN (u − uN ) + Lu − LN u,

and obtain, after integration, (the dependence on the space coordinate $x$ is skipped)

Using the stability property (143

), the norm ( $2 Lw$ ) of this equation implies

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

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 (60 ). 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:

If the operator $L$ is semibounded:

where $I$ is the identity operator.
In the parabolic case, if $L$ satisfies the coercivity condition (see also Chapter 6.5 of Canuto et al. [57 ]¹ ):

and the continuity condition:
In the hyperbolic case, if there exists a constant $C > 0$ such that

and if the operator verifies the negativity condition:

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

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”

grows at most exponentially with time. Since $uN (t) = eLNtuN (0)$ for any $uN (0)$ , we obtain

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

4.3.3.1 Heat equation

We first study the linear heat equation

with homogeneous Dirichlet boundary conditions

and initial condition

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

which are time ODEs (discussed in Section 4.1) with the initial conditions

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 (159 ) 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:

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

and integrating by parts twice, one gets the relation

By the properties of the Chebyshev weight

it is possible to show that

and thus that

Therefore, integrating relation (161

) over the time interval $[0,t]$ , one obtains

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. [57

]) that discrete and $L2 w$ -norms are uniformly equivalent, therefore:

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. [57

]) 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$ .

4.3.3.2 Advection equation

We now study the Legendre-tau approximation to the simple advection equation

with homogeneous Dirichlet boundary condition

and initial condition

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:

Equating coefficients of $PN$ on both sides of (172

), we get

daN τN = ----. dt

Applying the $L2w$ scalar product with $uN$ to both sides of Equation (172

), 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:

Finally, $aN (t)$ is bounded because it is determined in terms of ${ai}i=0...N− 1$ from the boundary condition (170

), 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$ .