Krylov SSP Integrating Factor Runge–Kutta WENO Methods

: Weighted essentially non-oscillatory (WENO) methods are especially efﬁcient for numerically solving nonlinear hyperbolic equations. In order to achieve strong stability and large time-steps, strong stability preserving (SSP) integrating factor (IF) methods were designed in the literature, but the methods there were only for one-dimensional (1D) problems that have a stiff linear component and a non-stiff nonlinear component. In this paper, we extend WENO methods with large time-stepping SSP integrating factor Runge–Kutta time discretization to solve general nonlinear two-dimensional (2D) problems by a splitting method. How to evaluate the matrix exponential operator efﬁciently is a tremendous challenge when we apply IF temporal discretization for PDEs on high spatial dimensions. In this work, the matrix exponential computation is approximated through the Krylov subspace projection method. Numerical examples are shown to demonstrate the accuracy and large time-step size of the present method.


Introduction
We consider the two-dimensional general form of the nonlinear hyperbolic partial differential equation (PDE) u t + f (u) x + g(u) y = 0, where u is the unknown function, f is the flux function in x spatial direction, and g is the flux function in y spatial direction. Nonlinear hyperbolic PDEs (1) are common mathematical models in applications from science and engineering. To numerically solve this time-dependent problem (1), we first need to choose an efficient nonlinearly stable spatial discretization. The solutions to hyperbolic PDEs may have discontinuities, sharp gradients, or discontinuous derivatives, and so forth. A straightforward high-order scheme for such kinds of problems will generate instability called the Gibbs phenomena. Thus, we desire schemes that not only have a uniformly high order of accuracy in smooth regions of the solution, but can also resolve singularities in an accurate and essentially nonoscillator (ENO) fashion. Weighted essentially nonoscillatory (WENO) schemes are a popular class of such schemes. WENO schemes [1,2] were designed based on the successful ENO schemes in [3][4][5] with additional advantages. For problems with physically complicated solution structures, such as turbulence systems in fluid dynamics, WENO methods with very highorder accuracy and refined computational grids can resolve the physical solution with mixed stable and unstable behavior [6]. The essence of the WENO schemes is to reconstruct a weighted combination of some local reconstructions based on different small stencils. The combination coefficients, which are also called nonlinear weights, are determined by the linear weights and the smoothness indicators. The linear weights are usually selected to increase the order of accuracy on each small stencil, and the smoothness indicators measure the smoothness of the reconstructed function in the relevant small stencils. Since they were developed, WENO schemes have been studied and applied extensively for solving problems with both singularities and a complex solution structure. Here are some examples. Fast sweeping WENO schemes and homotopy WENO schemes were designed in [7][8][9][10][11][12] to efficiently solve hyperbolic PDEs' steady-state problems. To solve stiff convectiondiffusion-reaction PDEs, high-order Krylov implicit integration factor WENO schemes were developed in [13][14][15][16]. In [17], sparse grid WENO methods were developed to solve high-dimensional hyperbolic PDEs.
Since WENO schemes were designed to handle problems with both complicated solution structures and discontinuities, they require more operations than many other schemes because of their sophisticated nonlinear properties and high-order accuracy. The computational cost rises greatly when the PDEs are on 2D or higher, or the number of grid points is large. Especially for solving long time simulations or solving steady-state problems, efficient large time-step sizes and high-order temporal discretization methods are crucial. Many innovatory high-order temporal schemes were designed, such as the total variation diminishing (TVD) Runge-Kutta (RK) schemes [4,[18][19][20], spectral deferrred correction (SDC) methods [21][22][23][24][25], high-order implicit-explicit (IMEX) multistep/RK methods [26][27][28][29][30][31][32][33] and their applications, hybrid methods of SDC and high-order RK schemes [34], SSP multistage two-derivative time-stepping schemes [35], and so forth. All of these numerical methods have their own advantages.
Challenges in designing large time-step size schemes for (1) also come from restrictions on local time-step sizes due to both linear stability and nonlinear stability. Since discontinuous solutions often arise in hyperbolic PDEs, the linear stability theory itself does not guarantee convergence. SSP temporal numerical schemes were developed [4,36] to assure the preservation of the nonlinear stability properties. That is, if the high-order spatial method (for example, a WENO scheme) coupled with forward Euler satisfies the nonlinear stability, then the same spatial discretization coupled with these higher-order time discretizaton methods will preserve the nonlinear stability. Exponential integrators such as integrating factor methods are time schemes which have large linear stability regions. They can achieve very large time-step sizes or Courant-Friedrichs-Lewy (CFL) numbers for problems with smooth solutions [13][14][15]37]. However, these schemes do not satisfy the SSP property. Hence, direct applications of these schemes to hyperbolic PDEs with discontinuous solutions may lead to nonlinear instability. Recently, integrating factor Runge-Kutta (IFRK) methods with a SSP property were developed in [38][39][40] to solve one-dimensional hyperbolic PDEs, which have explicit linear and nonlinear parts. These methods are exponential integrators which satisfy the SSP requirement. They can completely relax the time-step size restriction coming from the linear part of the high-order spatial discretization. However, the question of how to apply the method to multi-dimensional problems, which do not have explicit linear and nonlinear parts, is still open.
In this paper, we extend WENO methods with large time-stepping SSP integrating factor Runge-Kutta time discretization to solve general nonlinear multi-dimensional problems. When applying the integrating factor methods to multi-dimensional problems, evaluating the product of the matrix exponential and a vector efficiently is crucial and challenging. To overcome this difficulty, we utilize the Krylov subspace technique to approximately calculate such products. In [41], Gallopoulos and Saad used and analyzed the Krylov subspace methods to approximate the product of an exponential matrix and a vector. In [42,43], we applied the Krylov subspace method to implicit integration factor methods to solve high-dimensional reaction-diffusion PDEs. To solve hyperbolic PDEs, in this paper, we first use third-or fifth-order WENO schemes for the spatial discretization. In order to apply exponential integrators to the resulting nonlinear ODE system, we split it into linear and nonlinear parts. Then, for the time discretization, we combine the SSP IFRK methods [38] with the Krylov subspace approximation method to solve two-dimensional general nonlinear hyperbolic problems.
The organization of the paper is as follows. WENO methods with the Krylov SSP IFRK methods are derived and formulated in Section 2. The stability region for our schemes are analyzed in Section 3. In Section 4, we present some numerical experiments. Discussions and conclusions are given in Section 5.

