KRYLOV IMPLICIT INTEGRATION FACTOR WENO METHOD FOR SIR MODEL WITH DIRECTED DIFFUSION

SIR models with directed diffusions are important in describing the population movement. However, efficient numerical simulations of such systems of fully nonlinear second order partial differential equations (PDEs) are challenging. They are often mixed type PDEs with ill-posed or degenerate components. The solutions may develop singularities along with time evolution. Stiffness due to nonlinear diffusions in the system gives strict constraints in time step sizes for numerical methods. In this paper, we design efficient Krylov implicit integration factor (IIF) Weighted Essentially Non-Oscillatory (WENO) method to solve SIR models with directed diffusions. Numerical experiments are performed to show the good accuracy and stability of the method. Singularities in the solutions are resolved stably and sharply by the WENO approximations in the scheme. Unlike a usual implicit method for solving stiff nonlinear PDEs, the Krylov IIF WENO method avoids solving large coupled nonlinear algebraic systems at every time step. Large time step size computations are achieved for solving the fully nonlinear second-order PDEs, namely, the time step size is proportional to the spatial grid size as that for solving a pure hyperbolic PDE. Two biologically interesting cases are simulated by the developed scheme to study the finite-time blow-up time and location or discontinuity locations in the solution of the SIR model.

1. Introduction.There have been a lot of mathematical models studying the diffusion of biological populations since 1970s.These models can roughly be classified into two groups: One is based on the assumption that dispersal is due to a random motion of individuals [27,24,1]; the other is based on the assumption that dispersal is a response of population pressure, which is often called the directed spatial diffusion [9,22,15,5,4,3,16,18].In the biological world, some species have been observed migrating to avoid crowding rather than random motion.As early as recorded in 1971, arctic ground squirrels migrate from densely populated areas into sparsely populated areas, even when the latter provide a less favorable habitat [6,9].In 1983, MacCamy [15] proposed a model with directed spatial diffusion to avoid crowd in the setting of studying epidemic diseases.In this model, it is assumed that infectives and removed do not move while susceptibles move away from concentrations of infectives.Milner et al. [16,17,18] extended the model to a more general case: not only do each individual move away from crowd, but also susceptibles moves away from concentrations of infectives.Let S( x, t), I( x, t), and R( x, t) denote the density functions of susceptible, infected, and recovered human population at location x and time t.The model in [18] is formulated as for x ∈ Ω ⊂ R n and t ∈ (0, T ) for some T > 0, where P ( x, t) = S( x, t) + I( x, t) + R( x, t) is the density function of total population; α > 0 and γ > 0 are epidemiological parameters and k 1 ≥ 0 and k 2 ≥ 0 are diffusion parameters.The system is completed with nonnegative initial condition S( x, 0) = S 0 ( x), I( x, 0) = I 0 ( x), R( x, 0) = R 0 ( x) and no-flux boundary conditions The density-dependent cross-diffusion term ∇•(S∇I) has been used to model disease avoidence for studying pattern formation [2,28].
The mathematical models with random diffusion can be classified into diffusion partial differential equations hence the models are well-posed [27,1].However, most models with cross-diffusion are degenerate, sometimes even ill-posed [3].For the model 1, if we write the right-hand side of Eq 1 as a matrix times the column vector of the higher-order terms plus a term involving derivatives of lower order, then the three eigenvalues of the matrix are If k 1 > 0 and k 2 = 0, then we have two zero eigenvalues and one positive eigenvalue; If k 1 = 0 and k 2 > 0, then we have three zero eigenvalues; If k 1 > 0 and k 2 > 0, then we have one positive, one negative, and one zero eigenvalue.Li and Yip proved that the Eq. 1 is ill-posed when k 1 > 0 and k 2 > 0 in terms of long-term behavior [12].The well-posedness of Eq. 1 is still an open problem.But we conjecture that the Eq. 1 has a unique solution for some time T > 0. There are numerous efficient numerical methods for solving some mathematical models with spatial cross-diffusions [2,25,7].However, the degeneracy of Eq. 1 poses a challenge in simulating it.Zhao and Milner [18] proposed a split Runge-Kutta discontinuous Galerkin method to solve Eq. 1.However, the method has only first-order accuracy.In the special case k 1 = 0 and k 2 > 0, Eq. 1 might have a finite-time blow-up solution or there might develop a shock solution [17,18].
High order WENO schemes are a class of efficient numerical methods to stably resolve singularities while maintaining high order accuracy in smooth regions of the solution [23,30].It is natural to apply WENO schemes in solving this challenging model here.To perform efficient numerical simulations for this time dependent PDE system, an effective temporal numerical scheme is needed.The difficulty here is due to the stiffness of the system caused by nonlinear diffusion terms, since traditional implicit methods for solving stiff systems need to solve a large coupled nonlinear algebraic system at every time step, and often advanced preconditioning techniques have to be developed to improve the convergence of nonlinear algebraic solver.To avoid solving large nonlinear algebraic system at every time step and achieve large time step sizes, an efficient approach is the implicit integration factor (IIF) methods.IIF methods are a class of integration factor type time discretization methods, which perform the time evolution of stiff operators via evaluation of exponential functions of the corresponding matrices.IIF methods were originally developed in [19] for solving stiff one dimensional reaction-diffusion systems.A nice property of the methods is that the implicit terms are free of the exponential operation.The size of the nonlinear algebraic system resulting from the implicit treatment is independent of the number of spatial grid points.It only depends on the number of the original PDEs in the mathematical model itself, and no large nonlinear algebraic system needs to be solved at every time step.To overcome the difficulty in applying integration factor type methods to high dimensional problems, two different approaches were developed, i.e., the compact IIF methods [20,26] and the Krylov IIF methods [8,13].In [11], high order Krylov IIF WENO schemes were developed to efficiently solve fully nonlinear stiff advection-diffusion-reaction equations.
In this paper, we design the Krylov IIF WENO method to solve the SIR model Eq. 1.We focus on the two spatial dimension case in this paper.Denote X = [S, I, R] T , then Eq. 1 can be written in a matrix form where By forming the system Eq. 2 from the original model Eq. 1, the non-diffusion term H( X x , X y ) is separated from diffusion terms.It has a form of Hamilton-Jacobi type operators, so high order WENO schemes for Hamilton-Jacobi equations (e.g., [21,29]) can be naturally used to discretize H( X x , X y ).Then the Krylov IIF method for fully nonlinear stiff advection-diffusion-reaction equations in [11] can be smoothly adopted to solve the Hamilton-Jacobi-diffusion-reaction system Eq. 2.
The organization of the paper is as follows.In Section 2 we will introduce the Krylov IIF WENO numerical scheme for the SIR model Eq.2; In Section 3, we will demonstrate the order of accuracy of the proposed numerical method using extended systems with known analytic solutions; In Section 4, we will numerically solve the SIR model Eq. 2 for biologically interesting cases using the proposed numerical scheme; In Section 5, we will give some remarks and discussions.
2. The Krylov IIF WENO scheme for SIR model.
For discretizing the Hamilton-Jacobi term, we use the third order WENO scheme with Lax-Friedrichs flux [21,31], i.e., ), where for s = x, y.Here ρ(•) is the spectral radius operator.It is worth pointing out that ∂ H ∂ Xx and ∂ H ∂ Xy are 3 by 3 matrices, so the eigenvalues of them can be easily found analytically and spectral radius ρ can be computed efficiently.( X x ) ± i,j = [(S x ) ± i,j , (I x ) ± i,j , (R x ) ± i,j ] T and for u = S, I, R, a third order WENO approximation for u x at the grid point (i, j) when the wind "blows" from the left to the right is where + (u i+1,j − 2u i,j + u i−1,j ) 2 ; on the other hand, the third order WENO approximation for u x at the grid point (i, j) when the wind "blows" from the right to the left is where Similarly for approximations of derivatives in y-direction ] T , third order WENO approximations are computed as following: where where After the spatial discretizations, we consider the time discretization.Let the vector T denote the unknowns.For the boundary-point values, V i,0 , V = S, I, R, i = 1, 2, . . ., M − 1, we construct a cubic interpolation polynomial P 3 (y) of V at the grid points (x i , y 0 ), (x i , y 1 ), (x i , y 2 ), and (x i , y 3 ), and impose the no-flux boundary condition P 3 (y 0 ) = 0, which leads to To discretize the Hamilton-Jacobi term using the third order WENO scheme, we need ghost point value V i,−1 .The extrapolation technique is used.Namely, we construct a cubic interpolation polynomial of V at (x i , y 0 ), (x i , y 1 ), (x i , y 2 ) and (x i , y 3 ), and evaluate it at (x i , y −1 ) with V i,0 be replaced with the above equality, then we get Similarly, we obtain the values of S, I, R at the other boundary points.Now we apply the second-order IIF (IIF2) temporal discretization.We use the approach in [11] to deal with the nonlinear diffusion terms in the SIR model.Spatial discretizations of Eq. 1 leads to the following semi-discretized ordinary differential equation (ODE) system where H( U ) and F d ( U ) result from the spatial discretizations of the Hamilton-Jacobi term and the nonlinear diffusion terms respectively, and where E( U ) is the remainder.The Eq. 4 can be rewritten as Multiplying by e −C( U k )t and integrating Eq. 6 from t = t k to t = t k+1 , we get where The IIF schemes are obtained by approximating the integrand in Eq. 7 by Lagrange interpolation polynomials.See [11] for details.The second-order IIF scheme takes the following form: where . Rewrite the second order scheme Eq. 8, we get (9) Denoting the right hand side terms as RHS, at each grid point (x i , y j ) we have the following 3-equation nonlinear system Here (RHS) i,j,S corresponds to the component for S i,j in the right hand side vector of the Eq. 9, similarly for (RHS) i,j,I and (RHS) i,j,R .It is worth to note that this is just a local nonlinear system with 3 equations for S k+1 i,j , I k+1 i,j , and R k+1 i,j , for each pair of (i, j), 1 ≤ i ≤ M − 1, 1 ≤ j ≤ N − 1. Newton's method can be easily employed to solve the small nonlinear system.This shows an advantage of IIF schemes over many implicit schemes, namely, solving large coupled nonlinear algebraic system is avoided.
Evaluating of e C k ∆t k v is approximated by the Krylov subspace method for efficiently implementing IIF schemes to solve high spatial dimension problems [8].The sparse matrix C k is projected to the Krylov subspace The well-known Arnoldi algorithm generates an orthonormal basis where γ = v 2 , and e 1 denotes the first column of the m × m identity matrix I m .
In the paper, m = 25 is being used.
3. Order of convergence of the Krylov IIF WENO scheme.
3.1.One-dimensional cases.First we apply the proposed Krylov IIF WENO method to the following problem: where Ω = (0, π), Ω T = Ω × (0, T ), f S , f I , and f R are chosen so that the system has a solution The initial condition of the system is set by using the true solution with t = 0.The system is subject to no-flux boundary conditions.Throughout this paper, we take α = 1 and γ = 5.We set T = 1 for all the simulations in this section.One advantage of the Krylov IIF WENO scheme for solving the fully nonlinear Hamilton-Jacobidiffusion-reaction system is that the time step size ∆t k can be adaptively chosen as that for a pure hyperbolic PDE.The time step size ∆t k satisfies the following Courant-Friedrichs-Lewy (CFL) condition for the two-dimensional problem where CFL is a specified number, called the CFL number.α x and α y are defined in the Eq. 3. At every time step, α x and α y are updated using the Eq. 3, then the time step size ∆t k is determined by the above CFL condition.∆t k will be used to evolve the numerical solution to the next time step.In one-dimensional case, the above equality is reduced into ∆t k α x 1 ∆x = CFL.This indicates the high efficiency of the Krylov IIF WENO scheme for solving the problems with diffusion terms, i.e., large time step size computations (∆t = O(∆x)) are achieved.
3.1.1.Case 1: k 1 = 0.1 and k 2 = 0.001.As shown in [12], the Eq. 10 is ill-posed for its long-term behavior.To start the computation by the Krylov IIF2 WENO scheme, the numerical values at the first and the second time steps, i.e., U 0 and U 1 are needed.U 0 can be directly computed by using the initial condition of the PDE.Here for simplicity, U 1 is set by using the exact solution.For the problems with unknown exact solutions, U 1 is computed by using a second order Runge-Kutta method.The numerical solutions are obtained for Eq. 10 at T = 1.The CFL number is taken to be 0.2.Numerical results are reported in Table 1.From Table 1, it is clear that the desired second-order accuracy is obtained for the proposed Krylov IIF2 WENO scheme.It is worthy to note that if U 1 is computed with a second order Runge-Kutta method, a similar table with the same order of accuracy can be acquired.1. Numerical results of the one-dimensional system Eq. 10 for k 1 = 0.1 and k 2 = 0.001.π/N is the mesh size in the spatial direction.Here the constant CFL = 0.2.
3.1.2.Case 2: k 1 = 0.1 and k 2 = 0.In this case, the coefficient matrix of the second-order terms at the right hand side of Eq. 10 has two zero eigenvalues and one positive eigenvalue.The nonlinear term of Eq. 10 behaves like a porous media equation.We expect that for this well-posed case, the CFL number can take larger values than in the previous case.In addition, Table 2 shows the numerical results for CFL = 0.2 and Table 3 shows the numerical results for CFL = 0.5.It can be seen that both tables demonstrate the desired second order accuracy, while a smaller CFL number results in smaller numerical errors.Again, the time step size ∆t = O(∆x) shows that the Krylov IIF WENO scheme is an efficient numerical method for this type of problems which involve diffusion terms.
N L ∞ error L ∞ order L 1 error L  3. Numerical results of the one-dimensional system Eq. 10 for k 1 = 0.1 and k 2 = 0. CFL = 0.5.π/N is the mesh size in the spatial direction.
3.1.3.Case 3: k 1 = 0, k 2 = 0.1.When k 1 = 0, and k 2 > 0, the last two equations are reduced into ordinary differential equations.The proposed Krylov IIF2 WENO scheme can be still directly applied to this degenerate case.Table 4 shows the numerical results for this special case.The CFL number is taken to be 0.  4. Numerical results of the one-dimensional system Eq. 10 for k 1 = 0 and k 2 = 0.1.CFL = 0.1.π/N is the mesh size in the spatial direction.
3.2.Two-dimensional cases.In this section, we apply the Krylov IIF WENO method to the following two-dimensional problem: (11) This system is defined on the spatial domain (0, π) × (0, π).f S , f I , and f R are chosen so that the system has a solution S(x, y, t) = e −t (1 − cos(2x))(1 − cos(2y)), I(x, y, t) = e −t (1 + cos(x))(1 + cos(y)), R(x, y, t) = e −t (1 − cos(x))(1 − cos(y)).The initial condition of the system is set by using the true solution with t = 0.The system is subject to no-flux boundary conditions.The parameters α = 1 and γ = 5.Again for simplicity, U 1 is set by using the exact solution.Of course, a second-order Runge-Kutta method can also be used to compute U 1 to start the time evolution.We run the simulation till time T = 1 by the Krylov IIF2 WENO scheme.Table 5 shows the numerical results for the case where k 1 = 0.1 and k 2 = 0.001.The numerical results for the case where k 1 = 0.1 and k 2 = 0 are reported in Table 6, and those for the case where k 1 = 0 and k 2 = 0.1 are given in Table 7.We draw the same conclusion as that for the one-dimensional system.The proposed Krylov IIF2 WENO method reaches the desired second-order accuracy for solving the twodimensional system Eq 11.Large time step size ∆t = O(∆x) is achieved by using the efficient Krylov IIF2 WENO scheme to solve this system.
Remark.We have demonstrated that the proposed Krylov IIF2 WENO method has desired second-order accuracy for solving the model systems Eq. 10 and Eq.11.The method works very well for different cases where k 1 and k 2 take various values.In our numerical experiments, we also observe that the numerical method may suffer  7. Numerical results of the two-dimensional system Eq.11 for k 1 = 0 and k 2 = 0.1.CFL = 0.2.π/N is the mesh size in each of the spatial directions.
instability for some strongly ill-posed cases such as k 1 = 0.1 and k 2 = 0.1.In these strongly ill-posed cases, the coefficient matrix of the second order term has one relatively large negative eigenvalue, which gives a strongly backward heat equation component in the system.This may cause instability for the numerical method, especially in refined meshes.How to deal with these most difficult cases is one of our next research topics.
4. Application of the Krylov IIF WENO method to the SIR model.In this section, we apply the Krylov IIF2 WENO method to the SIR model with directed diffusion, i.e., the system Eq. 1. Two biologically interesting cases are studied: (i) k 1 = 0, k 2 > 0, corresponding to the directed dispersion avoiding infection; (ii) k 1 > 0, k 2 = 0, corresponding to the directed dispersion avoiding crowding.
4.1.Case 1: k 1 = 0, k 2 > 0 (avoiding infection).First let us consider the same example as in the Milner and Zhao's paper [18] for the one-dimensional case where k 2 = 0.1, α = 1, and γ = 5.Let Ω = (0, 1).The initial condition are chosen as S 0 (x) = 6, R 0 (x) = 0, for x ∈ Ω and The solution might break down with a formation of a shock wave, or a finite-time blow-up solution might exist by rewriting the system as a hyperbolic system [16].
The discontinuity occurs on the boundary of compact support of infected individuals, x = 0.8.Numerical simulations in [18] suggest that there might exist a finite-time blow-up solution at t = 0.7.Following the approach in [18], we compute the numerical solutions of S at t = 0.7 and at mesh points x = 0.8 − 1 M , M = 1280, 640, . . ., 10.Then, we compute the ratio r = S(0.8− 1 M )/S(0.8− 1 M/2 ).If r approaches a constant, we can conjecture that S takes the form S(x) = const (0.8 − x) α where α = ln(r)/ln (2).Table 8 shows the numerical results using the Krylov IIF WENO scheme for the system with the above initial condition.We can see that the value of r still has a little oscillation at t = 0.7 when M approaches 1280.The convergence of r has a better pattern at t = 0.706, as shown in Table 8.This indicates that the solution might blow-up at t = 0.706 instead of t = 0.7. Figure 1 shows the density profile of solutions of Eq. 1 at t = 0.706.We can see that the singularity of the solution has been resolved well by the Krylov IIF WENO 8. Numerical solution of S at t = 0.7 and t = 0.706 and x = 0.8 − 1/M for different M .All values are computed using the uniform mesh ∆x = 1/1280.Numerically verifying the finite-time blow-up is challenging.Hirota et al proposed a numerical method to estimate the blow-up time for the problems of partial differential equations [10].However, the method was applied to the test problems where the partial differential equation does have finite-time blow-up solutions.We conjecture that susceptibles might be finite-time blow-up at x = 0.8, the compact support boundary of the infective individuals, at time t = 0.706.On the other hand, if we do not assume that there is a finite-time blow-up solution, we will expect that a form of shock or discontinuity of solutions occurs.As infected individuals can not move away from their compact support region, susceptibles will stop moving further as soon as they reach a safe zone of infection free region.This will result to a situation where the population accumulation of susceptibles on the boundary of compact support of infected individuals.Numerical simulation shows the discontinuity of solution, particularly in I. one-dimensional case of Eq. 1 at t = 1.5.We can see that the discontinuities of the solution are resolved very well, even with a relatively coarse mesh N = 40.The WENO approximation plays a key role here in capturing the discontinuities of the solution sharply and stably.In the two-dimensional case, Milner and Zhao [18] also conjectured that a finitetime blow-up solution might exist.Let Ω = (0, 1) × (0, 1), S 0 (x, y) = 6, R 0 (x, y) = 0 on Ω, and I(x, y) is the tensor product of the function in Eq. 12. Again, we take the parameter values k 1 = 0 and k 2 = 0.1.Figure 3 shows the numerical solution of Eq. 1 for the two-dimensional case with the above initial conditions.From the figure, it is reasonable to believe that the finite-time blow-up solution occurs first at the intersection of domain boundary and the boundary of compact support of initial conditions of infected individuals.As that for the one-dimensional case, the singularities of the solution have been resolved very well by the Krylov IIF WENO scheme.In [18], it was observed that the finite-time blow-up might occur along the boundary of the initial compact of infected individuals.We believe that the higher-order numerical scheme in this paper can better capture the behavior of finite-time blow-up of solutions.Similarly as in the one-dimensional case, the numerical scheme in this paper allows us to run the simulation beyond t = 0.7.We observe accumulation of susceptibles on the boundary of compact support of initial infected population as t increases.
It is clear that the Krylov IIF WENO scheme can be directly applied to the Eq. 13 if we are only interested in the dynamics of total population.Here we still solve the original system Eq. 1 to find S, I and R by the Krylov IIF WENO method.The total population P is evaluated by adding S, I and R we obtain.
In the one-dimensional case, we use the same initial condition as that in the previous section.The parameter values are set as k 1 = 0.1 and k 2 = 0.At the steady-state, the infected individuals will vanish and the total population will approach a constant.This is verified by our numerical simulation.Figure 4 shows the density profiles of populations at t = 0 and at t = 20, respectively.At t = 20, the numerical solution reaches the steady-state.In this case, the CFL number can be taken as large as 0.6.In the two-dimensional case, we use the tensor product of the initial condition in the one-dimensional case.We take CFL = 0.2. Figure 5 shows the initial density of populations.The numerical solutions of Eq. 1 at t = 1, t = 10 and t = 25 are shown in Figure 6, Figure 7 and Figure 8, respectively.From the figures, we see that the total population P approaches a constant; the infected individuals vanish; and the susceptible and recovered individuals converge to their own steady-states.This observation is consistent with that in the one-dimensional case.

