Two-level schemes for the advection equation

1 Introduction

Mathematical models of continuum mechanics Batchelor ; LandauLifshic1986

describe the transport of scalar and vector quantities due to advection. In particular, the basic equation of hydrodynamics is the continuity equation. Advective transfer causes the fulfillment of conservation laws

Godunov ; LeVeque . In addition, there are properties of the positivity and monotonicity of the solution. Such important properties of the differential problem must be inherited when passing to the discrete problem Anderson ; Wesseling .

For spatially approximation, conservative approximations are constructed on basis of using the conservative (divergent) form of the advection equation. Most naturally such technology is implemented when using the integro-interpolation method (balance method) on regular and irregular grids

Samarskii1989 , in the control volume method LeVeque ; Versteeg . The construction of monotonic approximations is discussed in many papers (see, for example, kulikovskii2000mathematical ; HundsdorferVerwer2003 ; Kuzmin ). In MortonKellogg1996 ; SamarskiiVabischevich1999a standard linear approximations are considered for convection-diffusion problems.

Currently, the main computing technology for solving applied problems is the finite element method Guermond ; Larson . It is widely used in computational fluid dynamics Donea ; Zienkiewicz . Monotonization of the solution is achieved by using various linear and nonlinear variants of stabilization techniques. It should be noted that the standard formulation of the equations of continuous medium mechanics in the conservative or non-conservative form is poorly suited for applying finite element approximations, for which the Hilbert spaces are natural.

Separate attention deserves the problems of constructing and investigating approximations in time. When solving boundary value problems for partial differential equations, two-level schemes (

θ

-method, schemes with weights) are traditionally widely used HundsdorferVerwer2003 ; Ascher2008 ; LeVeque2007 . Research of this schemes can be based on the general theory of stability (well-posedness) of operator-difference schemes Samarskii1989 ; SamarskiiMatusVabischevich2002

. In particular, unimprovable (coinciding necessary and sufficient) stability conditions can be used, which are formulated as operator inequalities in a finite-dimensional Hilbert space. The achieved level of theoretical stability research allows us to abandon the heuristic methods widely used in computational fluid dynamics to study the stability of difference schemes: the von Neumann method for stability analysis, Fourier analysis, the principle of frozen coefficients, the consideration of the problem without taking into account the boundary conditions.

In this paper, the standard finite-element approximation in space is used for the nonstationary equation. The advection equation is written in the so-called symmetric form SamarskiiVabischevich1999a , when the advection operator is the half-sum of the advection operators in the conservative (divergent) and non-conservative (characteristic) forms. Thus the continuum mechanics equations are written using the SD (Square root from Density) variables VabForm ; form1 . In this case the corresponding conservation laws are a direct consequence of the skew-symmetry of the advection operator. The conservativeness property is related to the preservation of the norm of the solution of the non-stationary advection equation, with stability with respect to the initial data. This property takes place not only for the solution, but also for some of its transformations. In this case, we are talking about the property of multiconservativeness. The principal point is related with the fact that the most important skew-symmetry property of the advection operator is inherited for finite element approximations.

For the advection equation an explicit two-level scheme as well as all schemes with a weight lower than 0.5, are aboslutly unstable. At the same time, extensive computational practice aims us to use the explicit schemes in these problems. In such schemes, the conditional stability (the Courant-Friedrichs-Lewy (CFL) condition) with a time-step limited by the Courant number is provided in fact by refusing the skew-symmetry of the advection operator — the use of dissipative approximations. In addition, the standard explicit schemes are not conservative.

We construct conditionally stable schemes for the advection equation based on the principle of regularization of operator-difference schemes Samarskii1967 ; Vabishchevich2014 , when stability is provided by a small perturbation. The explicit second-order Runge-Kutta scheme Butcher2008 ; HairerWanner2010 is considered. In this context, the classical explicit Lax-Wendroff scheme lax1964difference is also considered as regularization of the second-order Runge-Kutta scheme — the perturbation of the advection operator squared. Using the stability criteria of two-level operator-difference schemes, the stability conditions of regularized schemes and the explicit Lax-Wendroff scheme are obtained. Effective computational implementation of explicit schemes using finite element approximation in space is provided by using the diagonalization procedures of the mass matrix (mass-lumping procedure) Thomee2006 ; Cohen .

Of greatest interest are implicit two-level schemes for the advection equation, which belong to the unconditionally stable class. The classical Crank-Nicolson scheme has a second-order of accuracy, is unconditionally stable and multiconservative. For the advection problems under the consideration, we can use a scheme of the fourth-order of accuracy, which is also unconditionally stable and multiconservative. A certain drawback of this scheme is associated with the need to use the lumping procedure. An implicit version of the Lax-Wendroff scheme is proposed, which is conditionally stable, but has a higher accuracy than the explicit Lax-Wendroff scheme and does not require the diagonalization of the mass matrix.

The paper is organized as follows. A model two-dimensional problem for the advection equation is formulated in Section 2. Approximation in space is constructed using Lagrangian finite elements, the main properties of the problem solution are noted. In Section 3, we consider known and new explicit difference schemes for the advection equations, and investigate the stability conditions. Central for this work is Section 4. Implicit schemes, their stability and conservatism are studied here. In Section 5, numerical experiments on the accuracy of explicit and implicit schemes are discussed for the model IBV problem. The results of the work are summarized in Section 6.

2 Problem statement

In a bounded two-dimensional domain $Ω$ , we consider the advection equation written in the symmetric form. A standard finite element approximation in space is used. The problem of constructing approximations in time is formulated in such a way that the approximate solution inherits the basic properties of the solution of the differential problem.

2.1 Differential problem

The Cauchy problem is considered in the domain $Ω$ ( $x = (x_{1}, x_{2}) \in Ω$ )

\frac{d w}{d t} + A w = 0, 0 < t \leq T,

(1)

w (0) = w^{0},

(2)

using notation $w (t) = w (x, t)$ . The operator of advection (convective transport) $A$ is assumed to be constant and is written in the symmetric form:

A w = \frac{1}{2} d i v (v w) + \frac{1}{2} v \cdot g r a d w .

(3)

Thus, we take the half-sum of the transfer operator in the divergent (conservative part with $d i v (v w)$ ) and non-divergent (characteristic part with $v \cdot g r a d w$ ) forms SamarskiiVabischevich1999a ; zang1991rotation . The convective transport is determined by the velocity of the medium $v (x)$ , and the no-permeability condition on the boundary of the domain

(v \cdot ν) = 0, x \in \partial Ω,

(4)

where $ν$ is the outer normal to the boundary $\partial Ω$ .