Krylov SSP Integrating Factor Runge-Kutta Methods
Consider two-dimensional hyperbolic conservation laws (1) with appropriate boundary conditions. Since high-order finite difference WENO methods are very efficient in resolving complex solution structures and keeping nonlinear stability for multi-dimensional problems due to its dimension-by-dimension implementation fashion [6], it is natural to apply high-order finite difference WENO schemes in the spatial discretization of (1).

Spatial Discretization: The Third and Fifth-Order WENO Schemes
We use a finite difference scheme to approximate the values of the hyperbolic terms f (u) x + g(u) y at a grid (x i , y j ) in a conservative fashion. To be specific, we use a conserva- to approximate f (u) x at (x i , y j ) along the line y = y j . Here, In the WENO schemes, the following Lax-Friedrichs flux splitting is used. where f + (u) and f − (u) are the positive and negative wind parts, respectively. Note that, for simplicity, we omit the index j for the y direction when the approximations for the x direction are described.
Let's first look at the third-order WENO (WENO3) method formulation. The stencil (in the 2D case) is given on the left of Figure 1. The numerical fluxf + i+1/2 depends on a left-biased stencil, namely, the value f (u k ) at three notes k = i − 1, i, i + 1, because the wind is positive; that is, f (u) ≥ 0 for the scalar case. Its formula is given bŷ where Here, d 0 = 2 3 and d 1 = 1 3 are called the linear weights, and To prevent a zero denominator from happening, we choose as a tiny positive number. From the Formula (6), we can see , which are both second-order numerical fluxes and are built on two substencils, (x i , x i+1 ) and (x i−1 , x i ), respectively. The coefficients ω i are determined by β i in each substencil.
Forf − i+1/2 , this is the case when f (u) < 0, or a negative wind. A WENO3 scheme can be constructed by using a right-biased stencil with values f (u i ), f (u i+1 ), and f (u i+2 ). whereω For the fifth-order WENO scheme (WENO5), the stencil is shown on the right of Figure 1. The numerical fluxf + i+1/2 depends on a left-biased stencil, namely, the five-point value Its formula is given bŷ where with d 0 = 3 10 , d 1 = 3 5 , and d 2 = 1 10 , and From the Formula (10), we can see that this numerical fluxf + 1+1/2 is a convex combination of three brackets, and each one of them is a third-order numerical flux and is built on three different substencils of three points each. The coefficients ω i are determined by β i in each substencil.
For the case of negative wind, that is, f (u) < 0, a fifth-order WENO approximation tof − i+1/2 can be constructed by using a right-biased stencil with numerical values whereω r =α withd 0 = 1 10 ,d 1 = 3 5 , andd 2 = 3 10 , and Combining Equations (3), (6), and (8), we have the numerical flux for the third-order WENO method Combining Equations (3), (10), (12), and using the following notations we have the numerical flux for the fifth-order WENO method Here, W0 + i means that f + is used for f in (15), and W0 − i means that f − is used for f in (15), and so on.
We apply the similar procedures in the y direction to approximate the numerical flux of g(u) y , then we have a nonlinear system where n and m are the numbers of grid points in x and y spatial directions, respectively.

