A new compact finite difference quasilinearization method for nonlinear evolution partial differential equations

Abstract This article presents a new method of solving partial differential equations. The method is an improvement of the previously reported compact finite difference quasilinearization method (CFDQLM) which is a combination of compact finite difference schemes and quasilinearization techniques. Previous applications of compact finite difference (FD) schemes when solving parabolic partial differential equations has been solely on discretizing the spatial variables and another numerical technique used to discretize temporal variables. In this work we attempt, for the first time, to use the compact FD schemes in both space and time. This ensures that the rich benefits of the compact FD schemes are carried over to the time variable as well, which improves the overall accuracy of the method. The proposed method is tested on four nonlinear evolution equations. The method produced highly accurate results which are portrayed in tables and graphs.


Introduction
Many researchers are now using higher order compact finite difference schemes in place of the conventional second order finite difference to solve differential equations [1][2][3][4][5]. This is because significant improvements to the accuracy of numerical solutions have been obtained by using fourth or sixth order compact FD schemes. Another advantage is that the high accuracy is obtained on coarser grids which ensures greater computational efficiency. Lele [6] presented various compact finite difference schemes for applications such as evaluating high order derivatives, interpolation and filtering.
A number of researchers have used compact FD schemes to solve partial differential equations arising in various fields. In all the previous work the compact FD schemes were only applied on space variables and other discretization techniques are used for the time variables. Quite a number of numerical techniques have been paired with the compact FD schemes in order to discretize time. For example [7] used a two-step predictor-corrector algorithm called McCormack method to find the solution of Burgers equation. Li and Visbal [3] used the classical fourth-order four-stage Runge-Kutta method (RK4). A lot of of researchers have combined the compact FD schemes with the third order total variation diminishing Runge-Kutta (TVD-RK3) scheme (see [4,[8][9][10]). Dlamini et al [11] combined the compact FD schemes with the Crank-Nicolson technique in time to solve unsteady boundary layer flow problems.
However these techniques limit the accuracy of the difference scheme to fourth order or less in time which means small mesh sizes have to be used in order to get desirable accuracy, and thus much more computational work is involved. The aim is to try the compact FD schemes (sixth order) on the time variable as well which will allow us to archive greater accuracy with a rough step size.
The purpose of this investigation is to explore the possibility of applying sixth order compact FD schemes in both space and time variables. We consider this on a few examples of nonlinear evolution PDEs. These are important equations arising in a number of fields of science and engineering. They are used to describe many complex nonlinear settings in applications such as vibration and wave propagation, fluid mechanics, plasma physics, quantum mechanics, nonlinear optics, solid state physics, chemical kinematics, physical chemistry, population dynamics, and many other areas of mathematical modelling. Understanding the behavior of these problems requires finding solutions of the differential equations and analyzing them. Most of the evolution equations are very complex to solve analytically and so we rely on approximate solutions. So the development of numerical solutions to solve such problems continues to be an active area of research.
Numerical experiments show that implementing the compact FD schemes on both space and time yield highly accurate results. This was experimented on the Fisher's, Burgers-Fisher, Burgers-Huxley and Fitzhugh-Nagumo equations. Before using the compact FD schemes we first linearize the PDEs using a quasilinearization technique developed by Bellman and Kalaba. The method is therefore called the compact finite difference quasilinearization method (CFDQLM). The rest of the paper can be summarized as follows, in section 2 we give a description of how to use the sixth order compact finite differences to discretize in space and in time as well as giving a description of how the CFDQLM is used to solve evolution partial differential equations. In section 3 we outline the examples considered in this work. We then present and discuss results in section 4 and in the last section we conclude.

Compact finite difference schemes
In this section we introduce the compact finite difference schemes and illustrate how we use them to solve nonlinear evolution partial differential equations. We consider nonlinear PDEs of the form @u @t D G u; @u @x ; defined on the the region t 2 OE0; T ; x 2 OEa; b with boundary condition u.a; t / D u a ; u.b; t / D u b ; In the derivation of the compact FD schemes we consider the function u.x; t / that depends on space variable x and temporal variable t . We consider one-dimensional uniform meshes on the regions OEa; b and OE0; T , with nodes x i (i D 1; 2; : : : ; N x ) and t j (j D 1; 2; : : : ; N t ) respectively, where and The distance between any two successive nodes is a constant x D x i x i 1 for the space variable and t D t j t j 1 for the temporal variable. We first show how we approximate the derivatives with respect to x using the compact FD schemes. Sixth order approximations of the first and second derivatives with respect to x, at interior nodes can be obtained using the following schemes (see [6] for details) u i C1;j u i 1;j 2x C 1 9 u i C2;j y i 2;j 4x ; (4) 2 11 u 00 i 1;j C u 00 i;j C 2 11 u 00 iC1;j D 12 11 where u i;j D u.x i ; t j / and the primes denote differentiation with respect to x. We apply the compact FD approximation for the first and second derivatives given by (4) and (5) respectively, at the interior nodes .i D 2; : : : ; N x 1/. Since we know boundary conditions at i D 1 and i D N x , the compact FD schemes must be adjusted for the nodes near the boundary points. In order to maintain the order O.h 6 / accuracy at the boundary points as in the interior points and to maintain the same tridiagonal format, we use the following one sided scheme at i D 2, and when i D N x 1 we use Similarly, for the second derivatives, we use at i D 2 and at i D N x 1. Using the above equations, the equations for approximating the first and second order derivatives can be expressed as where A x D : : : : : : ;j D OEu 0 2;j ; u 0 3;j ; : : : ; u 0 N x 2;j ; u 0 N x 1;j T ; U 00 ;j D OEu 00 2;j ; u 00 3;j ; : : : ; u 00 N x 2;j ; u 00 N x 1;j T From equations (10) and (11), U 0 and U 00 can be expressed as; where Next we show how we use the compact FD schemes to approximate the time derivatives. Approximation of the first derivative at interior points is given by where the dots denote differentiation with respect to time and t is the distance between successive nodes. We adjust the schemes for i D 1; 2; N t 1 and N t with the following one sided schemes respectively.
The equation for approximating the first time derivative is obtained by combining equations (15) -(19) and is given by where A t D A x and B t D 1 t From equation (20), P U can be expressed as e j;k u.x i ; t k /; i D 2; 3; : : : ; N x 1; j D 1; 2; : : : where E t D A 1 t B t and e j k are the elements of E t and U i; D OEu i;1 ; u i;2 ; : : : ; u i;N t . To solve the PDE (1) we start by linearizing it by using the quasilinearization method which was proposed by Bellman and Kalaba [12]. It is convenient to split the function G in (1) into its linear and nonlinear components and rewrite the governing equation in the form, LOEu; u 0 ; u 00 C N OEu; u 0 ; u 00 P u D 0; where the dot and primes denote the time and space derivatives, respectively. L is a linear operator and N is a non-linear operator. If we assume that the difference u rC1 u r and all its space derivatives is small, then we can approximate the non-linear operator N using the linear terms of the Taylor series and hence N OEu; u 0 ; u 00 N OEu r ; u 0 r ; u 00 where r and r C 1 denote previous and current iterations respectively. Equation (23) can be expressed as N OEu; u 0 ; u 00 N OEu r ; u 0 r ; u 00 r C 2 X kD0 k;r OEu r ; u 0 r ; u 00 r u .k/ rC1 2 X kD0 k;r OEu r ; u 0 r ; u 00 r u .k/ r (24) where k;r OEu r ; u 0 r ; u 00 r D @N @u .k/ OEu r ; u 0 r ; u 00 r : Substituting equation (24) into equation (22), we get where R r OEu r ; u 0 r ; u 00 r D 2 X kD0 k;r u .k/ r N OEu r ; u 0 r ; u 00 r : Substituting (21) and U rC1;j D OEu rC1;j .x 2 /; u rC1;j .x 3 /; : : : ; u rC1;j .x N x 1 / Since the initial condition is known, then we express equation (27) as where R j D R r OEU r;j ; U 0 r;j ; U 00 r;j C e i;1 U 1 LOEH x ; H xx ˆ1 ;r H x ˆ2 ;r H xx ; j D 2; 3; : : : ; N t : Equation (30) can be expressed as the following .N t 1/.N x 1/ .N t 1/.N x 1/ matrix system 2 6 6 6 6 4 X 2;2 X 2;3 X 2;N t X 3;2 X 3;3 X 3;N t : : : : : : : : : : : :

Numerical experiments
We apply the proposed algorithm to well-known nonlinear PDEs of the form (1) with exact solutions. In order to determine the level of accuracy of the CFDQLM approximate solution, at a particular time level, in comparison with the exact solution we report maximum error which is defined by where Q u.x r ; t / is the approximate solution and u.x r ; t / is the exact solution at the time level t .

Example 3.1 (Fisher's equation). We consider the Fisher's equation which represents a reactive-diffusive system and is encountered in chemical kinetics and population dynamics applications. The equation is given by
subject to the initial condition and exact solution [13] u.x; t / D 1 which is consistent with the quasilinearization discussed in the previous section. Therefore, the linearized equation can be expressed as Applying the compact FD schemes both in x and t, and initial condition, we get Equation (47) can be expressed as 2 6 6 6 6 4 X 2;2 X 2;3 X 2;N t X 3;2 X 3;3 X 3;N t : : : : : : : : : : : :

Example 3.2 (Burgers-Fisher equation).
We consider the generalized Burgers-Fisher equation [14] @u @t C˛u ı @u @x D with initial condition and exact solution where˛,ˇand ı are parameters. In this work the parameters are chosen to be˛DˇD ı D 1.

Example 3.3 (Burgers-Huxley equation). We consider the Burgers-Huxley equation
where˛;ˇ 0 are constant parameters, ı is a positive integer (set to be ı D 1 in this study) and 2 .0; 1/. The exact solution subject to the initial condition The general solution (57) is reported in [15,16].

Example 3.4 (Fitzhugh-Nagumo equation). The Fitzhugh-Nagumo equation is given by
with initial condition This equation has the exact solution [17] u.x; t / D 1 2 where˛is a parameter.

Results and discussion
In this section we present the results of the implementation of the CFDQLM on the nonlinear evolution equations mentioned above. In Tables 1 -4 we show the maximum errors between the exact solutions and the CFDQLM results for the Fisher's, Fishers-Burger, Burgers-Huxely and Fitzhugh-Nagumo equations respectively. The results were obtained with grid points N x D 30 for the space variable x. As shown in the tables, for the time variable t, the number of grid points N t varied from 10 to 50. It can be seen that after N t D 40 there is a small change in the error . This means that N t D 40 is sufficient to give sufficient accuracy. As shown by the results, the method gives highly accurate results of at least 10 12 which is achieved with remarkably few grid points in both x and t . The same observation is seen in Figure 5 which shows the errors for each of the equations considered. It can be seen that the average error is around 10 12 for the values of t considered. Another remarkable observation is that the method gives accurate results even for t > 1. A lot of methods loose significant accuracy when t > 1.
A graphical comparison between the CFDQLM results and the exact solutions is shown in Figures 1-4 for the Fisher's, Burgers-Fisher, Burgers-Huxely, Fitzhugh-Nagumo equations respectively. A good agreement is observed between the two sets of results. Figure 6 shows the convergence plots of the CFDQLM for the the four equations. It can be seen that the method reaches full convergence after 4 iterations. The convergence criteria is defined as the infinity norm of the difference between the solution at current iteration and previous iteration, ie, jjU rC1 U r jj 1 ; where r and r C 1 denote previous and current iterations. This shows that the method quickly converges to the true solution.

Conclusion
In this paper, we have extended the use of compact finite difference schemes to both space and time variable when solving parabolic differential equations. This idea was tested on nonlinear evolution equations, namely Fisher's, Burgers-Fisher, Burgers-Huxley and Fitzhugh-Nagumo equations. The quasilinearization technique was used to first linearize the equations. The results obtained show that the method is highly accurate. A noteworthy result is that the high accuracy is obtained using few grid points on both space and time variables. This makes the method more computationally efficient. For future work, we plan to extend the application of the method to: -systems of PDEs -PDEs with higher spatial domains -hyperbolic PDEs