Consistency, stability and convergence for linear problems 26

In this section we consider more generally an initial boundary value problem (IBVP) defined on a domainΩby an equation of the form

∂tu(t, x) =L[u(x, t)], (3.7) with initial datau(0, x) = u0(x)and appropriate conditions on the boundary

∂Ωof the domain. The unknown solutionu(t, x)is a function of timetand of the space coordinatex. For simplicity we carry on the discussion in the one-dimensional case, but most of it can be extended to higher (typically three) dimensions. In the same spirit, the unknownuis treated as a scalar function, but the discussion can be extended to systems of equations, whereuis actually a vector field (in which case it is written asu).Lis a differential operator acting onu, and we consider in general operators that depend at most on the second derivative in space ofu.

We define a computational mesh in space, composed by the set of points {xi}(often taken to be equally spaced), as well as succession of time steps{tⁿ}. Letuandu˜be the true, analytic solution to (3.7) and an approximate solution obtained via some numerical method, respectively.uⁿ_i andu˜ⁿ_i are the values of these functions att = tⁿ andx= xi (the precise sense in which these values represent the values of the corresponding function at a given point depends on the choice of the numerical method and will made explicit in discussing each of them).

The numerical method used to discretize the operatorLis denoted byL˜^∆. It depends on some discretization parameter∆, typically related to the grid spacing,i.e.the distance between two adjacent pointsxiandxi+1of the mesh.

We consider initially the case in which the operatorLis linear, which greatly reduces the complexity of the problem at hand. We define thepointwise errorof a numerical scheme as

Eⁿ_i = ˜uⁿ_i −uⁿ_i , (3.8) i.e. the difference of the numerical solution to the true one at a given point in time and space. Theglobal erroris the norm||E(t)||of the pointwise error over the computational domain. The usual choice for the norm is theL¹-norm, or sometimes theL²-norm. The use of the∞-norm, while in principle desirable, leads to unrealistically stringent conditions for discontinuous solutions.

A numerical method is said to beconvergentif

∆lim→0||E(t)||= 0 ∀t , (3.9) i.e.if at all times the global error vanishes as the discretization parameter (grid spacing) tends to zero.

Thelocal truncation errorof the numerical scheme is defined as

H_iⁿ= ˜L^∆[u]− L[u]. (3.10) The local truncation error is the difference between the original operatorLand its discrete versionL˜^∆, both applied to the true solutionu, and thus measures

how well the discrete version of the equation (3.7) approximates the original one locally. We say that a numerical scheme isconsistentif

∆→0lim ||H(t)||= 0 ∀t , (3.11) i.e.if the local truncation error converges to zero in the continuum limit, for all possible initial datau0. In particular a scheme is of orderpif||H(t)||=O(∆^p).

Finally, a scheme is said to bestableif the norm of the local truncation error is limited:

sup

u6=0

||H(t)||

||u|| ≤C ∀t , (3.12)

where the constantCdoes not depend onu.

For linear equations, the Lax-Richtmeyer equivalence theorem guarantees that if a scheme is stable and consistent, then it is convergent (Lax and Richt-myer,1956;Richtmyer and Morton,1994). Furthermore, if the scheme is of orderp, thenE(t) =O(∆^p),i.e.the global error converges to zero with thep-th power of the discretization parameter.

3.1.3 Non-linear stability

To study the stability and convergence of non-linear equations it turns out to be important to consider the properties of both the spatial discretization and the time discretization. In the previous section we worked with semi-discrete schemes, where time was continuous. This is generally a valid approximation because the error term associated with the time discretization is often negligi-ble compared to the one arising from the space discretization. In this section we will consider instead fully discrete schemes,i.e. where both the time and space dependence of the solution are discretized.

We introduce therefore a family of evolution operators{T^s}, depending on the positive real parameters(i.e. the time). They form a semi-group,i.e. T^s◦ T^t=T^t+s, and are such that

u(t) =T^t(u0), (3.13)

i.e. the initial datau0is evolved to the timetby the operatorT^t(seeKružkov (1970) for details). The discrete version of the operators is denoted byT˜s^∆, which depends yet again on the single parameter∆since the time discretiza-tion is usually linked to the space one by a stability condidiscretiza-tion,i.e.the Courant-Friedrichs-Lewy condition (see the next section3.1.4).

With these definitions, the fully discrete analogue of equation (3.7) can be written as

u(t+ ∆t) =T∆t^∆[u(t)], (3.14) and we can translate the definitions of the previous section in the fully discrete case: the truncation error is

H(t) = ˜T∆t^∆[u(t)]− T^∆t[u(t)] ; (3.15) a scheme is consistent ifH tends to zero as∆tends to zero for some choice of a norm; and in particular if||H(t)||=O(∆^r)the scheme is said to be of order

r; finally, the scheme is linearly stable if sup

u6=0

||H(t)||

||u|| ≤C ∀t , (3.16)

for someCconstant.

The Lax-Wendroff theorem (Lax and Wendroff,1960) guarantees that if the solution to a non-linear conservation law obtained with a consistent and con-servative scheme,i.e.such that

Z

Ω

T˜s^∆[v]dx= Z

Ω

vdx (3.17)

for anyv ∈L¹(Ω), converges in theL¹-norm, then the solutionuit converges to is a weak solution of the equation. Therefore if a convergence condition can be found, the Lax-Wendroff theorem will guarantee that the solution found is a weak solution (but not necessarily an entropic one, this has to be proven by different means).

We introduce now the concept of total variation of a functionv(x),TV(v), as

TV(v) = sup

→0

1

Z

Ω|v(x)−v(x−)|dx . (3.18) Ifvis differentiable this is equivalent to

TV(v) = Z

Ω|v⁰(x)|dx , (3.19) and this last expression can also be used more generally for distributions if the derivative is interpreted as a distribution derivative.

A scheme is said to beTV-stable if for all initial datau0 with finite total variation, there exist two positive constantsCand∆0such that

TV( ˜Ts^∆[u0])≤0 ∀∆<∆0. (3.20) ATV-stable scheme, if consistent and conservative, is convergent and the Lax-Wendroff theorem applies.

Most TV-stable schemes are actually so-called total variation diminishing (TVD) schemes,i.e.such that

TV( ˜T∆t^∆[u])≤TV(u). (3.21) A way to ensure a scheme to be TVD is to require it to bemonotone:

ifu≥valmost everywhere, thenT˜t^∆[u]≥T˜t^∆[v]. (3.22) Crandall and Tartar(1980) andCrandall and Majda(1980) proved that mono-tone schemes are TVD and converge to weak (and entropic solutions). How-everHarten et al.(1976) showed that monotone schemes are at most first order accurate.

To achieve high order accuracy without sacrificing linear stability, non-monotone TVD methods in one dimension have been formulated; but Good-man and LeVeque (1985) proved that while such schemes exist in the one-dimensional case, in higher dimensions they are again limited to first order

accuracy. Therefore even weaker stability conditions have been considered, such as schemes satisfying as maximum principle or monotonicity preserving schemes,i.e. schemes that (at least in one-dimension) cannot generate spuri-ous extrema. In practice, many schemes commonly employed, including the ones presented here, have not been provenTV-stable or TVD, although nu-merical evidence seem to confirm that they converge to the correct entropic solutions. Furthermore for systems of equations no scheme has been proven stable or convergent for generic initial data, even in the one-dimensional case (Leveque,2002).