Temporal Scheme: SSP Integrating Factor Runge-Kutta Schemes
After the spacial discretizations, a semi-discretized nonlinear ODE system (17) is obtained. We write it in the vector form Here, U is the vector of the solution variables, and F( U) is the approximation for the nonlinear convection terms by high WENO schemes. Each component of F( U) is a nonlinear function of numerical solution values on approximation stencils of WENO schemes. In order to apply exponential integrators to the nonlinear system (18), we factor out the linear part of the system as was done in [13] for the fully nonlinear convection-diffusion PDEs. Specifically, at every time-step interval [t k , t k+1 ] with time-step size ∆t = t k+1 − t k , we perform the splitting: where C k U is the linear part with being the Jacobian matrix (or an approximation of the Jacobian matrix if its direct computation is difficult) and U k being the numerical solution of U(t k ). F N is the nonlinear remainder of the system (18) after the splitting For example, for 1D problems, if c ij is the i th row and j th column element of the based on the 1D version of (17).
For the third-order WENO method, from Equation (14), we have where ∂ω r ∂u j and ∂ω r ∂u j , r = 0, 1 can be easily obtained based on the definitions of (7) and (9).
For all j except j = i − 2, i − 1, i, i + 1, i + 2, and i + 3, we have Remark 1. C k is a n × n sparse matrix for a one-dimensional problem. For two-dimensional problems, the dimension-by-dimension procedure is followed. A nm × nm sparse matrix C k is obtained. However, we only store the nonzero elements and their locations in the matrix C k for the Krylov approximation procedure.
Then, for system (19), we apply the integration factor approach in [13,14], with e −C k t as the integration factor, and the resulting integration factor equations of the system is Denoting Y = e −C k t U, then (21) can be written in terms of Y as The next step is to apply Runge-Kutta methods to (22) to approximate Y(t); then the relation U = e C k t Y is used to obtain the exponential integrator time schemes for U. To obtain the SSP exponential integrator, we apply the method in [38]. As pointed out in [38], application of a SSP Runge-Kutta method to derive an exponential integrator time scheme is not enough to ensure the preservation of the SSP property. In order to obtain SSP exponential integrator schemes for U based on SSP Runge-Kutta methods, the system (22) must be solved by SSP Runge-Kutta methods with non-decreasing abscissas. For example, the third-order SSP Runge-Kutta method with non-decreasing abscissas is used to solve (22) and leads to the following third-order SSP exponential integrator:

Integrating Factor Methods Based on Krylov Subspace Approximation
The efficiency of the SSP IF Runge-Kutta methods such as (23) largely relies on the methods to compute the product of the matrix exponential and a vector, e C k ∆t U k . When we use a specific spatial discretization to solve a PDE, a large matrix C k such as (20) will be generated. Although C k is sparse, the exponential matrix e C k ∆t is not. For one-dimensional problems, e C k ∆t can be directly computed, since the size of the matrix C k is manageable. However, for 2D problems, as an illustration, if we use a finite difference method on a n × n rectangular mesh, then we will obtain matrices C k with size n 2 × n 2 . Therefore, it cannot be afforded to compute and store e C k ∆t directly with regard to both computational cost and computer storage, due to the huge size of the matrix C k . Fortunately, if we take a close look at (23), we will see that we do not need e C k ∆t itself-all we need are the products of e C k ∆t and some vectors in (23). An outstanding option for doing this job is the Krylov subspace approximation, which can provide both accuracy and efficiency. Now, we implement the Krylov subspace technique to find approximations to the products of some exponential matrices and vectors in the SSP IF Runge-Kutta scheme (23) and derive the new Krylov SSP IF Runge-Kutta methods. First, we review the procedure of the Krylov subspace methods to approximate e A∆t v, following the literature (e.g., [41,44]).
First, we form a Krylov subspace K M based on the large sparse matrix A and the vector v in the following way Here, M is the dimension of the subspace K M , and it is much smaller than the size of A. We denote V M = [v 1 , v 2 , v 3 , · · · , v M ] as an orthonormal basis of the Krylov subspace K M , which is generated by the Arnoldi algorithm [45] in the steps described below (Algorithm 1).
Step 2.4 If h j+1,j = 0, then Set M = j; Stop and exit; else The vectors h i,j produced from the above Arnoldi algorithm consist of an M × M upper-Hessenberg matrix, denoted by H M . Because the vectors in V M are orthogonal, it is easy to obtain The above equation shows that H M is the projection of A to the Krylov subspace K M , with respect to the basis V M . Note that the dimension of H M is M, which is much smaller than that of the matirx A. What's more, because V M is orthonormal, the vector V M V T M e A∆t v is the orthogonal projection of e A∆t v on K M , that is, V M V T M e A∆t v is the best approximation to e A∆t v from K M . Thus, we have where e 1 is the first column of the identity matrix I M with size M × M. From (25), we obtain the approximation e A∆t v ≈ v 2 V M e H M ∆t e 1 .
Therefore, the large e A∆t matrix exponential computation is substituted by a small e H M ∆t matrix exponential computation. The much smaller e H M ∆t problem can be calculated by a scaling and squaring algorithm with a Padé approximation. The computational cost of this algorithm is only O(M 2 ), see [41,44,46]. After implementing the Krylov subspace approximation (27) to (23), we obtain the Krylov SSP IF Runge-Kutta scheme.
We choose the value of M to be small, but still large enough to ensure that the truncation errors of the numerical methods (23) are much larger than the errors of Krylov subspace approximations. In this paper, we use M = 25. From our numerical experiments, we observe that the proposed Krylov SSP IF RK methods have already achieved clear accuracy orders for such a small value of M. What's more, we do not need to increase the size of M when we refine the spatial-temporal resolution.

Linear Stability Analysis
To discuss the linear stability of (23), we consider the test model problem where Lu is the linear part, and Nu is the linearization of the nonlinear part F N in (19). Based on applying the third-order SSP IF Runge-Kutta method (23) to the test problem (28), we plot the stability region boundaries, which are some curves for different values of ∆tL. This is done by following the similar analysis for ETD methods in [47,48]. Note that the linear stability analysis here is different from that in [38].
Applying (23) to the Equation (28), we obtain u k+1 = e ∆tL u k 59 128 (29) Let λ = N∆t, and substituting u k = e ikθ into the above equation, we have We plot the stability regions on the complex plane of λ in Figure 2, for ∆tL = −0.5, −5, and −10. The interior of the closed curves is the stable region. The size of the stability region is an increasing function of |∆tL|. As L → −∞, the stability region approaches the whole complex plane. This is consistent with the property of the method, because the linear part is integrated exactly. Therefore, the larger the linear part is, the larger the stability region is.

Numerical Experiments
In this section, numerical examples are presented to demonstrate the stability, accuracy, and large CFL numbers for the WENO spatial discretizations coupled with the third-order Krylov SSP integrating factor Runge-Kutta method for 1D and 2D scalar hyperbolic equations. Both the linear advection equation and nonlinear Burgers' equations are tested. For all numerical experiments in this paper, the dimension of the Krylov subspace M is taken as 25.

Example 1. Consider the one-dimensional Burgers' equation
with periodic boundary conditions. The initial condition is given by u(x, 0) = 0.3 + 0.7 sin πx. The exact solution can be obtained by using the method of characteristics and Newton's method, namely, where x 0 is the solution to the nonlinear equation The computation is carried up to time T = 0.5 π 2 , which is before shocks are formed. Here, CFL=α x ∆t ∆x , where ∆t and ∆x are the time-step size and spatial step size, respectively, and α x is defined in (5). Different CFL numbers are tried. For CFL=4 and CFL= 10, the L ∞ and L 1 errors and order of accuracy for the WENO3 method coupled with the proposed third-order Krylov SSP IF Runge-Kutta method are reported in Table 1. We can see that the CFL number can be as large as 10 in the new method while still keeping the third order of accuracy. When we implement the fifth-order WENO with the third-order SSP IF Runge-Kutta method, in order to keep the consistency of truncation errors in the spatial and temporal directions, we are using ∆t = (∆x) 5/3 . The results are reported in Table 2. We observe that the numerical accuracy orders are actually higher than five. We also notice that for n = 80, WENO5 has larger errors than WENO3. This is because WENO schemes are highly nonlinear, and they have asymptotical convergence properties. We also implement the proposed WENO schemes with SSP IF Runge-Kutta methods till time T = 1, when shocks are formed. In Figure 3, the exact solution and the numerical solutions of the third-order WENO (top) and the fifth-order WENO (middle) are displayed. It is clearly shown in the Figure that the numerical shock is well-resolved for both WENO methods. We also plot the errors in log scale for both WENO3 and WENO5 at T = 1 (bottom). We can see that WENO5 has much smaller errors than WENO3 in smooth regions at this time.

