- Email: [email protected]
Some splitting methods for hyperbolic PDEs
Applied Numerical Mathematics 146 (2019) 361–378
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Some splitting methods for hyperbolic PDEs Rasool Hosseini a , Mehdi Tatari a,b,∗ a b
Department of Mathematical Sciences, Isfahan University of Technology, Isfahan, 84156-83111, Iran Institute for Research in Fundamental Sciences (IPM), P.O. Box: 19395-5746, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 10 January 2019 Received in revised form 18 May 2019 Accepted 4 July 2019 Available online 9 July 2019 Keywords: Advection equation Burgers’ equation Maximum principle Damping of solution
a b s t r a c t In this work, a new splitting technique is implemented for solving hyperbolic PDEs. As the main result, the new methods preserve the maximum principle unconditionally or with a mild condition on discretization parameters in comparison with well known methods. Damping of numerical solution in time evolution is investigated. For numerical solution of the Burgers’ equation, as a nonlinear problem, an iterative method based on the new splitting technique is presented. Efficiency of the new methods is examined via numerical examples. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Consider the advection equation
ut + au x = 0. This equation describes the passive advection of some scalar field u (x, t ) carried along by a flow of constant speed a. This equation must surely be the simplest of all partial differential equations. Yet to approximate it well on a fixed (x, t )-mesh is a far from trivial problem that is still under active discussion in the numerical analysis literature [23]. The exact solution of this equation is constant along the characteristic lines associated to this equation. Since finding characteristic lines are complicated especially for system of hyperbolic equations, usually use of numerical methods is preferred. Some well-known numerical methods for solving this equation are upwind and Lax-Wendroff methods. In numerical methods for solving advection equation, preservation of maximum principle is an important subject. The upwind method preserves this property under the CFL condition while Lax-wendroff always does not preserve the property. As an other important issue in this kind of equations, damping of the solution should be studied. In this paper some new methods are presented and their behaviour for solving advection equation are discussed. Burgers’ equation or Bateman-Burgers’ equation is a nonlinear hyperbolic partial differential equation occurring in various areas of applied mathematics, such as fluid mechanics, dispersive water and traffic flow. Also the Burgers’ equation appears often as a simplification of a more complicated model like Navier Stokes equations. The equation was first introduced by Harry Bateman [1] and later studied by Johannes Martinus Burgers’ [4]. For a given field u (x, t ) and kinematic viscosity , the general form of Burgers’ equation (also known as viscous Burgers’ equation or advective form of the Burgers’ equation) in one space dimension is the dissipative system:
*
Corresponding author. E-mail addresses: [email protected] (R. Hosseini), [email protected] (M. Tatari).
https://doi.org/10.1016/j.apnum.2019.07.005 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.
362
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
ut + uu x = u xx . When
= 0, Burgers’ equation becomes the inviscid Burgers’ equation
ut + uu x = 0, which is a prototype for conservation equations that can develop discontinuities as shock waves. The inviscid Burgers’ equation is a conservation equation, more generally a first order nonlinear hyperbolic equation. The solution of this equation with the initial condition u (x, 0) = f (x) can be constructed by many of numerical methods like the method of characteristics [19,28]. In recent years, computing the solution of Burgers’ equation has attracted a lot of attention and various methods have been introduced by researchers to solve this equation. Although the exact solution of this equation with some simple initial conditions is presented by Hopf [18], Benton and Platzman [2] and Cole [8], but in most situations, this equation can only be solved numerically. To see a complete overview of numerical methods for solving Burgers’ equation see [19]. Splitting methods introduced by Trotter [27] and Strang [25] are one of the most widely used methods among geometric integrators, which are often used to overcome the complexity of computing problems, constructing very accurate high order numerical algorithms and extension of the region of stability of methods. These ideas are used to overcome the computational complexity which is arisen in numerical solution of higher dimensional problems [10,12,24]. In [11] a second order accurate operator splitting method with unconditional stability is presented to solve a system of nonlinear hyperbolic PDEs that is obtained from the simulation of flows of isothermal compressible natural gas over transmission pipelines. The idea of the splitting for multidimensional linear hyperbolic systems and for Burgers’ equation as a nonlinear hyperbolic equation is introduced in [16]. In [20] a second-order operator splitting Fourier spectral method is presented for solving fractional in space reaction-diffusion equations. Splitting approaches for Schrödinger equation are given in [6,9]. Also, iterative splitting methods are one of the suitable methods for solving nonlinear equations [14,15]. The idea of splitting can be applied to wide range of problems [3,17]. For more information about splitting methods see [5,22,29]. In [7], it is shown that Saul’yev-type schemes can be derived from the exponential splitting of the semi-discretized equation and the splitting idea in which it is used is to rewrite the resulting matrix from the semi-discretization as a sum of matrices that have only a single, nonvanishing 2 × 2 matrix along the diagonal. Some splitting methods are introduced in [21] and recent developments are presented. New additive and iterative splitting schemes and algorithms are constructed in [13]. In this work some splitting methods are outlined for hyperbolic equations. Then a new idea which is different from existing schemes is presented. It is based on splitting of semi-discretized equations along the diagonal. The reminder of this paper is organized as follows. In Section 2, we present preliminaries about splitting methods. In Section 3, some proper splitting methods are implemented for solving advection and Burgers’ equations. Numerical results are investigated in Section 4 and finally a conclusion is drawn in Section 5. 2. Splitting methods In this section, basic concepts of some splitting methods that are used in the next section are introduced [3,17]. 2.1. Preliminary Consider the differential equation
dy(t ) dt
= A (y(t )) + B (y(t )),
y(0) = y0 ,
(2.1)
where y is a vector or possibly a matrix and A and B are operators on y. In a splitting method for solving (2.1) numerically, it is needed to solve exactly (or numerically) the following two subproblems
dy(t ) dt
= A (y(t )),
dy(t ) dt
= B (y(t )).
(2.2)
ϕtA and ϕtB the exact t flows for each subproblems in (2.2). This means that, for example, the solution of dy(t ) subproblem dt = A (y(t )) in time t is obtained by y (t ) = ϕtA ( y (t 0 )). In this way, we have the first order Lie-Trotter splitting Denote by
methods that one is the adjoint of the other
h∗ = ϕhB ◦ ϕhA , h = ϕhA ◦ ϕhB , where ◦ (composition operator) is an operation that takes two functions f and g and produces a function h such that h(x) = g ( f (x)). By composing the Lie-Trotter splitting method and its adjoint with a half time step, the Strang splitting formulae
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
363
h[s] = ϕhA/2 ◦ ϕhB ◦ ϕhA/2 ,
h[s] = ϕhB/2 ◦ ϕhA ◦ ϕhB/2 , will be obtained that are second order time symmetric methods. If Problem (2.1) can be decomposed into more than two parts as
dy(t ) dt
= A [1] (y(t )) + · · · + A [m] (y(t )),
and the exact flows
ϕh[i] correspond to each split part
dy(t ) dt
= A [i ] (y(t )), i = 1, ..., m be available, then the scheme
χh = ϕh[m] ◦ · · · ◦ ϕh[2] ◦ ϕh[1] , and its adjoint map
χh∗ = χ−−h1 = ϕh[1] ◦ ϕh[2] ◦ · · · ◦ ϕh[m] , are first order approximations (Lie-Trotter type) and the compositions
χh/2 ◦ χh∗/2
χh∗/2 ◦ χh/2 ,
or
provides second order methods (Strang type). 2.2. New splitting methods Here we introduce a new splitting technique for system of ordinary differential equations. For this, first consider the following preposition which will be used later. Preposition 2.1. Any strictly triangular matrix A ∈ Rn×n (a triangular matrix with zero on main diagonal) is nilpotent, in other words there is an integer m ≤ n such that Am = 0. Proof. By induction on n.
2
Now consider a linear system of ordinary differential equations
dx(t ) dt
= Ax(t ),
x(t 0 ) = x0 ,
(2.3)
with coefficients matrix A = (ai j ) ∈ Rn×n . By splitting matrix A into its 2n − 1 diagonals
A = L1 + · · · + Ln−1 + D + U1 + · · · + Un−1 , System (2.3) can be split as follows
dx(t ) dt
= L1 x(t ) + · · · + Ln−1 x(t ) + Dx(t ) + U1 x(t ) + · · · + Un−1 x(t ).
The exact flow of the system
dx(t ) dt
= Dx(t ),
is
x(t ) = e (Dt ) x0 = D(t )x0 ,
D(t ) = diag(e (a11 t ) , . . . , e (ann t ) ),
where D is a diagonal matrix that include main diagonal of matrix A. Let
Ek (t ) = e (Uk t ) =
∞ tj (U K ) j , k = 1, . . . , n − 1, j! j =1
∞ tj Sk (t ) = e (Lk t ) = (L K ) j , k = 1, . . . , n − 1. j! j =1
By Preposition 2.1 matrices Uk and Lk are nilpotent and Ek (t ) and Sk (t ) are calculable finitely, such that Ek (t ) is an upper triangular matrix with 1 on the main diagonal and
364
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
j−k i
s=1 ai +(s−1)k,i +sk
(Ek )i j =
t
j −i k
for
i ( j− )! k
j > i & j − i = mk, m = 1, . . . ,
n−1
k
,
(2.4)
and Sk (t ) is a lower triangular matrix with 1 on the main diagonal and
i− j k
(Sk )i j =
s=1 a j +sk, j +(s−1)k t
i− j k
for
( i −k j )!
j < i & i − j = mk, m = 1, . . . ,
n−1 k
.
(2.5)
dx(t ) = Uk x(t ) and = Lk x(t ) are x(t ) = Ek (t )x0 and x(t ) = Sk (t )x0 respecdt dt tively for k = 1, . . . , n − 1. Now let L(t ) = S1 (t ) · · · Sn−1 (t ) and U(t ) = E1 (t ) · · · En−1 (t ), then the Lie-Trotter splitting method for the system (2.3) is Therefore, the exact flows of the systems
dx(t )
x(t ) = L(t ) D(t )U(t )x0 ,
(2.6)
where the adjoint map of this method is
x(t ) = U(t ) D(t )L(t )x0 .
(2.7)
By composing (2.6) and (2.7) with half time step, we have the Strang splitting formula in the form
x(t ) =L(
=L(
t t t t t t )D( )U( )U( ) D( )L( )x0 2 t 2
) D(
2 t 2
2
)U(t ) D(
2 t 2
)L(
2 t 2
(2.8)
2
)x0 .
(2.9)
This idea specially is useful for solving systems that A is a band matrix. These systems appear in spatial discretization of partial differential equations. 3. New splitting methods and hyperbolic PDEs In this section, pure advection equation,
ut + au x = 0,
(3.1)
is considered, where u is the advection velocity and nonzero a ∈ R is advection coefficient. 3.1. An upwind type splitting method Consider nodes x j in [c , d] with x j = c + j x, j = 1, ..., J where x =
d−c J
, J ∈ N . The approximate solution is then
denoted by
U j ≈ u (x j , t ). As an upwind type method, we use a backward difference in space if a is positive and a forward difference if a is negative [23]
⎧ −a ⎪ U j +1 − U j if a < 0 ⎨ dU j x = . ⎪ dt ⎩ −a U − U if a > 0 j j −1 x Set
ν =a
t . Consider scheme (3.2) as x
dU dt
(3.2)
= AU,
(3.3)
where
A( J −1)×( J −1) = upperbidiag( or
a −a , ), x x
if a < 0,
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
A( J −1)×( J −1) = lowerbidiag(
a −a , ), x x
if a > 0.
Now assume a > 0, write A = D + N, where D = diag( follows
dU dt
365
−a a ) and N = lowerdiag( ). Therefore, System (3.3) can be split as x x
= DU + NU.
By using the splitting procedure explained in previous section, we have the upwind type splitting method for the Problem (3.1) as
U(t ) = exp(−ν )EU0 ,
(3.4)
where E is a lower triangular matrix with elements
⎧ i− j ⎪ ⎨ (ν ) e i j = (i − j )! ⎪ ⎩ 0
j≤i
(3.5)
.
otherwise
This method preserves the maximum principle property, more precisely we have
E∞ = 1 + ν + ≤
∞ (ν )i i =0
i!
(ν )2 (ν ) J − 2 + ··· + 2! ( J − 2)! = exp(ν ),
and if we set G = exp(−ν )E, then
Gn ∞ ≤ exp(−nν )En∞ , and this concludes that
U(nt )∞ ≤ U(0)∞ . Damping of the solution is an important feature of numerical methods for solving advection equation. More precisely, since there is no damping in the analytical solution, we would like to have G∞ = 1. In the next theorem, damping of the new method is investigated. Theorem 3.1. Maximum damping in time step n + 1 of upwind type splitting method (3.4) is bounded by
−
ν J −1 ( J − 1)!
U(nt )∞ .
Proof. We know
G∞ = exp(−ν )E∞ = exp(−ν )(exp(ν ) − =1−
exp(ξ )ν J −1
exp(ξ − ν )ν J −1
( J − 1)!
( J − 1)!
)
,
where 0 ≤ ξ ≤ ν . Therefore, maximum damping of amplitude matrix is
−
ν J −1 ( J − 1)!
.
Thus, maximum damping in time step n + 1 is
−
ν J −1 ( J − 1)!
U(nt )∞ . 2
366
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
If a < 0, the proposed splitting method will be as
U(t ) = exp(ν )SU0 ,
(3.6)
where S is an upper triangular matrix with
⎧ (−ν ) j −i ⎪ ⎪ ⎨ ( j − i )! si j = ⎪ ⎪ ⎩ 0
j≥i (3.7)
. otherwise
Here, results are similar to the case a > 0. Remark 3.2. For a > 0, notice that elements on each diagonal of G are a multiplier of elements on its main diagonal. More precisely, element on main diagonal of G are exp(−ν ), and elements on the next diagonals are etc. Similar strategy should be considered when a < 0.
ν exp(−ν ),
ν 2 exp(−ν ) 2!
and
3.2. A second order time splitting method As we discussed above, the upwind type method resulting from the splitting scheme is first order accurate. One simple way to achieve second order accuracy on the advection equation is to use an artificial diffusion term u xx with coefficient > 0. Applying this directly to the (3.1) and use a second order discretization on space, gives the second order scheme
dU j dt
=
−a U j +1 − U j −1 + U j +1 − 2U j + U j −1 . 2 x
(3.8)
Consider this scheme as
dU dt
= AU,
(3.9)
where
A( J −1)×( J −1) = tridiag(
x2
+
a
,
−2
,
2x x2 x2
−
a 2 x
).
a a 2 Now Write A = M + D + N, where M = upperdiag( 2 + 2 ), D = diag( − ) and N = lowerdiag( x2 − ). Therefore, x x x2 2x System (3.9) can be split as follows
dU dt
= MU + DU + NU.
Let μ = t2 and ν = axt . By using the explained splitting procedure and introducing F(t ) = e (−μ) S and F(t ) = e (−μ) H, x the Lie-Trotter type of proposed splitting method for the Problem (3.8) is
U(t ) = F(t ) F(t )U(0),
(3.10)
and the Strang type formula is
U(t ) = F(
t t )F(t )F( )U(0), 2
(3.11)
2
where H is an upper triangular matrix with elements
⎧ 1 ( j −i ) ⎪ ⎪ ⎨ (μ − 2 ν ) hi j = ( j − i )! ⎪ ⎪ ⎩ 0
j≥i
,
otherwise
and S is a lower triangular matrix with elements
⎧ 1 (i − j ) ⎪ ⎪ ⎨ (μ + 2 ν ) si j = (i − j )! ⎪ ⎪ ⎩ 0
i≥ j otherwise
.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
367
Theorem 3.3. For small enough and without any restriction on t, the methods (3.10) and (3.11) preserve the maximum principle if
x ≤
2
|a|
.
Proof. For stability analysis, we can consider the method (3.10) in the general form
U(nt ) = Gn U(0), where
G = e (−2μ) SH. Now if
1
μ + ν ≥ 0
1
μ − ν ≥ 0,
and
2
(3.12)
2
we have n
n
Gn ∞ ≤ e (−2nμ) e (nμ+ 2 ν ) e (nμ− 2 ν ) = 1. But if x ≤ 2|a| , the conditions (3.12) are satisfied and the proof is complete. The proof for the method (3.11) is similar.
2
3.3. An upwind iterative splitting method for Burgers equation Consider the initial value problem for the one dimensional viscid Burgers’ equation
ut + uu x = u xx
x ∈ R,
u (x, 0) = u 0 (x)
x∈R
t > 0,
>0
.
(3.13)
In order to solve this nonlinear equation, we will use an iterative procedure simultaneous with upwind type of splitting method. For this, consider the space discretization of equation (3.13) as below:
(k+1)
dU j
=
dt
⎧ (k) −U j (k+1) (k+1) ⎪ ⎪ U j +1 − U j ⎪ x ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (k+1) (k+1) (k+1) (k) ⎪ + U j +1 if U j < 0 ⎨ + x2 U j −1 − 2U j
, k = 0, 1 , 2 , · · · , m (k) ⎪ −U j ⎪ (k+1) (k+1) ⎪ ⎪ U − U ⎪ j j −1 x ⎪ ⎪ ⎪ ⎪ ⎪ (k) ⎩ + U (k+1) − 2U (k+1) + U (k+1) if U j ≥ 0 j −1 j j +1 x2
(3.14)
where U(k) denotes the k−th iteration. We choose the answer of previous time step as starting vector U(0) in iteration procedure. For a tolerance δ , if U(k+1) − U(k) ∞ < δ , we stop the procedure and then the approximation at the time step is considered as U(k+1) . For solving the equation in each iteration by diagonal splitting method, consider the scheme (3.14) in the following form
dU( K +1) dt
= A(k) U(k+1) ,
(3.15)
(k)
where A( J −1)×( J −1) is a tridiagonal matrix whose j −th row is (k)
(k)
Uj Uj −2 ( 2, + , − ) x x2 x x x2
(k)
if U j < 0 and
(
x2
(k)
+
(k)
Uj −2 − , ) x x2 x x2
Uj
,
(k)
if U j ≥ 0. Now let A(k) = M(k) + D(k) + N(k) , where M(k) , D(k) and N(k) are three matrices which include three diagonals of matrix A . By using splitting process explained in Section 2.2, the solution U(k+1) of Equation (3.15) is defined as (k)
368
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
U(k+1) = e −2μ H(k) S(k) U0 ,
(3.16)
(N(k) t )
where S(k) = e is a lower triangular matrix computed by formula (2.5) and H(k) = e matrix computed by (2.4).
(M(k) t )
is an upper triangular
4. Numerical examples As comparative methods, we consider from [23] the upwind and Lax-Wendroff methods for one dimension advection Equation (3.1). The simplest scheme for solving (3.1) is upwind scheme as below
⎧ U j +1 − U j ⎪ ⎪ ⎨ −a U nj +1 − U nj x =
⎪ t ⎪ ⎩ −a U j − U j −1 x
if a < 0 (4.1)
. if a > 0
In this method an explicit computation shows that the upwind scheme preserves the maximum principle under the CFL condition ν ≤ 1 (a > 0). On the other hand the upwind method is first order accurate and also first order convergence along a refinement path which satisfies the CFL condition everywhere [23]. Furthermore, we consider the Lax-Wendroff method for the Problem (3.1) as below
U nj +1 − U nj
t
= −a
U nj+1 − U nj−1
+
2 x
1 2
a2 t
U n − 2U n + U n j +1 j j −1
x2
(4.2)
or
U nj +1 =
1 2
1
ν (1 + ν )U nj−1 + (1 − ν 2 )U nj − ν (1 − ν )U nj+1 . 2
(4.3)
This method is of order 2 in time and space and does not preserve the maximum principle. On the other hand, we will use the MATLAB function pdepe for further investigation of efficiency of the proposed methods. The function pdepe solves initial-boundary value problems for small systems of PDEs in one space variable of the form
c (x, t , u , u x ) ut = x(−m) xm f (x, t , u , u x )
x
+ s(x, t , u , u x ),
where f (x, t , u , u x ) is a flux, s(x, t , u , u x ) is a source term and m must be 0, 1, or 2. In this function, the ODEs resulting from discretization in space are integrated to obtain approximate solutions at times specified in the input array. The time integration is done with an ODE solver that selects both the time step and formula dynamically [26]. 4.1. Example 1 Consider the one dimensional advection equation
ut + au x = 0, in (−2, 10) with initial condition u (x, 0) = exp (−(x − 7)2 ) and the boundary conditions u (10, t ) = 0. Let a = −1, the upwind method for this equation is stable if
|ν | ≤ 1 . Here we compare upwind type splitting method (3.4) and upwind method (4.1). In Fig. 1, we have x = 0.01, t = 0.018 and |ν | = 1.8. It is clear that the upwind method is not stable (|ν | = 1.8 is outside of region of stability), but the proposed method is stable and has an acceptable quality compared with the result of applying pdepe to the desired problem. 4.2. Example 2 Here we solve the one dimensional advection equation
ut + au x = 0, by the Lie-Trotter splitting method (3.10) and Strang splitting method (3.11) in domain (−2, 10) with the boundary conditions u (−2, t ) = u (10, 0) = 0. Let a = 1 and the initial condition as a square pulse as below:
u (x, 0) =
4 if − 1.2 ≤ x ≤ 0 0
otherwise
.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
369
Fig. 1. Numerical solution of one dimensional advection equation with initial condition u (x, 0) = exp (−(x − 7)2 ) with upwind method (4.1), upwind type of splitting method (3.4) and pdepe function with |ν | = 1.8 at time t = 2.
Fig. 2 shows the solution of Lax-Wendroff (4.3), Lie-Trotter splitting method (3.10) and pdepe function with x = 0.01, t = 0.018 (ν = 1.8) and = 12 a2 t at time t = 1. As the figures show, the proposed splitting method, unlike the pdepe function, is stable and preserves the maximum principle with ν = 1.8, but the Lax-Wendroff method is not stable. In Fig. 3, we take x = 0.02, t = 0.01 (ν = 0.5 < 1) and = 12 a2 t = 0.005. As expected, the Lax-Wendroff method does not preserve the maximum principle, which is the same in pdepe function, but the situation in the proposed method is much better. In order to test accuracy of the Strang splitting method (3.11), we used this method to solve the above problem with = 12 a2 t and smooth initial condition u (x, 0) = exp(−(x − 2)2 ). Table 1 shows the approximate error and temporal convergence rate of this method in different time steps with x = 0.05 at t = 0.5. This table shows that, as expected, the proposed Strang splitting method (3.11) has a second order accuracy in time. The numerical evaluation of order of time convergence is done as below E t 1
p≈
ln( E
t 2
)
t1 ln( t )
,
(4.4)
2
where the error associated with step size t is measured as the difference between this solution, ut , and the solution u2t for time step size 2t, i.e.,
370
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 2. Numerical solution of one dimensional advection equation with a square pulse initial condition with Lax-Wendroff method (4.3), Lie-Trotter splitting method (3.10) and pdepe with = 12 a2 t and ν = 1.8 > 1 at time t = 1.
E t = ut − u2t ∞ ,
(4.5)
and ut and u2t are numerical solutions obtained from step sizes t and 2t respectively. 4.3. Example 3 Consider the two dimensional advection equation
u t + a 1 u x + a 2 u y = 0, in domain (−1, 1) × (−1, 1) with initial condition u (x, 0) = sin(
(4.6)
π
x) sin(
π
y ). 2 2 Let a1 = −1, a2 = 1, t = 0.05 and x = y = 0.05. Fig. 4 shows the results obtained from applying the second order Strang splitting method (3.11) on the desired equation in different times. It can be seen that the initial wave moves to the corner of the domain without any oscillations.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
371
Fig. 3. Numerical solution of one dimensional advection equation with a square pulse initial condition with Lax-Wendroff method (4.3), the Strang splitting method (3.11) and pdepe function with = 12 a2 t and ν = 0.5 at time t = 1. Table 1 Approximation of error and temporal convergence rate of Strang splitting method (3.11) in Example 2 with x = 0.05 at t = 0.5.
t
E t
Convergence order
2.5e−2 1.25e−2 6.25e−3 3.125e−3 1.5625e−3 7.8125e−4
1.15622e−002 4.2488e−003 1.51662e−003 5.25546e−004 1.5666e−004 4.2857e−005
– 1.4443 1.48626 1.52906 1.74625 1.87008
4.4. Example 4 In this example we use upwind iterative splitting method (3.16) for solving the viscid Burgers’ equation
ut + uu x = u xx ,
x ∈ [a, b],
t > 0,
> 0.
Set δ = 10e−4. The results show that the convergence of this iteration process is very fast, so far as the number of iterations does not exceed 2 per time step.
372
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 4. Numerical solution of Example 3 in different times. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)
Fig. 5 shows the numerical solution of viscid Burgers’ equation in [−0.5, 3] with Gaussian initial function u (x, 0) = 2
e (−2(x−1) ) , x = 0.1, t = 0.01 and = 0.01 after 100 time steps. In Fig. 6, we test the proposed method with Gaussian initial function 2 u (x, 0) = e (−2(x−1) ) ,
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 4. (continued)
cases (a) and (b), Riemann initial function
u (x, 0) =
uL = 1 uR = 0
x<0 x≥0
cases (c) and (d) and piecewise continuous initial function
373
374
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 5. Numerical solution of viscid Burgers’ equation with upwind type iterative splitting method and Gaussian initial function at t = 1.
⎧ x<0 ⎨ u=1 1−x 0≤x≤1 u (x, 0) = ⎩ uR = 0 x≥1 cases (e) and (f) with different values for
. Here we have used x = 0.1 and t = 0.001.
4.5. Example 5 Here we want to apply proposed splitting methods on linear system of first order hyperbolic equations in the form
Ut +AUx = 0,
(4.7)
where U(x, t) is a vector in Rn and A is a matrix with real eigenvalues and a full set of eigenvectors. For example consider the system
u v
+ t
a11 a21
a12 a22
u v
= 0,
we can rewrite (4.8) in the form
u v
=− t
a11 0
0 a22
(4.8)
x
u v
− x
0 a12 0 0
u v
− x
0 a21
0 0
u v
x
and solve any of this subsystems by using the splitting ideas. Two hyperbolic equations corresponding to first subsystem
ut = −a11 u x v t = −a22 v x
can be solved by proposed splitting methods introduced in Section 3. Now consider the equations
ut = −a12 v x vt = 0
corresponding to the second subsystem. Equation v t = 0 implies that v is constant during a time step and the second equation
ut = −a12 v x can be solved by using finite difference methods like
u ∗∗ j =
−a12 t ∗ v j +1 − v ∗j −1 + u ∗j , 2 x
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
375
Fig. 6. Numerical solution of upwind type iterative splitting method on viscid Burgers’ equation with Gaussian initial function and = 0.1 (a), Gaussian initial function and = 0.01 (b), Riemann initial function and = 0.1 (c), Riemann initial function and = 0.01 (d), piecewise continuous initial function and = 0.1 (e) and piecewise continuous initial function and = 0.01 (f) with t = 0.001 and x = 0.1.
376
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 7. Numerical solution of system (4.8) with proposed splitting method with coefficients matrix (4.9).
where u ∗ and v ∗ are the solutions obtained from solving the first subsystem. The third system can also be solved in the same way. Let x = 0.01, t = 0.001 and
A=
−1 1 1
0
.
(4.9)
Fig. 7 shows the solutions u and v of the system (4.8) in domain (−1, 1) × (−1, 1) with proposed method and with the initial conditions
u (x, 0) = exp(−80(x + 0.33)2 ) v (x, 0) = exp(−80(x − 0.66)2 ), and homogenous Dirichlet boundary conditions. On the other hand Fig. 8 shows the solutions u and v of the system (4.8) with
A=
0 1 1 0
(4.10)
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
377
Fig. 8. Numerical solution of system (4.8) with proposed splitting method with coefficients matrix (4.10).
and with conditions and parameters same as before. It should be noted that this system is derived from the familiar second order wave equation w tt = w xx using the variables u = w t and v = − w x . 5. Conclusion In the current research, some numerical methods are obtained for solving hyperbolic partial differential equations. An upwind time splitting method with unconditionally maximum principle preserving property is presented, also damping of its solution is studied. By adding an artificial diffusion to the right hand of the advection equation a second order method with relatively mild condition for stability is derived. These results have important role in the numerical treatment of hyperbolic partial differential equations. As a nonlinear hyperbolic equation, Burgers’ equation is solved using the new ideas and an iterative procedure is implemented. Numerical examples validate the theoretical results and efficiency of the proposed methods. Splitting of the differential equation in certain parts and solving these parts exactly may be the main reason of the excellence of the present methods over the other existing numerical methods. The ides is used for solving two dimensional and system of hyperbolic partial differential equations.
378
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Acknowledgement The authors are very grateful to reviewers for carefully reading the paper and for their constructive comments and suggestions. This research was in part supported by a grant from IPM (No. 98650416). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
H. Bateman, Some recent researches on the motion of fluids, Mon. Weather Rev. 43 (1915) 163–170. G.P.E. Benton, A table of solutions of the one-dimensional Burgers equation, Q. Appl. Math. 30 (1972) 195–212. S. Blanes, F. Casas, A Concise Introduction to Geometric Numerical Integration, CRC Press, 2016. J.M. Burgers, A mathematical model illustrating the theory of turbulence, Adv. Appl. Mech. 1 (1948) 171–199. J. Cai, H. Liang, C. Zhang, Efficient high-order structure-preserving methods for the generalized Rosenau-type equation with power law nonlinearity, Commun. Nonlinear Sci. Numer. Simul. 59 (2018) 122–131. M. Caliari, A. Ostermann, C. Piazzola, A splitting approach for the magnetic Schrödinger equation, J. Comput. Appl. Math. 316 (2017) 74–85. S.A. Chin, Understanding Saul’yev-type unconditionally stable schemes from exponential splitting, Numer. Methods Partial Differ. Equ. 30 (6) (2014) 1961–1983. J.D. Cole, On a quasilinear parabolic equation occurring in aerodynamics, Q. Appl. Math. 9 (1951) 225–236. M. Dehghan, A. Taleei, A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients, Comput. Phys. Commun. 181 (2010) 43–51. J. Douglas Jr., H.H. Rachford Jr., On the numerical solution of the heat conduction problems in two and three variables, Trans. Am. Math. Soc. 82 (1956) 421–439. S.A. Dyachenko, A. Zlotnik, A.O. Korotkevich, M. Chertkov, Operator splitting method for simulation of dynamic flows in natural gas pipeline networks, Phys. D: Nonlinear Phenom. 361 (2017) 1–11. ´ E.G. Dyakonov, Difference schemes of second order accuracy with a splitting operator for parabolic equations without mixed partial derivatives, Ž. Vyˇcisl. Mat. Mat. Fiz., Mosc. 4 (1964) 935–941. J. Geiser, Additive and iterative splitting methods for multiscale and multiphase coupled problems, J. Coupled Syst. Multiscale Dyn. 4 (4) (2016) 271–291. J. Geiser, Iterative Splitting Methods for Differential Equations, CRC Press, 2017. J. Geiser, J.L. Huesob, E. Martinez, New versions of iterative splitting methods for the momentum equation, J. Comput. Appl. Math. 309 (2017) 359–370. P. González de Alaiza Martínez, M.E. Vázquez-Cendón, Operator-Splitting on Hyperbolic Balance Laws, Springer International Publishing, Cham, 2014, pp. 279–287. E. Hairer, Ch. Lubich, G. Wanner, Geometric Numerical Integration, second edition, Springer Series in Computational Mathematics, 2000. E. Hopf, The partial differential equation ut + uu x = μu xx , Commun. Pure Appl. Math. 3 (1950) 201–230. R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1990. H.G. Lee, A second-order operator splitting Fourier spectral method for fractional-in-space reaction-diffusion equations, J. Comput. Phys. 333 (2018) 345–403. S. MacNamara, G. Strang, Operator splitting, in: Splitting Methods in Communication, Imaging, Science and Engineering, 2016, pp. 95–114. G. Marchuk, Some applications of splitting-up methods to the solution of mathematical physics problems, Apl. Mat. 13 (1968) 103–132. K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, second edition, Cambridge University Press, 2005. D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 28–41. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. TheMathWorks Inc. (2017). H.F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10 (1959) 545–551. G. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, 1974. J. Xie, Z. Zhang, D. Liang, A conservative splitting difference scheme for the fractional-in-space Boussinesq equation, Appl. Numer. Math. 143 (2019) 61–74, https://doi.org/10.1016/j.apnum.2019.03.013.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Some splitting methods for hyperbolic PDEs Rasool Hosseini a , Mehdi Tatari a,b,∗ a b
Department of Mathematical Sciences, Isfahan University of Technology, Isfahan, 84156-83111, Iran Institute for Research in Fundamental Sciences (IPM), P.O. Box: 19395-5746, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 10 January 2019 Received in revised form 18 May 2019 Accepted 4 July 2019 Available online 9 July 2019 Keywords: Advection equation Burgers’ equation Maximum principle Damping of solution
a b s t r a c t In this work, a new splitting technique is implemented for solving hyperbolic PDEs. As the main result, the new methods preserve the maximum principle unconditionally or with a mild condition on discretization parameters in comparison with well known methods. Damping of numerical solution in time evolution is investigated. For numerical solution of the Burgers’ equation, as a nonlinear problem, an iterative method based on the new splitting technique is presented. Efficiency of the new methods is examined via numerical examples. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Consider the advection equation
ut + au x = 0. This equation describes the passive advection of some scalar field u (x, t ) carried along by a flow of constant speed a. This equation must surely be the simplest of all partial differential equations. Yet to approximate it well on a fixed (x, t )-mesh is a far from trivial problem that is still under active discussion in the numerical analysis literature [23]. The exact solution of this equation is constant along the characteristic lines associated to this equation. Since finding characteristic lines are complicated especially for system of hyperbolic equations, usually use of numerical methods is preferred. Some well-known numerical methods for solving this equation are upwind and Lax-Wendroff methods. In numerical methods for solving advection equation, preservation of maximum principle is an important subject. The upwind method preserves this property under the CFL condition while Lax-wendroff always does not preserve the property. As an other important issue in this kind of equations, damping of the solution should be studied. In this paper some new methods are presented and their behaviour for solving advection equation are discussed. Burgers’ equation or Bateman-Burgers’ equation is a nonlinear hyperbolic partial differential equation occurring in various areas of applied mathematics, such as fluid mechanics, dispersive water and traffic flow. Also the Burgers’ equation appears often as a simplification of a more complicated model like Navier Stokes equations. The equation was first introduced by Harry Bateman [1] and later studied by Johannes Martinus Burgers’ [4]. For a given field u (x, t ) and kinematic viscosity , the general form of Burgers’ equation (also known as viscous Burgers’ equation or advective form of the Burgers’ equation) in one space dimension is the dissipative system:
*
Corresponding author. E-mail addresses: [email protected] (R. Hosseini), [email protected] (M. Tatari).
https://doi.org/10.1016/j.apnum.2019.07.005 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.
362
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
ut + uu x = u xx . When
= 0, Burgers’ equation becomes the inviscid Burgers’ equation
ut + uu x = 0, which is a prototype for conservation equations that can develop discontinuities as shock waves. The inviscid Burgers’ equation is a conservation equation, more generally a first order nonlinear hyperbolic equation. The solution of this equation with the initial condition u (x, 0) = f (x) can be constructed by many of numerical methods like the method of characteristics [19,28]. In recent years, computing the solution of Burgers’ equation has attracted a lot of attention and various methods have been introduced by researchers to solve this equation. Although the exact solution of this equation with some simple initial conditions is presented by Hopf [18], Benton and Platzman [2] and Cole [8], but in most situations, this equation can only be solved numerically. To see a complete overview of numerical methods for solving Burgers’ equation see [19]. Splitting methods introduced by Trotter [27] and Strang [25] are one of the most widely used methods among geometric integrators, which are often used to overcome the complexity of computing problems, constructing very accurate high order numerical algorithms and extension of the region of stability of methods. These ideas are used to overcome the computational complexity which is arisen in numerical solution of higher dimensional problems [10,12,24]. In [11] a second order accurate operator splitting method with unconditional stability is presented to solve a system of nonlinear hyperbolic PDEs that is obtained from the simulation of flows of isothermal compressible natural gas over transmission pipelines. The idea of the splitting for multidimensional linear hyperbolic systems and for Burgers’ equation as a nonlinear hyperbolic equation is introduced in [16]. In [20] a second-order operator splitting Fourier spectral method is presented for solving fractional in space reaction-diffusion equations. Splitting approaches for Schrödinger equation are given in [6,9]. Also, iterative splitting methods are one of the suitable methods for solving nonlinear equations [14,15]. The idea of splitting can be applied to wide range of problems [3,17]. For more information about splitting methods see [5,22,29]. In [7], it is shown that Saul’yev-type schemes can be derived from the exponential splitting of the semi-discretized equation and the splitting idea in which it is used is to rewrite the resulting matrix from the semi-discretization as a sum of matrices that have only a single, nonvanishing 2 × 2 matrix along the diagonal. Some splitting methods are introduced in [21] and recent developments are presented. New additive and iterative splitting schemes and algorithms are constructed in [13]. In this work some splitting methods are outlined for hyperbolic equations. Then a new idea which is different from existing schemes is presented. It is based on splitting of semi-discretized equations along the diagonal. The reminder of this paper is organized as follows. In Section 2, we present preliminaries about splitting methods. In Section 3, some proper splitting methods are implemented for solving advection and Burgers’ equations. Numerical results are investigated in Section 4 and finally a conclusion is drawn in Section 5. 2. Splitting methods In this section, basic concepts of some splitting methods that are used in the next section are introduced [3,17]. 2.1. Preliminary Consider the differential equation
dy(t ) dt
= A (y(t )) + B (y(t )),
y(0) = y0 ,
(2.1)
where y is a vector or possibly a matrix and A and B are operators on y. In a splitting method for solving (2.1) numerically, it is needed to solve exactly (or numerically) the following two subproblems
dy(t ) dt
= A (y(t )),
dy(t ) dt
= B (y(t )).
(2.2)
ϕtA and ϕtB the exact t flows for each subproblems in (2.2). This means that, for example, the solution of dy(t ) subproblem dt = A (y(t )) in time t is obtained by y (t ) = ϕtA ( y (t 0 )). In this way, we have the first order Lie-Trotter splitting Denote by
methods that one is the adjoint of the other
h∗ = ϕhB ◦ ϕhA , h = ϕhA ◦ ϕhB , where ◦ (composition operator) is an operation that takes two functions f and g and produces a function h such that h(x) = g ( f (x)). By composing the Lie-Trotter splitting method and its adjoint with a half time step, the Strang splitting formulae
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
363
h[s] = ϕhA/2 ◦ ϕhB ◦ ϕhA/2 ,
h[s] = ϕhB/2 ◦ ϕhA ◦ ϕhB/2 , will be obtained that are second order time symmetric methods. If Problem (2.1) can be decomposed into more than two parts as
dy(t ) dt
= A [1] (y(t )) + · · · + A [m] (y(t )),
and the exact flows
ϕh[i] correspond to each split part
dy(t ) dt
= A [i ] (y(t )), i = 1, ..., m be available, then the scheme
χh = ϕh[m] ◦ · · · ◦ ϕh[2] ◦ ϕh[1] , and its adjoint map
χh∗ = χ−−h1 = ϕh[1] ◦ ϕh[2] ◦ · · · ◦ ϕh[m] , are first order approximations (Lie-Trotter type) and the compositions
χh/2 ◦ χh∗/2
χh∗/2 ◦ χh/2 ,
or
provides second order methods (Strang type). 2.2. New splitting methods Here we introduce a new splitting technique for system of ordinary differential equations. For this, first consider the following preposition which will be used later. Preposition 2.1. Any strictly triangular matrix A ∈ Rn×n (a triangular matrix with zero on main diagonal) is nilpotent, in other words there is an integer m ≤ n such that Am = 0. Proof. By induction on n.
2
Now consider a linear system of ordinary differential equations
dx(t ) dt
= Ax(t ),
x(t 0 ) = x0 ,
(2.3)
with coefficients matrix A = (ai j ) ∈ Rn×n . By splitting matrix A into its 2n − 1 diagonals
A = L1 + · · · + Ln−1 + D + U1 + · · · + Un−1 , System (2.3) can be split as follows
dx(t ) dt
= L1 x(t ) + · · · + Ln−1 x(t ) + Dx(t ) + U1 x(t ) + · · · + Un−1 x(t ).
The exact flow of the system
dx(t ) dt
= Dx(t ),
is
x(t ) = e (Dt ) x0 = D(t )x0 ,
D(t ) = diag(e (a11 t ) , . . . , e (ann t ) ),
where D is a diagonal matrix that include main diagonal of matrix A. Let
Ek (t ) = e (Uk t ) =
∞ tj (U K ) j , k = 1, . . . , n − 1, j! j =1
∞ tj Sk (t ) = e (Lk t ) = (L K ) j , k = 1, . . . , n − 1. j! j =1
By Preposition 2.1 matrices Uk and Lk are nilpotent and Ek (t ) and Sk (t ) are calculable finitely, such that Ek (t ) is an upper triangular matrix with 1 on the main diagonal and
364
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
j−k i
s=1 ai +(s−1)k,i +sk
(Ek )i j =
t
j −i k
for
i ( j− )! k
j > i & j − i = mk, m = 1, . . . ,
n−1
k
,
(2.4)
and Sk (t ) is a lower triangular matrix with 1 on the main diagonal and
i− j k
(Sk )i j =
s=1 a j +sk, j +(s−1)k t
i− j k
for
( i −k j )!
j < i & i − j = mk, m = 1, . . . ,
n−1 k
.
(2.5)
dx(t ) = Uk x(t ) and = Lk x(t ) are x(t ) = Ek (t )x0 and x(t ) = Sk (t )x0 respecdt dt tively for k = 1, . . . , n − 1. Now let L(t ) = S1 (t ) · · · Sn−1 (t ) and U(t ) = E1 (t ) · · · En−1 (t ), then the Lie-Trotter splitting method for the system (2.3) is Therefore, the exact flows of the systems
dx(t )
x(t ) = L(t ) D(t )U(t )x0 ,
(2.6)
where the adjoint map of this method is
x(t ) = U(t ) D(t )L(t )x0 .
(2.7)
By composing (2.6) and (2.7) with half time step, we have the Strang splitting formula in the form
x(t ) =L(
=L(
t t t t t t )D( )U( )U( ) D( )L( )x0 2 t 2
) D(
2 t 2
2
)U(t ) D(
2 t 2
)L(
2 t 2
(2.8)
2
)x0 .
(2.9)
This idea specially is useful for solving systems that A is a band matrix. These systems appear in spatial discretization of partial differential equations. 3. New splitting methods and hyperbolic PDEs In this section, pure advection equation,
ut + au x = 0,
(3.1)
is considered, where u is the advection velocity and nonzero a ∈ R is advection coefficient. 3.1. An upwind type splitting method Consider nodes x j in [c , d] with x j = c + j x, j = 1, ..., J where x =
d−c J
, J ∈ N . The approximate solution is then
denoted by
U j ≈ u (x j , t ). As an upwind type method, we use a backward difference in space if a is positive and a forward difference if a is negative [23]
⎧ −a ⎪ U j +1 − U j if a < 0 ⎨ dU j x = . ⎪ dt ⎩ −a U − U if a > 0 j j −1 x Set
ν =a
t . Consider scheme (3.2) as x
dU dt
(3.2)
= AU,
(3.3)
where
A( J −1)×( J −1) = upperbidiag( or
a −a , ), x x
if a < 0,
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
A( J −1)×( J −1) = lowerbidiag(
a −a , ), x x
if a > 0.
Now assume a > 0, write A = D + N, where D = diag( follows
dU dt
365
−a a ) and N = lowerdiag( ). Therefore, System (3.3) can be split as x x
= DU + NU.
By using the splitting procedure explained in previous section, we have the upwind type splitting method for the Problem (3.1) as
U(t ) = exp(−ν )EU0 ,
(3.4)
where E is a lower triangular matrix with elements
⎧ i− j ⎪ ⎨ (ν ) e i j = (i − j )! ⎪ ⎩ 0
j≤i
(3.5)
.
otherwise
This method preserves the maximum principle property, more precisely we have
E∞ = 1 + ν + ≤
∞ (ν )i i =0
i!
(ν )2 (ν ) J − 2 + ··· + 2! ( J − 2)! = exp(ν ),
and if we set G = exp(−ν )E, then
Gn ∞ ≤ exp(−nν )En∞ , and this concludes that
U(nt )∞ ≤ U(0)∞ . Damping of the solution is an important feature of numerical methods for solving advection equation. More precisely, since there is no damping in the analytical solution, we would like to have G∞ = 1. In the next theorem, damping of the new method is investigated. Theorem 3.1. Maximum damping in time step n + 1 of upwind type splitting method (3.4) is bounded by
−
ν J −1 ( J − 1)!
U(nt )∞ .
Proof. We know
G∞ = exp(−ν )E∞ = exp(−ν )(exp(ν ) − =1−
exp(ξ )ν J −1
exp(ξ − ν )ν J −1
( J − 1)!
( J − 1)!
)
,
where 0 ≤ ξ ≤ ν . Therefore, maximum damping of amplitude matrix is
−
ν J −1 ( J − 1)!
.
Thus, maximum damping in time step n + 1 is
−
ν J −1 ( J − 1)!
U(nt )∞ . 2
366
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
If a < 0, the proposed splitting method will be as
U(t ) = exp(ν )SU0 ,
(3.6)
where S is an upper triangular matrix with
⎧ (−ν ) j −i ⎪ ⎪ ⎨ ( j − i )! si j = ⎪ ⎪ ⎩ 0
j≥i (3.7)
. otherwise
Here, results are similar to the case a > 0. Remark 3.2. For a > 0, notice that elements on each diagonal of G are a multiplier of elements on its main diagonal. More precisely, element on main diagonal of G are exp(−ν ), and elements on the next diagonals are etc. Similar strategy should be considered when a < 0.
ν exp(−ν ),
ν 2 exp(−ν ) 2!
and
3.2. A second order time splitting method As we discussed above, the upwind type method resulting from the splitting scheme is first order accurate. One simple way to achieve second order accuracy on the advection equation is to use an artificial diffusion term u xx with coefficient > 0. Applying this directly to the (3.1) and use a second order discretization on space, gives the second order scheme
dU j dt
=
−a U j +1 − U j −1 + U j +1 − 2U j + U j −1 . 2 x
(3.8)
Consider this scheme as
dU dt
= AU,
(3.9)
where
A( J −1)×( J −1) = tridiag(
x2
+
a
,
−2
,
2x x2 x2
−
a 2 x
).
a a 2 Now Write A = M + D + N, where M = upperdiag( 2 + 2 ), D = diag( − ) and N = lowerdiag( x2 − ). Therefore, x x x2 2x System (3.9) can be split as follows
dU dt
= MU + DU + NU.
Let μ = t2 and ν = axt . By using the explained splitting procedure and introducing F(t ) = e (−μ) S and F(t ) = e (−μ) H, x the Lie-Trotter type of proposed splitting method for the Problem (3.8) is
U(t ) = F(t ) F(t )U(0),
(3.10)
and the Strang type formula is
U(t ) = F(
t t )F(t )F( )U(0), 2
(3.11)
2
where H is an upper triangular matrix with elements
⎧ 1 ( j −i ) ⎪ ⎪ ⎨ (μ − 2 ν ) hi j = ( j − i )! ⎪ ⎪ ⎩ 0
j≥i
,
otherwise
and S is a lower triangular matrix with elements
⎧ 1 (i − j ) ⎪ ⎪ ⎨ (μ + 2 ν ) si j = (i − j )! ⎪ ⎪ ⎩ 0
i≥ j otherwise
.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
367
Theorem 3.3. For small enough and without any restriction on t, the methods (3.10) and (3.11) preserve the maximum principle if
x ≤
2
|a|
.
Proof. For stability analysis, we can consider the method (3.10) in the general form
U(nt ) = Gn U(0), where
G = e (−2μ) SH. Now if
1
μ + ν ≥ 0
1
μ − ν ≥ 0,
and
2
(3.12)
2
we have n
n
Gn ∞ ≤ e (−2nμ) e (nμ+ 2 ν ) e (nμ− 2 ν ) = 1. But if x ≤ 2|a| , the conditions (3.12) are satisfied and the proof is complete. The proof for the method (3.11) is similar.
2
3.3. An upwind iterative splitting method for Burgers equation Consider the initial value problem for the one dimensional viscid Burgers’ equation
ut + uu x = u xx
x ∈ R,
u (x, 0) = u 0 (x)
x∈R
t > 0,
>0
.
(3.13)
In order to solve this nonlinear equation, we will use an iterative procedure simultaneous with upwind type of splitting method. For this, consider the space discretization of equation (3.13) as below:
(k+1)
dU j
=
dt
⎧ (k) −U j (k+1) (k+1) ⎪ ⎪ U j +1 − U j ⎪ x ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ (k+1) (k+1) (k+1) (k) ⎪ + U j +1 if U j < 0 ⎨ + x2 U j −1 − 2U j
, k = 0, 1 , 2 , · · · , m (k) ⎪ −U j ⎪ (k+1) (k+1) ⎪ ⎪ U − U ⎪ j j −1 x ⎪ ⎪ ⎪ ⎪ ⎪ (k) ⎩ + U (k+1) − 2U (k+1) + U (k+1) if U j ≥ 0 j −1 j j +1 x2
(3.14)
where U(k) denotes the k−th iteration. We choose the answer of previous time step as starting vector U(0) in iteration procedure. For a tolerance δ , if U(k+1) − U(k) ∞ < δ , we stop the procedure and then the approximation at the time step is considered as U(k+1) . For solving the equation in each iteration by diagonal splitting method, consider the scheme (3.14) in the following form
dU( K +1) dt
= A(k) U(k+1) ,
(3.15)
(k)
where A( J −1)×( J −1) is a tridiagonal matrix whose j −th row is (k)
(k)
Uj Uj −2 ( 2, + , − ) x x2 x x x2
(k)
if U j < 0 and
(
x2
(k)
+
(k)
Uj −2 − , ) x x2 x x2
Uj
,
(k)
if U j ≥ 0. Now let A(k) = M(k) + D(k) + N(k) , where M(k) , D(k) and N(k) are three matrices which include three diagonals of matrix A . By using splitting process explained in Section 2.2, the solution U(k+1) of Equation (3.15) is defined as (k)
368
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
U(k+1) = e −2μ H(k) S(k) U0 ,
(3.16)
(N(k) t )
where S(k) = e is a lower triangular matrix computed by formula (2.5) and H(k) = e matrix computed by (2.4).
(M(k) t )
is an upper triangular
4. Numerical examples As comparative methods, we consider from [23] the upwind and Lax-Wendroff methods for one dimension advection Equation (3.1). The simplest scheme for solving (3.1) is upwind scheme as below
⎧ U j +1 − U j ⎪ ⎪ ⎨ −a U nj +1 − U nj x =
⎪ t ⎪ ⎩ −a U j − U j −1 x
if a < 0 (4.1)
. if a > 0
In this method an explicit computation shows that the upwind scheme preserves the maximum principle under the CFL condition ν ≤ 1 (a > 0). On the other hand the upwind method is first order accurate and also first order convergence along a refinement path which satisfies the CFL condition everywhere [23]. Furthermore, we consider the Lax-Wendroff method for the Problem (3.1) as below
U nj +1 − U nj
t
= −a
U nj+1 − U nj−1
+
2 x
1 2
a2 t
U n − 2U n + U n j +1 j j −1
x2
(4.2)
or
U nj +1 =
1 2
1
ν (1 + ν )U nj−1 + (1 − ν 2 )U nj − ν (1 − ν )U nj+1 . 2
(4.3)
This method is of order 2 in time and space and does not preserve the maximum principle. On the other hand, we will use the MATLAB function pdepe for further investigation of efficiency of the proposed methods. The function pdepe solves initial-boundary value problems for small systems of PDEs in one space variable of the form
c (x, t , u , u x ) ut = x(−m) xm f (x, t , u , u x )
x
+ s(x, t , u , u x ),
where f (x, t , u , u x ) is a flux, s(x, t , u , u x ) is a source term and m must be 0, 1, or 2. In this function, the ODEs resulting from discretization in space are integrated to obtain approximate solutions at times specified in the input array. The time integration is done with an ODE solver that selects both the time step and formula dynamically [26]. 4.1. Example 1 Consider the one dimensional advection equation
ut + au x = 0, in (−2, 10) with initial condition u (x, 0) = exp (−(x − 7)2 ) and the boundary conditions u (10, t ) = 0. Let a = −1, the upwind method for this equation is stable if
|ν | ≤ 1 . Here we compare upwind type splitting method (3.4) and upwind method (4.1). In Fig. 1, we have x = 0.01, t = 0.018 and |ν | = 1.8. It is clear that the upwind method is not stable (|ν | = 1.8 is outside of region of stability), but the proposed method is stable and has an acceptable quality compared with the result of applying pdepe to the desired problem. 4.2. Example 2 Here we solve the one dimensional advection equation
ut + au x = 0, by the Lie-Trotter splitting method (3.10) and Strang splitting method (3.11) in domain (−2, 10) with the boundary conditions u (−2, t ) = u (10, 0) = 0. Let a = 1 and the initial condition as a square pulse as below:
u (x, 0) =
4 if − 1.2 ≤ x ≤ 0 0
otherwise
.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
369
Fig. 1. Numerical solution of one dimensional advection equation with initial condition u (x, 0) = exp (−(x − 7)2 ) with upwind method (4.1), upwind type of splitting method (3.4) and pdepe function with |ν | = 1.8 at time t = 2.
Fig. 2 shows the solution of Lax-Wendroff (4.3), Lie-Trotter splitting method (3.10) and pdepe function with x = 0.01, t = 0.018 (ν = 1.8) and = 12 a2 t at time t = 1. As the figures show, the proposed splitting method, unlike the pdepe function, is stable and preserves the maximum principle with ν = 1.8, but the Lax-Wendroff method is not stable. In Fig. 3, we take x = 0.02, t = 0.01 (ν = 0.5 < 1) and = 12 a2 t = 0.005. As expected, the Lax-Wendroff method does not preserve the maximum principle, which is the same in pdepe function, but the situation in the proposed method is much better. In order to test accuracy of the Strang splitting method (3.11), we used this method to solve the above problem with = 12 a2 t and smooth initial condition u (x, 0) = exp(−(x − 2)2 ). Table 1 shows the approximate error and temporal convergence rate of this method in different time steps with x = 0.05 at t = 0.5. This table shows that, as expected, the proposed Strang splitting method (3.11) has a second order accuracy in time. The numerical evaluation of order of time convergence is done as below E t 1
p≈
ln( E
t 2
)
t1 ln( t )
,
(4.4)
2
where the error associated with step size t is measured as the difference between this solution, ut , and the solution u2t for time step size 2t, i.e.,
370
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 2. Numerical solution of one dimensional advection equation with a square pulse initial condition with Lax-Wendroff method (4.3), Lie-Trotter splitting method (3.10) and pdepe with = 12 a2 t and ν = 1.8 > 1 at time t = 1.
E t = ut − u2t ∞ ,
(4.5)
and ut and u2t are numerical solutions obtained from step sizes t and 2t respectively. 4.3. Example 3 Consider the two dimensional advection equation
u t + a 1 u x + a 2 u y = 0, in domain (−1, 1) × (−1, 1) with initial condition u (x, 0) = sin(
(4.6)
π
x) sin(
π
y ). 2 2 Let a1 = −1, a2 = 1, t = 0.05 and x = y = 0.05. Fig. 4 shows the results obtained from applying the second order Strang splitting method (3.11) on the desired equation in different times. It can be seen that the initial wave moves to the corner of the domain without any oscillations.
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
371
Fig. 3. Numerical solution of one dimensional advection equation with a square pulse initial condition with Lax-Wendroff method (4.3), the Strang splitting method (3.11) and pdepe function with = 12 a2 t and ν = 0.5 at time t = 1. Table 1 Approximation of error and temporal convergence rate of Strang splitting method (3.11) in Example 2 with x = 0.05 at t = 0.5.
t
E t
Convergence order
2.5e−2 1.25e−2 6.25e−3 3.125e−3 1.5625e−3 7.8125e−4
1.15622e−002 4.2488e−003 1.51662e−003 5.25546e−004 1.5666e−004 4.2857e−005
– 1.4443 1.48626 1.52906 1.74625 1.87008
4.4. Example 4 In this example we use upwind iterative splitting method (3.16) for solving the viscid Burgers’ equation
ut + uu x = u xx ,
x ∈ [a, b],
t > 0,
> 0.
Set δ = 10e−4. The results show that the convergence of this iteration process is very fast, so far as the number of iterations does not exceed 2 per time step.
372
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 4. Numerical solution of Example 3 in different times. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)
Fig. 5 shows the numerical solution of viscid Burgers’ equation in [−0.5, 3] with Gaussian initial function u (x, 0) = 2
e (−2(x−1) ) , x = 0.1, t = 0.01 and = 0.01 after 100 time steps. In Fig. 6, we test the proposed method with Gaussian initial function 2 u (x, 0) = e (−2(x−1) ) ,
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 4. (continued)
cases (a) and (b), Riemann initial function
u (x, 0) =
uL = 1 uR = 0
x<0 x≥0
cases (c) and (d) and piecewise continuous initial function
373
374
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 5. Numerical solution of viscid Burgers’ equation with upwind type iterative splitting method and Gaussian initial function at t = 1.
⎧ x<0 ⎨ u=1 1−x 0≤x≤1 u (x, 0) = ⎩ uR = 0 x≥1 cases (e) and (f) with different values for
. Here we have used x = 0.1 and t = 0.001.
4.5. Example 5 Here we want to apply proposed splitting methods on linear system of first order hyperbolic equations in the form
Ut +AUx = 0,
(4.7)
where U(x, t) is a vector in Rn and A is a matrix with real eigenvalues and a full set of eigenvectors. For example consider the system
u v
+ t
a11 a21
a12 a22
u v
= 0,
we can rewrite (4.8) in the form
u v
=− t
a11 0
0 a22
(4.8)
x
u v
− x
0 a12 0 0
u v
− x
0 a21
0 0
u v
x
and solve any of this subsystems by using the splitting ideas. Two hyperbolic equations corresponding to first subsystem
ut = −a11 u x v t = −a22 v x
can be solved by proposed splitting methods introduced in Section 3. Now consider the equations
ut = −a12 v x vt = 0
corresponding to the second subsystem. Equation v t = 0 implies that v is constant during a time step and the second equation
ut = −a12 v x can be solved by using finite difference methods like
u ∗∗ j =
−a12 t ∗ v j +1 − v ∗j −1 + u ∗j , 2 x
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
375
Fig. 6. Numerical solution of upwind type iterative splitting method on viscid Burgers’ equation with Gaussian initial function and = 0.1 (a), Gaussian initial function and = 0.01 (b), Riemann initial function and = 0.1 (c), Riemann initial function and = 0.01 (d), piecewise continuous initial function and = 0.1 (e) and piecewise continuous initial function and = 0.01 (f) with t = 0.001 and x = 0.1.
376
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Fig. 7. Numerical solution of system (4.8) with proposed splitting method with coefficients matrix (4.9).
where u ∗ and v ∗ are the solutions obtained from solving the first subsystem. The third system can also be solved in the same way. Let x = 0.01, t = 0.001 and
A=
−1 1 1
0
.
(4.9)
Fig. 7 shows the solutions u and v of the system (4.8) in domain (−1, 1) × (−1, 1) with proposed method and with the initial conditions
u (x, 0) = exp(−80(x + 0.33)2 ) v (x, 0) = exp(−80(x − 0.66)2 ), and homogenous Dirichlet boundary conditions. On the other hand Fig. 8 shows the solutions u and v of the system (4.8) with
A=
0 1 1 0
(4.10)
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
377
Fig. 8. Numerical solution of system (4.8) with proposed splitting method with coefficients matrix (4.10).
and with conditions and parameters same as before. It should be noted that this system is derived from the familiar second order wave equation w tt = w xx using the variables u = w t and v = − w x . 5. Conclusion In the current research, some numerical methods are obtained for solving hyperbolic partial differential equations. An upwind time splitting method with unconditionally maximum principle preserving property is presented, also damping of its solution is studied. By adding an artificial diffusion to the right hand of the advection equation a second order method with relatively mild condition for stability is derived. These results have important role in the numerical treatment of hyperbolic partial differential equations. As a nonlinear hyperbolic equation, Burgers’ equation is solved using the new ideas and an iterative procedure is implemented. Numerical examples validate the theoretical results and efficiency of the proposed methods. Splitting of the differential equation in certain parts and solving these parts exactly may be the main reason of the excellence of the present methods over the other existing numerical methods. The ides is used for solving two dimensional and system of hyperbolic partial differential equations.
378
R. Hosseini, M. Tatari / Applied Numerical Mathematics 146 (2019) 361–378
Acknowledgement The authors are very grateful to reviewers for carefully reading the paper and for their constructive comments and suggestions. This research was in part supported by a grant from IPM (No. 98650416). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
H. Bateman, Some recent researches on the motion of fluids, Mon. Weather Rev. 43 (1915) 163–170. G.P.E. Benton, A table of solutions of the one-dimensional Burgers equation, Q. Appl. Math. 30 (1972) 195–212. S. Blanes, F. Casas, A Concise Introduction to Geometric Numerical Integration, CRC Press, 2016. J.M. Burgers, A mathematical model illustrating the theory of turbulence, Adv. Appl. Mech. 1 (1948) 171–199. J. Cai, H. Liang, C. Zhang, Efficient high-order structure-preserving methods for the generalized Rosenau-type equation with power law nonlinearity, Commun. Nonlinear Sci. Numer. Simul. 59 (2018) 122–131. M. Caliari, A. Ostermann, C. Piazzola, A splitting approach for the magnetic Schrödinger equation, J. Comput. Appl. Math. 316 (2017) 74–85. S.A. Chin, Understanding Saul’yev-type unconditionally stable schemes from exponential splitting, Numer. Methods Partial Differ. Equ. 30 (6) (2014) 1961–1983. J.D. Cole, On a quasilinear parabolic equation occurring in aerodynamics, Q. Appl. Math. 9 (1951) 225–236. M. Dehghan, A. Taleei, A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients, Comput. Phys. Commun. 181 (2010) 43–51. J. Douglas Jr., H.H. Rachford Jr., On the numerical solution of the heat conduction problems in two and three variables, Trans. Am. Math. Soc. 82 (1956) 421–439. S.A. Dyachenko, A. Zlotnik, A.O. Korotkevich, M. Chertkov, Operator splitting method for simulation of dynamic flows in natural gas pipeline networks, Phys. D: Nonlinear Phenom. 361 (2017) 1–11. ´ E.G. Dyakonov, Difference schemes of second order accuracy with a splitting operator for parabolic equations without mixed partial derivatives, Ž. Vyˇcisl. Mat. Mat. Fiz., Mosc. 4 (1964) 935–941. J. Geiser, Additive and iterative splitting methods for multiscale and multiphase coupled problems, J. Coupled Syst. Multiscale Dyn. 4 (4) (2016) 271–291. J. Geiser, Iterative Splitting Methods for Differential Equations, CRC Press, 2017. J. Geiser, J.L. Huesob, E. Martinez, New versions of iterative splitting methods for the momentum equation, J. Comput. Appl. Math. 309 (2017) 359–370. P. González de Alaiza Martínez, M.E. Vázquez-Cendón, Operator-Splitting on Hyperbolic Balance Laws, Springer International Publishing, Cham, 2014, pp. 279–287. E. Hairer, Ch. Lubich, G. Wanner, Geometric Numerical Integration, second edition, Springer Series in Computational Mathematics, 2000. E. Hopf, The partial differential equation ut + uu x = μu xx , Commun. Pure Appl. Math. 3 (1950) 201–230. R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1990. H.G. Lee, A second-order operator splitting Fourier spectral method for fractional-in-space reaction-diffusion equations, J. Comput. Phys. 333 (2018) 345–403. S. MacNamara, G. Strang, Operator splitting, in: Splitting Methods in Communication, Imaging, Science and Engineering, 2016, pp. 95–114. G. Marchuk, Some applications of splitting-up methods to the solution of mathematical physics problems, Apl. Mat. 13 (1968) 103–132. K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, second edition, Cambridge University Press, 2005. D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 28–41. G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517. TheMathWorks Inc. (2017). H.F. Trotter, On the product of semi-groups of operators, Proc. Am. Math. Soc. 10 (1959) 545–551. G. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, 1974. J. Xie, Z. Zhang, D. Liang, A conservative splitting difference scheme for the fractional-in-space Boussinesq equation, Appl. Numer. Math. 143 (2019) 61–74, https://doi.org/10.1016/j.apnum.2019.03.013.