The standard continuity equation has a divergent form:

\frac{\partial ϱ}{\partial t} + d i v (ϱ v) = 0,

(5)

where $ϱ$ is the density. From (4), (5), by direct integration, we obtain the law of mass conservation:

m (t) = c o n s t, m = \int_{Ω} ϱ (x, t) d x .

When orienting to finite elemental approximations in space, it is convenient to use of SD (Square root from Density) variables for the equations of hydrodynamics VabForm . Let $w = ϱ^{1 / 2}$ , then equation (5) is written in the form (1), (3).

In the Hilbert space $H = L_{2} (Ω)$ , we define the scalar product and norm in the standard way:

(w, u) = \int_{Ω} w (x) u (x) d x, ∥ w ∥ = (w, w)^{1 / 2} .

Under constraint (4), the convective transfer operator is skew-symmetric:

A = - A^{*},

(6)

and therefore, in particular,

(A u, u) = 0.

Using the scalar product of equation (1) and $w$ , taking into account (2), (6), we obtain

∥ w (t) ∥ = ∥ w^{0} ∥, 0 < t \leq T .

(7)

This relation with respect to the continuity equation (5) expresses the law of mass conservation. The property (7) for the solution of the the Cauchy problem (1), (2) is associated with the non-dissipativity of the system. The equation (1) for (6) has many other conservation laws. If the constant operator $B$ is permutable with $A$ , then we have the equality

∥ B w (t) ∥ = ∥ B w^{0} ∥, 0 < t \leq T .

(8)

In particular, equality (8) holds for $B = f (A)$ .

2.2 Finite element approximation in space

To solve numerically the problem (1)–(3), we employ finite element approximations in space (see, e.g., Thomee2006 ; brenner2008mathematical ). For (3), we define the bilinear form

a (w, u) = \int_{Ω} (\frac{1}{2} d i v (v w) u + \frac{1}{2} v \cdot g r a d w u) d x .

By (4), we have

a (w, u) = - a (u, w), a (w, w) = 0.

Define the subspace of finite elements $V^{h} \subset H^{1} (Ω)$ and the discrete operator $A$ as

(A w, u) = a (w, u), \forall w, u \in V^{h} .

The operator $A$ acts on the finite dimensional space $V^{h}$ and, similarly to (8), is skew-symmetric:

A = - A^{*} .

(9)

In finite element approximation, the main property of the advection operator, namely, its skew-symmetry, is inherited and, as a consequence, there is energy neutrality (the equality $(A w, w) = 0$ ).

The Cauchy problem (1), (2) is associated with the problem

\frac{d u}{d t} + A u = 0, 0 < t \leq T,

(10)

u (0) = u^{0},

(11)

for $u (t) \in V^{h}$ , where $u^{0} = P w^{0}$ with $P$ denoting $L_{2}$ -projection onto $V^{h}$ . Similarly (7), for the solution of the problem (10), (11), the conservative property is established:

∥ u (t) ∥ = ∥ u^{0} ∥, 0 < t \leq T .

(12)

The multiconservative property (see (8)) is associated with the equality

∥ B u (t) ∥ = ∥ B u^{0} ∥, 0 < t \leq T,

(13)

provided that

\frac{d}{d t} B = B \frac{d}{d t}, B A = A B .

When the problem (10), (11) is approximated in time, it is necessary to focus on the fulfillment of the properties (12) and (13) at separate time levels.

2.3 Approximation in time

Let, for simplicity, $τ$ be a step of a uniform grid in time such that $y^{n} = y (t^{n}), t^{n} = n τ$ , $n = 0, 1, . . ., N, N τ = T$ . To solve numerically the problem (10), (11), we use explicit and implicit two-level schemes, when the solution at the new level $y^{n + 1}$ is determined by the previously found solution $y^{n}$ , $n = 0, 1, . . ., N - 1$ .

When approximating in time, we focus, first of all, on the scheme stability. With respect to our problem, stability will be ensured by the following estimate of the solution at each time level:

∥ y^{n + 1} ∥ \leq exp (μ τ) ∥ y^{n} ∥, μ = c o n s t, n = 0, 1, . . ., N - 1.

(14)

The most favorable situation is associated with the fulfillment of the equality

∥ y^{n + 1} ∥ = ∥ y^{n} ∥, n = 0, 1, . . ., N - 1.

(15)

In this case, not only stability is ensured, but also the conservatism of the solution takes place (see (12)). Multiconservatism (see (13)) is due to the fulfillment of the equality

∥ B y^{n + 1} ∥ = ∥ B y^{n} ∥, n = 0, 1, . . ., N - 1,

(16)

for some operators $B$ .

Conservative approximations are related to the property of time reversibility for the problem (9)–(11). The computational algorithm is reversible in time if we calculate the solution at the level $t^{n + 1}$ and then change the transfer velocity ( $v (x)$ ) to the opposite (by $- v (x)$ ), then at the next step in time we get the solution that coincides with the solution at the level $t^{n}$ .

3 Explicit schemes

The explicit scheme for the problem (9)–(11) has the first-order error in time and is absolutely unstable. A conditionally stable scheme can be constructed on the basis of the regularization principle of difference schemes. We also consider a scheme of second-order approximation, namely, the Lax-Wendroff scheme.

3.1 Necessary and sufficient conditions for the stability of schemes with weights

Strict results on the stability of difference schemes in finite-dimensional Hilbert spaces are obtained in the theory of stability (well-posedness) of operator-difference schemes Samarskii1989 ; SamarskiiMatusVabischevich2002 . With respect to the subject of our investigation, we give the best possible conditions for the stability of a standard two-level scheme with weights ( $θ$ -method) to be unimprovable (matching necessary and sufficient) for the solution of the problem (10), (11).

When passing from the level $t^{n}$ to the level $t^{n + 1}$ , the approximate solution of the problem (10), (11) is determined from the equation

\frac{y^{n + 1} - y^{n}}{τ} + C (θ y^{n + 1} + (1 - θ) y^{n}) = 0, n = 0, 1, . . ., N - 1.

(17)

with some constant operator $C$ . In the simplest case, $C = A$ . The initial condition is

y^{0} = u^{0} .

(18)

If $θ = 0$ , we have then explicit scheme, for $θ = 1$ , we obtain the fully implicit scheme and for $θ = 0.5$ , the Crank-Nicolson scheme appears. The stability criterion is formulated as follows.

Lemma 1

Condition

(C y, y) + (θ - \frac{1}{2}) τ ∥ C y ∥^{2} \geq 0

(19)

is necessary and sufficient for the stability of the scheme (17), (18) in $H$ , and for the solution the level-wise inequality (14) holds with $μ = 0$ .