Example 2. Consider the two-dimensional linear advection equation
with periodic boundary conditions. The initial condition is given by u(x, y, 0) = sin(x + y). The exact solution is u(x, y, t) = sin(x + y − 2t). For two-dimensional cases, the proposed Krylov method is used. The CFL number in this case is defined by CFL= ∆t α x ∆x + α y ∆y , where α x = max u | f (u)|, and α y = max u |g (u)|.
The computation is carried up to T = 1. The L ∞ and L 1 errors and order of accuracy for the third-order WENO method coupled with the third-order Krylov SSP IF RK method are reported in Table 3 for different CFL numbers. We see that the CFL number can be as large as five in the proposed new method, while still keeping the third order of accuracy. Again, when we implement the fifth-order WENO with the third-order Krylov SSP Runge-Kutta method, in order to keep the consistency of truncation errors in the spatial and temporal directions, we are using ∆t = min((∆x) 5/3 , (∆y) 5/3 ).
The results are reported in Table 4. We observe that the numerical accuracy orders are actually higher than five.
with periodic boundary conditions. The initial condition is given by u(x, y, 0) = 0.3 + 0.7 sin(π(x + y)/2). The exact solution can be obtained by using the method of characteristics and Newton's method.
This problem has a smooth solution before time 0.5 π 2 ; after that, shocks will form. We first calculate up to T = 0.5 π 2 with M = 25 at which L ∞ and L 1 errors are measured. Numerical errors and order of accuracy for the third-order WENO with the third-order SSP IF Runge-Kutta method are recorded in Table 5. We see that the desired third order of accuracy is achieved with the CFL number as large as five. Again, when we implement the fifth-order WENO with the third-order SSP Runge-Kutta method, in order to keep the consistency of truncation errors in the spatial and temporal directions, we use ∆t = min((∆x) 5/3 , (∆y) 5/3 ). The results are reported in Table 6. We see that a fifth-order of accuracy is achieved. We also implement the Krylov SSP IF Runge-Kutta WENO methods till time T = 1, when shocks are formed. In Figure 4, the numerical solutions of the third-order WENO (top) and the fifth-order WENO (bottom) are displayed. We can see that the numerical shocks are resolved very well for both WENO methods.

Discussions and Conclusions
In this paper, we developed a class of efficient high-order numerical methods by extending the work in [38] to solve general nonlinear 2D hyperbolic equations. High-order WENO schemes are incorporated with a large time-stepping third-order SSP integrating factor Runge-Kutta method. The reason we chose SSP RK methods as the time discretization is because they are the most efficient time-method with WENO spatial discretization for solving pure hyperbolic PDEs. A detailed discussion can be found in [20]. It will also be interesting to compare the method proposed here with other methods, such as semi-implicit schemes, which will be our future work.
Due to the large size of matrix exponentials rising from two-dimensional problems, we used the Krylov subspace projection method to efficiently approximate the matrix exponential operator. Numerical examples show that we can use large CFL numbers to achieve the desired accuracy of the methods. Numerical examples also verify that the error of the Krylov subspace approximation does not affect the accuracy orders of the third-order integrating factor SSP Runge-Kutta time discretization or the high-order WENO spatial discretizations, which suggests that the Krylov subspace approximation error is already much smaller than the temporal and spatial truncation errors of the numerical schemes.
The linear stability analysis in the paper is only for the time scheme. Because WENO schemes are highly nonlinear schemes, the complete stability analysis for the fully discretized schemes, including both time and spatial directions, is still an open problem, which will be studied in future research.
In this first paper, we only tested the proposed method for scalar equations. In future work, we will extend the method to solve more complicated hyperbolic PDE systems and higher-dimensional problems.