Conclusion and discussions.
Mathematical models with directed diffusions are important in describing the population movement.In particular, the population will move away from crowd and also away from infection during epidemics.However, such models are usually difficult to simulate since the describing system of partial differential equations is often degenerated or sometimes even not well-posed.For example, the leading coefficient matrix of Eq. 1 might have one positive eigenvalue, one negative eigenvalue and one zero eigenvalue when k 1 > 0 and k 2 > 0. Therefore, the system might behave as a mixture of heat equation, backward heat equation, and hyperbolic system.Numerically simulating such models is often very challenging.
In this paper, we design a Krylov implicit integration factor WENO method to solve the SIR models with directed diffusions.The method is demonstrated to have second-order accuracy by using test problems with smooth solutions, for the general case (k 1 > 0 and k 2 > 0) and two biologically interesting cases: (i) avoid infection (k 1 = 0, k 2 > 0) and (ii) avoid crowd (k 1 > 0 and k 2 = 0).Numerical simulations show the high efficiency of the Krylov IIF WENO method in solving the SIR models with directed diffusion.Singularities in the solution are resolved stably and sharply by the WENO approximations in the scheme.Unlike a usual implicit method for solving stiff nonlinear PDEs, the Krylov IIF WENO method avoids solving large coupled nonlinear algebraic systems at every time step.Only a system of three nonlinear equations need to be solved at every grid point for the SIR model in this paper.These small nonlinear systems can be solved easily by the Newton's method.Large time step size computations are achieved, namely, the time step size ∆t = O(∆x) for solving the fully nonlinear second-order PDEs in this paper.The time step size ∆t is adaptively determined at every time step by the CFL condition as that for a pure hyperbolic PDE.
The numerical results of the Krylov IIF WENO method applied to Eq. 1 show the consistent results as that obtained in [18].Furthermore, it suggests that the finite-time blow-up might occur on the special points in the spatial domain rather than along a curve for the two-dimensional case.High order numerical methods are important in capturing the details of the system, such as finding the finite-time blow-up time and location or discontinuity locations in this case.
Mathematical modeling population's spatial movement is biologically important.Other cross-diffusion models were studied in the literatures (see the introduction).High order numerical methods will be efficient to simulate such models.A possible future work is to apply the Krylov IIF WENO scheme to other cross-diffusion mathematical models.Human spatial movement typically happens in two dimensions.The proposed Krylov IIF WENO scheme can be easily applied to three dimensional applications, such as fish movement in the contaminated water system.The spatial discretizations are based on finite difference schemes, hence dimensionby-dimension approach can be directly applied to multidimensional problems (e.g., three dimensional problems).Although for three dimensional problems or even higher dimensional problems, we will obtain a ODE system with much larger size than that for two dimensional problems, the Krylov IIF scheme can still solve the system efficiently, as that shown in [14].For problems with other boundary conditions, corresponding algebraic equations can be derived for the boundary points based on that specific boundary condition.Then the derived algebraic equations will be used to form the ODEs for grid points adjacent to the boundary points.Hence we will obtain a slightly different ODE system, which can be solved directly by the proposed Krylov IIF scheme.

Figure 5 .
Figure 5.Initial density profiles of population at t = 0.

Table 2 .
Numerical results of the one-dimensional system Eq. 10 for k 1 = 0.1 and k 2 = 0. CFL = 0.2.π/N is the mesh size in the spatial direction.N L ∞ error L ∞ order L 1 error L 1 order L 2 error L

Table 5 .
Numerical results of the two-dimensional system Eq.11 for k 1 = 0.1 and k 2 = 0.001.CFL = 0.4.π/N is the mesh size in each of the spatial directions.N L ∞ error L ∞ order L 1 error L 1 order L 2 error L

Table 6 .
Numerical results of the two-dimensional system Eq.11 for k 1 = 0.1 and k 2 = 0. CFL = 0.6.π/N is the mesh size in each of the spatial directions.
N L ∞ error L ∞ order L 1 error L 1 order L 2 error L scheme.