The proof of this statement is given, for example, in the book SamarskiiGulin1973 . The condition (19) can be written in the form of the operator inequality

C + (θ - \frac{1}{2}) τ C^{*} C \geq 0.

3.2 Explicit schemes of the first and second order of approximation

It is natural to start with the simplest explicit scheme for the problem (10), (11). In this case, $C = A, θ = 0$ in (17), i.e.

\frac{y^{n + 1} - y^{n}}{τ} + A y^{n} = 0, n = 0, 1, . . ., N - 1.

(20)

For the explicit scheme, the stability criterion (19) takes the form

(C y, y) - \frac{τ}{2} ∥ C y ∥^{2} \geq 0.

(21)

In the case of a skew-symmetric operator $C$ , inequality (21) is not satisfied for any $τ > 0$ . Therefore, the explicit scheme (18), (20) is absolutely unstable.

Among two-level schemes, explicit schemes of the second-order approximation are deserved separate consideration. Instead of (20) we will use the explicit second-order Runge-Kutta scheme:

\frac{y^{n + 1} - y^{n}}{τ} + A y^{n} - \frac{τ}{2} A^{2} y^{n} = 0, n = 0, 1, . . ., N - 1.

(22)

In this case, we have

C = A - \frac{τ}{2} A^{2} .

Taking into account (9), we get

\begin{matrix} C^{*} C & = (A^{*} + \frac{τ}{2} A^{*} A) (A + \frac{τ}{2} A^{*} A) = \frac{τ^{2}}{4} (A^{*} A)^{2} + A^{*} A . \end{matrix}

In view of this, the stability condition (21) takes the form

- \frac{τ^{3}}{8} ∥ A^{2} y ∥^{2} \geq 0.

Again, we cannot specify $τ > 0$ such that the stability of the scheme (22) holds.

Thus, we can formulate the following result.

Theorem 2

The explicit scheme of the first-order approximation (18), (20) and the explicit scheme of the second-order approximation (18), (22) are absolutly unstable in $H$ .

3.3 Conditionally stable non-standard scheme

We note the possibility of constructing a conditionally stable scheme based on the correction of the explicit scheme (18), (20). We consider more weak stability requirements, allowing the growth of the solution norm in accordance with (14) with the choice of some positive constant $μ > 0$ .

The construction of the considering approximation in time is based (see, for example, afanas2013unconditionally ; VabExact ) on the solution representation of the problem (10), (11) in the form

u (t) = exp (μ t) v (t) .

(23)

For $v (t)$ , from (10), (11), (23), we have the Cauchy problem

\frac{d v}{d t} + (A + μ I) v = 0,

(24)

v (0) = u^{0},

(25)

where $I$ is the identity operator.

For the approximate solution of the problem (24), (25), the explicit difference scheme is used

\frac{v^{n + 1} - v^{n}}{τ} + (A + μ I) v^{n} = 0, n = 0, 1, . . ., N - 1,

(26)

v^{0} = u^{0} .

(27)

Taking into account relation

y^{n} = exp (μ t^{n}) v^{n}, n = 0, 1, . . ., N

scheme (26), (27) corresponds to the use of

\frac{exp (- μ τ) u^{n + 1} - u^{n}}{τ} + (A + μ I) u^{n} = 0, n = 0, 1, . . ., N - 1,

(28)

with the initial condition (18). Such schemes belong to the class of non-standard mickens1994nonstandard ; mickens2002nonstandard . The level-wise estimate

∥ v^{n + 1} ∥ \leq ∥ v^{n} ∥

for the solution of the difference scheme (26), (27) corresponds to the estimate (14) for the scheme (18).

Theorem 3

The explicit scheme (18), (28) is stable in $H$ for

τ \leq \frac{2 μ}{μ^{2} + ∥ A ∥^{2}}, μ > 0.

(29)

In these conditions, for the problem solution the estimate (14) is satisfied.

[Proof.] It suffices to formulate the stability conditions of the scheme (26), (27). In this case $C = A + μ I$ and the stability criterion (21) gives

μ I - \frac{τ}{2} (μ^{2} I + A^{*} A) \geq 0.

This inequality will be satisfied with time-step constraints This inequality will be satisfied with constraints of time step (29).

For the advection problems under the consideration, an acceptable time step is associated with the Courant condition, when

τ \leq O (∥ A ∥^{- 1}) .

The time step limitation (29) is much more stiff: $τ \leq O (∥ A ∥^{- 2})$ . Therefore, the explicit scheme (18), (28) is not suitable for computational practice.

3.4 Regularized schemes

The regularization principle for difference schemes Samarskii1989 ; Vabishchevich2014 provides great opportunities in constructing difference schemes of a prescribed quality. The standard approach to the construction of stable schemes on the basis of the regularization principle is associated with the introduction of additional terms (regularizers) in operators of an origional (generating) difference scheme.

Instead of the explicit scheme (20), we will use the regularized scheme

\frac{y^{n + 1} - y^{n}}{τ} + A y^{n} + \frac{τ}{2} D y^{n} = 0, n = 0, 1, . . ., N - 1.

(30)

with the regularizer $D = D^{*} > 0$ . Let us formulate the stability conditions for this scheme.

For the scheme (30), we have

C = A + \frac{τ}{2} D .

Taking into account the skew-symmetry of the operator $A$ , we obtain

\begin{matrix} C^{*} C & = (A^{*} + \frac{τ}{2} D) (A + \frac{τ}{2} D) = \frac{τ^{2}}{4} D^{2} + A^{*} A . \end{matrix}

The stability condition (21) gives

\frac{τ^{2}}{4} D^{2} \leq D - A^{*} A .

(31)

Therefore, we can rely on the conditional stability of the scheme (30) for $D \geq A^{*} A$ .

Let $λ_{max}$

be the maximum eigenvalue of the spectral problem

D^{2} = λ (D - A^{*} A),

(32)

and $λ_{max} > 0$ . Then the inequality (31) will be satisfied for

τ \leq τ_{0}, τ_{0} = \frac{2}{λ_{max}^{1 / 2}} .

(33)

Thus, we can formulate the stability conditions.

Theorem 4

The explicit scheme (18), (30) for $D = D^{*} \geq A^{*} A$ is stable in $H$ if the time step satisfies the constraints (33), where $λ_{max}$ is the maximum eigenvalue of the spectral problem (32).

We highlight some interesting possibilities for the choice of the regularizer $D$ . The simplest version is related to the regularizer

D = β A^{*} A, β > 1.

(34)

The conditions (32), (33) give the following restrictions on the time step:

τ \leq 2 \frac{(β - 1)^{1 / 2}}{β} \frac{1}{∥ A ∥} .

(35)

Thus, the explicit scheme (18), (30), (34) is conditionally stable under the time step constraints of Courant type (see (35)). In the special case, for $β = 2$ , we have

τ \leq \frac{1}{∥ A ∥} .

The regularized scheme (18), (30), (34) has the first-order approximation in time. It is natural to expect a higher accuracy when choosing $D \approx A^{*} A$ , i.e. when the scheme (30) is close to the explicit second-order Runge-Kutta scheme (22). Here we can indicate two variants of such schemes.

The first variant is related with choosing the regularizer $D$ in the form

D = (1 + β τ) A^{*} A, β > 0.

(36)

The stability condition (21) for this case takes the form

\frac{τ^{2}}{4} (1 + β τ)^{2} (A^{*} A)^{2} \leq β τ A^{*} A .

It will be satisfied for

τ \leq 4 β \frac{1}{∥ A ∥^{2}} .

(37)

The restrictions (37) for the scheme (18), (30), (36) are substantially more strong than the restrictions (35) for the scheme (18), (30), (34). This disadvantage is compensated by the fact that we have the second-order approximation in time instead of the first order.

Let us now consider the second variant of the regularized scheme of the second-order approximation. The most well-known modification of the absolutely unstable explicit second-order Runge-Kutta scheme (18), (22) is the Lax-Wendroff scheme lax1964difference ; leveque2002finite . This scheme is actually based on replacing the operator $- A^{2} = A^{*} A$ in the explicit second-order Runge-Kutta scheme by a self-adjoint non-negative operator $Q$ close to it.

We define the new bilinear form

q (w, u) = \int_{Ω} (\frac{1}{2} d i v (v w) + \frac{1}{2} v \cdot g r a d w) (\frac{1}{2} d i v (v u) + \frac{1}{2} v \cdot g r a d u) d x .

By virtue of this definition, we have

q (w, u) = q (u, w), q (w, w) \geq 0.

Define the discrete operator $Q$ as

(Q w, u) = q (w, u), \forall w, u \in V^{h} .

(38)

The regularized scheme (18), (30) for $D = Q$ is a finite element version of the explicit Lax-Wendroff scheme. The stability conditions for $D \geq A^{*} A$ on a finite element space $V^{h}$ are formulated in the theorem 4.

3.5 Computational implementation of explicit schemes

The explicit schemes under the consideration are explicit only in form, according to the time approximation. The computational implementation of explicit finite-element approximations, unfortunately, is connected with the solution of systems of linear equations on a new level in time. The standard approach is related to the correction of approximations based on mass lumping (see, for example, Thomee2006 ).

The use of certain schemes for the Cauchy problem (10), (11

) is based on the solution of matrix problems for finding an approximate solution at a new time level. We associate with the differential-operator equation the corresponding system of ordinary differential equations.

We consider a standard quasi-uniform triangulation of the domain $Ω$ into triangles (or tetrahedra in 3-D). Let $x_{i}, i = 1, 2, . . ., N_{h}$ be vertices of this triangulation. We introduce the finite dimensional space $V^{h} \subset H^{1} (Ω)$ of continuous functions that are liner over each finite element, see, e.g. Thomee2006 . As a nodal basis we take the standard hat function $χ_{i} (x) \in V^{h}, i = 1, 2, . . ., N_{h}$ . Then for $v \in V^{h}$ , we have the representation

v (x) = N_{h} \sum i = i v_{i} χ_{i} (x) .

where $v_{i} = v (x_{i}), i = 1, 2, . . ., N_{h}$ .

From (10), we obtain the equation

M \frac{d z}{d t} + K z = 0,

(39)

where $z (t)$ is the vector of unknowns $z_{i} (t), i = 1, 2, . . ., N_{h}$ . Here $M = (m_{i j})$ is the mass matrix and $K = (k_{i j})$ is the stiffness matrix, and

m_{i j} = (χ_{i}, χ_{j}), k_{i j} = a (χ_{i}, χ_{j}), i, j = 1, 2, . . ., N_{h} .

If we set $u = M^{1 / 2} z$ , then equation (39) is written in the form (9), (10) for

A = M^{- 1 / 2} K M^{- 1 / 2}, M = M^{*} > 0, K = - K^{*} .

(40)

Thus, we can construct approximations in time for the equation (39) on the basis of approximations for equation (10), taking into account that $u = M^{1 / 2} z$ and (40).

For example, the explicit scheme (20) corresponds to the scheme

M \frac{z^{n + 1} - z^{n}}{τ} + K z^{n} = 0, n = 0, 1, . . ., N - 1.

(41)

It is necessary to solve the problem

M z^{n + 1} = M z^{n} - τ K z^{n}

at a new level in time. Because of this, the scheme (41) is explicit in terms of time approximation, but implicit in terms of the computational implementation, since the symmetric matrix $M$ is off-diagonal.

In order to provide an explicit computational implementation, various diagonalization procedures for the mass matrix are used:

M ⟶ ˜ M .

In the simplest mass lumping procedure Thomee2006 we have

˜ M = d i a g {{˜ m}_{1}, {˜ m}_{2}, . . ., {˜ m}_{N_{h}}}, {˜ m}_{i} = N_{h} \sum j = 1 m_{i j}, i = 1, 2, . . ., N_{h} .

Instead of (39), we seek an approximate solution of the Cauchy problem for the equation

˜ M \frac{d z}{d t} + K z = 0.

This equation corresponds to setting

A = {˜ M}^{- 1 / 2} K {˜ M}^{- 1 / 2}, ˜ M = {˜ M}^{*} > 0, K = - K^{*},

(42)

in the problem (9)–(11).

4 Implicit schemes

Unconditionally stable schemes for the advection equation are built on the basis of implicit approximations. It seems reasonable to employ the Crank-Nicolson scheme, which has the second order of accuracy in time and has conservative properties. In addition, schemes of higher accuracy are outlined, which have the fourth-order approximation. The implicit version of the Lax-Wendroff scheme is highlighted.

4.1 The Crank-Nicolson scheme

Among the implicit schemes for the problem (10), (11), the most important is the symmetric scheme whith $C = A, θ = 0.5$ in (17). In particular, it has the best SM (Spectral Mimetic) properties in the class of schemes with the skew-symmetric operator Vabischevich2011 . In this case, the solution is determined from the equation

\frac{y^{n + 1} - y^{n}}{τ} + A \frac{y^{n + 1} + y^{n}}{2} = 0, n = 0, 1, . . ., N - 1.

(43)

By virtue of (19), the scheme (18), (43) is absolutely stable and approximates the problems (10), (11) with the second order in $τ$ .

The Crank-Nicolson scheme is nondissipative. To show this fact, it is sufficient to multiply the scalar equation (43) by $τ (y^{n + 1} + y^{n})$ . Taking into account the skew-symmetry of the operator $A$ , we have

∥ y^{n + 1} ∥^{2} - ∥ y^{n} ∥^{2} = 0, n = 0, 1, . . ., N - 1.

Taking into account (18), we arrive at (15). Also there is the multiconservative property, see (13). For any constant operator $B$ , which commutes with $A$ , (16) is satisfied. We immediately see that the scheme (18), (43) generates a computational algorithm that is time reversible. The result of our consideration is the following theorem.

Theorem 5

The scheme (18), (43) provides a time-reversible computational algorithm and is unconditionally stable in $H$ . The scheme is conservative and multiconservative in the sense of the fulfillment of the equalities (15) and (16) for all operators $B$ that are constant and permutable with $A$ .

4.2 Scheme of higher accuracy order

The Crank-Nicolson scheme (43) corresponds to the use of the following Pade approximant for the exponential:

exp (- z) = \frac{1 - \frac{1}{2} z}{1 + \frac{1}{2} z} + O (z^{3}) .

When considering more accurate approximations, we can consider

exp (- z) = \frac{1 - \frac{1}{2} z + \frac{1}{12} z^{2}}{1 + \frac{1}{2} z + \frac{1}{12} z^{2}} + O (z^{5}) .

The corresponding difference scheme ( $z = τ A$ ) has the form

(I + \frac{1}{12} τ^{2} A^{2}) \frac{y^{n + 1} - y^{n}}{τ} + A \frac{y^{n + 1} + y^{n}}{2} = 0, n = 0, 1, . . ., N - 1.

(44)

We set

E = I + \frac{1}{12} τ^{2} A^{2},

and apply scalar multiplication of equation (44) by $τ (y^{n + 1} + y^{n})$ . This leads to the equality

(E y^{n + 1}, y^{n + 1}) = (E y^{n}, y^{n}), n = 0, 1, . . ., N - 1.

For $E > 0$ , this equality, taking into account $E A = A E$ , implies the stability estimate, the conservativeness and multiconservative properties. The positivity of the operator $E$ is satisfied for at not very large steps in time. Really,

E = I - \frac{1}{12} τ^{2} A^{*} A \geq I - \frac{1}{12} τ^{2} ∥ A ∥^{2} I .

Thus $E > 0$ for

τ < \frac{2 \cdot 3^{1 / 2}}{∥ A ∥} .

(45)

The restrictions (45) can be removed by modifying the proof. We write the scheme (44) in the form

S_{1} y^{n + 1} = S_{2} y^{n}, n = 0, 1, . . ., N - 1,

(46)

with notation

S_{1} = I + P + \frac{1}{3} P^{2}, S_{2} = I - P + \frac{1}{3} P^{2}, P = - P^{*},

(47)

and $P = \frac{τ}{2} A$ . The first question is related with the proof of the existence $S_{1}^{- 1}$ . Necessary and sufficient condition for the existence of the inverse operator $S_{1}^{- 1}$ is (see, for example, (Anderson, ; Lusternik, )) the fulfillment of the inequality

∥ S_{1} y ∥ \geq δ ∥ y ∥, δ > 0,

wherein $∥ S_{1}^{- 1} ∥ \leq δ^{- 1}$ .

Lemma 6

For the operator $S_{1}$ , which is defined according to (47), the inequality is satisfied

∥ S_{1} y ∥ \geq ∥ y ∥ .

(48)

[Proof.] The inequality (48) is equivalent to the inequality

(S_{1}^{*} S_{1} y, y) \geq (y, y) .

We have

\begin{matrix} S_{1}^{*} S_{1} & = (I + \frac{1}{3} P^{2} - P) (I + \frac{1}{3} P^{2} + P) = {(I + \frac{1}{3} P^{2})}^{2} - P^{2} = I - \frac{1}{3} P^{2} + \frac{1}{9} P^{4} = I + \frac{1}{3} P^{*} P + \frac{1}{9} P^{*} P P^{*} P \geq I \end{matrix}

taking into account the skew-symmetry of the operator $P$ .

Taking into account lemma 6, we can write (46) as follows

y^{n + 1} = S y^{n}, S = S_{1}^{- 1} S_{2}, n = 0, 1, . . ., N - 1.

(49)

For the transition operator from one level in time to another, the following statement holds.

Lemma 7

The operator $S = S_{1}^{- l} S_{2}$ under the conditions (47) is unitary, i.e.

S^{*} S = I .

(50)

[Proof.] We have

S = {(I + \frac{1}{3} P^{2} + P)}^{- 1} (I + \frac{1}{3} P^{2} - P),

S^{*} = (I + \frac{1}{3} P^{2} + P) {(I + \frac{1}{3} P^{2} - P)}^{- 1} .

Taking into account the permutability of the operator factors on the right-hand side, we get the equality (50).

The noted unitarity property (50) of the operator $S$ allows (49) to come to the conservativity property (15). Similarly, the multiconservative property is established. The result of our consideration is the following statement.

Theorem 8

The fourth-order accuracy scheme (18), (44) provides a time-reversible computational algorithm and is unconditionally stable in $H$ . The scheme (18), (44) is conservative and multiconservative in the sense of the equality (15) and (16) for all operators $B$ that are constant and permutable with $A$ .

In the computational implementation, we solve the Cauchy problem for equation (39). Because of this, the representation (40) is used in equation (10). The scheme (44) corresponds to the scheme

\begin{matrix} (M + \frac{1}{12} τ^{2} K M^{- 1} K) \frac{z^{n + 1} - z^{n}}{τ} + & K \frac{z^{n + 1} + z^{n}}{2} = 0, n = 0, 1, . . ., N - 1. \end{matrix}

(51)

Thus, for the transition to a new time level, it is necessary to solve equation $R z^{n + 1} = r^{n}$ with the operator

R = M + \frac{1}{2} τ K + \frac{1}{12} τ^{2} K M^{- 1} K .

This makes the scheme (18), (44) practically useless.

As in the case of explicit schemes, the diagonalizing procedure of the mass matrix saves the situation. In this case, instead of (40), we use the representation (42), and the scheme (45) is replaced by the scheme

\begin{matrix} (˜ M + \frac{1}{12} τ^{2} K {˜ M}^{- 1} K) \frac{z^{n + 1} - z^{n}}{τ} + & K \frac{z^{n + 1} + z^{n}}{} 2 = 0, n = 0, 1, . . ., N - 1. \end{matrix}

(52)

The scheme (52) is stable under constraint (45) taking into account the fact that the representation (42) holds.

4.3 Implicit Lax-Wendroff scheme

On the base of the implicit scheme (44), we can construct an implicit version of the Lax-Wendroff scheme. Taking into account the notation introduced above, instead of (44), we will use the scheme

(I - \frac{1}{12} τ^{2} Q) \frac{y^{n + 1} - y^{n}}{τ} + A \frac{y^{n + 1} + y^{n}}{2} = 0, n = 0, 1, . . ., N - 1.

(53)

As in the case of the explicit scheme, the construction is based on replacing the operator $- A^{2} = A^{*} A$ by the operator $Q$ , which is defined by (38).

We rewrite the scheme (53) in the form

(I + \frac{1}{12} τ^{2} A^{2}) \frac{y^{n + 1} - y^{n}}{τ} + A \frac{y^{n + 1} + y^{n}}{2} - \frac{1}{12} τ^{2} R \frac{y^{n + 1} - y^{n}}{τ} = 0,

where

R = Q - A^{*} A .

The operator $R$ explicitly highlights the approximation error in space due to the replacement of $- A^{2}$ by $Q$ .

For a positive constant selfadjoint operator $E$ , we define the Hilbert space $H_{E}$ , where the scalar product and norm are defined as follows

(u, v)_{E} = (E u, v), ∥ u ∥_{E} = (u, u)_{E}^{1 / 2} .

Analogously to theorem 8, the following result is formulated.

Theorem 9

The implicit Lax-Wendroff scheme (18), (53) provides a time-reversible computational algorithm and is stable for the constraint

τ < \frac{2 \cdot 3^{1 / 2}}{∥ Q ∥^{1 / 2}},

(54)

and for the solution, the following properties takes place

∥ y^{n + 1} ∥_{E} = ∥ y^{n} ∥_{E}, E = I - \frac{1}{12} τ^{2} Q, n = 0, 1, . . ., N - 1.

(55)

The condition (54) ensures the positivity of the operator $E$ . To prove the equality (55), it is sufficient to multiply equation (53) by $τ (y^{n + 1} + y^{n})$ . Thus, we have the conservative property not in $H$ (see (15)), but only in $H_{E}$ (see (55)). For the scheme (53), we cannot prove the multiconservative property.

5 Numerical experiments

The possibilities of the explicit and implicit schemes under the consideration for solving the advection equation are illustrated by the numerical results for the model problem.

5.1 Model problem

We will consider the problem (1)-(3) in the unit square:

Ω = {x | x = (x_{1}, x_{2}), 0 < x_{α} < 1, α = 1, 2} .

The initial condition is taken in the form

w^{0} (x) = 2 \cdot 10^{3} x_{1}^{2} (1 - x_{1})^{4} x_{2}^{2} (1 - x_{2})^{4} .

The components of the velocity $v = (v_{1}, v_{2})$ are defined by the stream function

ψ (x) = \frac{1}{π} sin (π x_{1}) sin (π x_{2}),

so that

v_{1} = \frac{\partial ψ}{\partial x_{2}}, v_{2} = - \frac{\partial ψ}{\partial x_{1}} .

The calculations are performed for $T = 5$ .

The computational code was implemented using the FEniCS fenics numerical framework. The finite element approximation in space is based on the use of continuous $P_{1}$ Lagrange element, namely, piecewise-linear elements. A uniform grid is used for spatial domain. The grid with the step $h = 0.02$ is shown in Fig. 1.

Figure 1: Grid with

h = 0.02

The accuracy of different approximations in time will be estimated by a reference solution. It was obtained using the scheme under the consideration with an essentially small time step: $τ = 2^{- 11} 10^{- 2}$ . The time-evolution of the solution on the grid with $h = 0.01$ (the main grid in space) is illustrated in Fig. 2. For the relative error of the approximate solution, we have

ε (t) = \frac{∥ y - ¯ y ∥}{∥ ¯ y ∥},

where $¯ y$ is the reference solution.

Figure 2:
The solution of the problem at different time-moments.

5.2 Explicit schemes

Prevously the explicit regularized schemes (18), (30) have been outlined. For the choice of a regularizer in the form (34), stability takes place with constraints (35) on the time step. It should be noted that we must orient on schemes with the operator (42).

We present the results of calculations for the regularized scheme with $β = 1$ . The constraints for the step are related with the norm of the operator $A$ . Taking into account its skew-symmetry property, we have

∥ A ∥ = max i | λ_{i} |,

where $λ_{i}, i = 1, 2, . . ., N_{h}$ are the eigenvalues of the operator $A$ :

A φ = λ φ .

Taking into account the representation (40), this spectral problem corresponds to the spectral problem

K ψ = λ M ψ .

Similarly, using the procedure of mass lumping with allowance for (42), we get the spectral problem

K ψ = λ ˜ M ψ .

To solve the spectral problems with symmetrical matrices, we use the SLEPc library (Scalable Library for Eigenvalue Problem Computations) hernandez2005slepc ). We apply the Krylov-Schur algorithm, a variation of the Arnoldi method, proposed by stewart2002krylov . The calculation results of the norm of the operator $A$ according to (40) and (42) on different grids are presented in the table 1. It should be noted that $∥ A ∥ = O (h^{- 1})$ and the norm of the operator $A$ decreases approximately by two times when using mass lumping, because the maximum of permissible time step (see (35)) increases approximately twice.

space grid ( $h$ )	$∥ A ∥$ for $M$	$∥ A ∥$ for $˜ M$
0.02	1.05288993e+02	5.59579462e+01
0.01	2.16001186e+02	1.14622718e+02
0.005	4.37491174e+02	2.31964151e+02

Table 1: The norm of the operator

A

The time-history of the error for various time steps for the regularized scheme (18), (30), (34) with (42) and $β = 1$ are shown in the Fig. 3. Convergence with first order in $τ$ is observed.

Figure 3: Error time-history of the regularized scheme.

Using the mass lumping procedure and taking into account (38), the explicit Lax-Wendroff scheme is written as

˜ M \frac{z^{n + 1} - z^{n}}{τ} + K z^{n} + \frac{τ}{2} G z^{n} = 0, n = 0, 1, . . ., N - 1.

(56)

It corresponds to the regularized scheme (30), whereh the operator $A$ is defined according to (42), and for the operator $D$ , we have

D = {˜ M}^{- 1 / 2} G {˜ M}^{- 1 / 2} .

(57)

The stability conditions of the scheme (30) are formulated in theorem 4. The inequality $D \geq A^{*} A$ will be satisfied if the constant $η \geq 1$ in the inequality $D \geq η A^{*} A$ . We consider the spectral problem

A^{*} A φ = λ D φ,

then $η = λ_{m a x}^{- 1}$ . Taking into account (42), (57), this spectral problem corresponds the spectral problem

K^{*} {˜ M}^{- 1} K ψ = λ G ψ .

(58)

The constraints on time step (33) is related to the spectral problem (32). In the case (42), (56), it is equivalent to the spectral problem

G {˜ M}^{- 1} G ψ = λ (G - K^{*} {˜ M}^{- 1} K) ψ .

(59)

The calculation results for the determination of $η$ and $τ_{0}$ from the numerical solution of the spectral problems (58), (59) are given in table 2. These data demonstrate the conditional stability of the Lax-Wendroff scheme, moreover $τ_{0} = O (h)$ .

space grid ( $h$ )	$η$	$τ_{0}$
0.02	1.00098795e+00	1.73477111e-02
0.01	1.00025320e+00	8.47323207e-03
0.005	1.00006414e+00	4.17705891e-03

Table 2: Parameters of the explicit Lax-Wendroff scheme

We write the explicit Lax-Wendroff scheme (56) in the form

˜ M \frac{z^{n + 1} - z^{n}}{τ} + K z^{n} - \frac{τ}{2} K^{*} {˜ M}^{- 1} K z^{n} + \frac{τ}{2} R z^{n} = 0, n = 0, 1, . . ., N - 1,

Here the term $R z^{n}$ with

R = G - K^{*} {˜ M}^{- 1} K

distinguishes this scheme from the explicit Runge-Kutta scheme of the second-order accuracy in time. In our problem, the convergence of the approximate solution with the second order in time is manifested (see Fig. 4), first of all, in the initial time interval, at which the smoothness of the solution is large enough and the influence of the term $R z^{n}$ is insignificant. The first order of accuracy due to the term $R z^{n}$ begins to appear at $t \approx 1$ .

Figure 4: Error time-history of the Lax-Wendroff scheme.

5.3 Implicit schemes

In the class of implicit schemes, the Crank-Nicolson scheme (18), (43) seems to be the basic one. Its accuracy in solving the model problem is illustrated in Fig. 5. In comparison with the explicit regularized scheme (see Fig. 3), a higher accuracy of the approximate solution is observed, convergence is approximately of the second order in $τ$ . In addition, this scheme stands out among all those considered schemes due to the fact that it is unconditionally stable.

Figure 5: Error time-history of the Crank-Nicolson scheme.

Similar results for the scheme of higher order of accuracy are shown in Fig. 6. The approximate solution using the lumping procedure is found from the equation (52). This scheme is absolutly stable and demonstrates the high accuracy of the approximate solution for substantially large time steps.

Figure 6: Error time-history of the scheme of higher accuracy.

The accuracy of the implicit version of the Lax-Wendroff scheme is illustrated in Fig. 7. As for the explicit scheme (see Fig. 4), there is a higher accuracy for the initial time, when the solution is smooth in space. Further, the effect of $Q - A^{*} A$ occures that results in the difference between the implicit Lax-Wendroff scheme (53) and the fourth-order accuracy scheme (44). The scheme is stable for $τ < τ_{0}$ , where $τ_{0}$ corresponds to the right side (54). Numerical results for $∥ Q ∥$ and $τ_{0}$ , which are obtained from the solution of the spectral problem for the operator $∥ Q ∥$ , are given in table 3 and show that $∥ Q ∥ = O (h^{- 2})$ and $τ_{0} = O (h)$ .

space grid ( $h$ )	$∥ Q ∥$	$τ_{0}$
0.02	3.22933843e+04	1.92767512e-02
0.01	1.33745164e+05	9.47221570e-03
0.005	5.44513748e+05	4.69446600e-03

Table 3: Parameters of the implicit Lax-Wendroff scheme

Figure 7: Accuracy of the implicit Lax-Wendroff scheme.

6 Conclusions

In numerical simulation of advection processes, it is necessary to focus on writing the equation in a symmetric form (the half-sum of the advection operator in the divergent (conservative) form and the advection operator in the non-divergent (characteristic) form). In this case the advection operator is skew-symmetric. For the solution of the non-stationary advection equation, the multiconservative property holds, which is associated with the fulfillment of the set of conservation laws.
The stability of known and new two-level difference schemes is investigated for the approximate solution of the Cauchy problem for the advection equation. Investigation of stability is carried out on the basis of the general theory of stability (well-posedness) of operator-difference schemes. The standard finite-element approximation is used in space, which ensures the skew-symmetry of the discrete advection operator.
The standard explicit scheme for the advection equation belongs to the class of absolutely unstable ones. The explicit regularized scheme of the first order of accuracy is proposed, and the stability conditions are formulated. The stability conditions are established for the finite-element version of the classical explicit Lax-Wendroff scheme. For finite-element approximation in space, it is necessary to use various diagonalization algorithms of the mass matrix.
For solving advection problems, the Crank-Nicolson scheme seems to be very attractive. It demonstrates the second order of accuracy and belongs to the class of unconditionally stable schemes. In addition, it demonstrates many conservation laws. The possibilities of using schemes with higher order of accuracy, which is also unconditionally stable and multiconservative, are also considered. The implicit version of the Lax-Wendroff scheme is proposed, and conditions for its stability are formulated.
The properties of the explicit schemes are illustrated by the solution of the model two-dimensional problem for the advection equation. The focus is on studing the accuracy of various approximations in time.

Acknowledgements

This work is supported by the mega-grant of the Russian Federation Government (# 14.Y26.31.0013).

References

(1) G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 2000.
(2) L. D. Landau, E. Lifshitz, Fluid Mechanics, Butterworth-Heinemann, 1987.
(3) S. K. Godunov, E. I. Romenskii, Elements of Continuum Mechanics and Conservation Laws, Springer, 2003.
(4) R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
(5) J. Anderson, Computational Fluid Dynamics: The Basics with Applications, McGraw-Hill, 1995.
(6) P. Wesseling, Principles of Computational Fluid Dynamics, Springer, 2010.
(7) A. A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001.
(8) H. Versteeg, W. Malalasekra, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Prentice Hall, 2007.
(9) A. G. Kulikovskii, N. V. Pogorelov, A. Y. Semenov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Taylor & Francis, 2000.
(10) W. H. Hundsdorfer, J. G. Verwer, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, Springer Verlag, 2003.
(11) D. Kuzmin, A Guide to Numerical Methods for Transport Equations, University Erlangen-Nuremberg, 2010.
(12) K. W. Morton, R. B. Kellogg, Numerical Solution of Convection-Diffusion Problems, Chapman & Hall London, 1996.
(13) A. A. Samarskii, P. N. Vabishchevich, Numerical Methods for Solving Convection-Diffusion, URSS, Moscow, 1999.
(14) A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004.
(15) M. G. Larson, F. Bengzon, The Finite Element Method: Theory, Implementation and Applications, Springer, 2013.
(16) J. Donea, A. Huerta, Finite Element Methods for Flow Problems, Wiley, 2003.
(17) O. C. Zienkiewicz, R. L. Taylor, N. P., The Finite Element Method for Fluid Dynamics, Butterworth-Heinemann, 2013.
(18) U. M. Ascher, Numerical Methods for Evolutionary Differential Equations, Society for Industrial Mathematics, 2008.
(19) R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, Society for Industrial Mathematics, 2007.
(20) A. A. Samarskii, P. P. Matus, P. N. Vabishchevich, Difference Schemes with Operator Factors, Kluwer Academic, Dordrecht, 2002.
(21) P. N. Vabishchevich, On the form of the hydrodynamics equations, in: West-East High Speed Flow Field conference, 19-22, November 2007, Moscow, Russia, Moscow, 2007, pp. 1–9.
(22) A. Churbanov, P. Vabishchevich, Numetical methods for solving convection-diffusion problems, in: C. Zhao (Ed.), Focus on Porous Media Research, Nova Science Publishers, 2013, pp. 1–83.
(23) A. A. Samarskii, Regularization of difference schemes, Zh. Vychisl. Mat. Mat. Fiz. 7 (1967) 62–93, in Russian.
(24) P. N. Vabishchevich, Additive Operator-Difference Schemes: Splitting Schemes, de Gruyter, Berlin, 2014.
(25) J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, 2008.
(26) E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer Verlag, 2010.
(27) P. D. Lax, B. Wendroff, Difference schemes for hyperbolic equations with high order of accuracy, Communications on Pure and Applied Mathematics 17 (3) (1964) 381–398.
(28) V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Verlag, Berlin, 2006.
(29) G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Springer, 2002.
(30) T. A. Zang, On the rotation and skew-symmetric forms for incompressible flow simulations, Applied Numerical Mathematics 7 (1) (1991) 27–40.
(31) S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, New York, 2008.
(32) A. A. Samarskii, A. V. Gulin, Stability of Difference Schemes, Nauka, Moscow, 1973, in Russian.
(33) N. M. Afanas’eva, P. N. Vabishchevich, M. V. Vasil’eva, Unconditionally stable schemes for convection-diffusion problems, Russian Mathematics 57 (3) (2013) 1–11.
(34) P. N. Vabishchevich, Fundamental mode exact schemes for nonstationary problems, arXiv (arxiv.org/abs/1705.07010) (2017) 1–17.
(35) R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1994.
(36) R. E. Mickens, Nonstandard finite difference schemes for differential equations, The Journal of Difference Equations and Applications 8 (9) (2002) 823–847.
(37) R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002.
(38) P. N. Vabishchevich, Two-level schemes of higher approximation order for time-dependent problems with skew-symmetric operators, Computational Mathematics and Mathematical Physics 51 (6) (2011) 1050–1060.
(39) L. Lusternik, V. Sobolev, Elements of Functional Analysis, Wiley, 1975.
(40) A. Logg, K. A. Mardal, G. Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Springer, 2012.
(41) V. Hernandez, J. E. Roman, V. Vidal, SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Transactions on Mathematical Software (TOMS) 31 (3) (2005) 351–362.
(42) G. W. Stewart, A Krylov–Schur algorithm for large eigenproblems, SIAM Journal on Matrix Analysis and Applications 23 (3) (2001) 601–614.

Two-level schemes for the advection equation

Splitting Schemes for Some Second-Order Evolution Equations

A compact fourth-order implicit-explicit Runge-Kutta type scheme for numerical solution of the Kuramoto-Sivashinsky equation

Modern theory of hydraulic fracture modeling with using explicit and implicit schemes

Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials

Two-component domain decomposition scheme with overlapping subdomains for parabolic equations

Parallel Type Decomposition Scheme for Quasi-Linear Abstract Hyperbolic Equation

A stable semi-implicit algorithm

1 Introduction

2 Problem statement

2.1 Differential problem

2.2 Finite element approximation in space

2.3 Approximation in time

3 Explicit schemes

3.1 Necessary and sufficient conditions for the stability of schemes with weights

Lemma 1

3.2 Explicit schemes of the first and second order of approximation

Theorem 2

3.3 Conditionally stable non-standard scheme

Theorem 3

3.4 Regularized schemes

Theorem 4

3.5 Computational implementation of explicit schemes

4 Implicit schemes

4.1 The Crank-Nicolson scheme

Theorem 5

4.2 Scheme of higher accuracy order

Lemma 6

Lemma 7

Theorem 8

4.3 Implicit Lax-Wendroff scheme

Theorem 9

5 Numerical experiments

5.1 Model problem

5.2 Explicit schemes

5.3 Implicit schemes

6 Conclusions

Acknowledgements

References

Two-level schemes for the advection equation

Related Research

Splitting Schemes for Some Second-Order Evolution Equations

A compact fourth-order implicit-explicit Runge-Kutta type scheme for numerical solution of the Kuramoto-Sivashinsky equation

Modern theory of hydraulic fracture modeling with using explicit and implicit schemes

Stable explicit schemes for simulation of nonlinear moisture transfer in porous materials

Two-component domain decomposition scheme with overlapping subdomains for parabolic equations

Parallel Type Decomposition Scheme for Quasi-Linear Abstract Hyperbolic Equation

A stable semi-implicit algorithm

1 Introduction

2 Problem statement

2.1 Differential problem

2.2 Finite element approximation in space

2.3 Approximation in time

3 Explicit schemes

3.1 Necessary and sufficient conditions for the stability of schemes with weights

Lemma 1

3.2 Explicit schemes of the first and second order of approximation

Theorem 2

3.3 Conditionally stable non-standard scheme

Theorem 3

3.4 Regularized schemes

Theorem 4

3.5 Computational implementation of explicit schemes

4 Implicit schemes

4.1 The Crank-Nicolson scheme

Theorem 5

4.2 Scheme of higher accuracy order

Lemma 6

Lemma 7

Theorem 8

4.3 Implicit Lax-Wendroff scheme

Theorem 9

5 Numerical experiments

5.1 Model problem

5.2 Explicit schemes

5.3 Implicit schemes

6 Conclusions

Acknowledgements

References