Numerical Method Using Cubic B-Spline for a Strongly Coupled Reaction-Diffusion System

In this paper, a numerical method for the solution of a strongly coupled reaction-diffusion system, with suitable initial and Neumann boundary conditions, by using cubic B-spline collocation scheme on a uniform grid is presented. The scheme is based on the usual finite difference scheme to discretize the time derivative while cubic B-spline is used as an interpolation function in the space dimension. The scheme is shown to be unconditionally stable using the von Neumann method. The accuracy of the proposed scheme is demonstrated by applying it on a test problem. The performance of this scheme is shown by computing and error norms for different time levels. The numerical results are found to be in good agreement with known exact solutions.


Introduction
This study is concerned with the numerical solution of strongly coupled reaction-diffusion system using cubic B-spline. Reactiondiffusion system arises in the study of biology, chemistry and population dynamics. It can be used to describe a mathematical model of a class of chemical exchange reaction that arises in the transport of ground water in an aquifer [1]. The mathematical formulation of the problem is: with the initial conditions and the boundary conditions u x (0,t)~a 1 (t), u x (1,t)~b 1 (t), v x (0,t)~a 2 (t), v x (1,t)~b 2 where u~u(x,t) and v~v(x,t) are the concentrations of the two substances in the interaction and the constants a,b,c are such that aw0,bw0, c=0. The following consistency conditions hold a 1 (0)~u' 0 (0), b 1 (0)~u' 0 (1), a 2 (0)~v' 0 (0), b 2 (0)~v' 0 (1): The global solutions of such a system of equations have attracted the attention of several researchers [2][3][4][5][6]. Researchers have also investigated the existence, uniqueness and boundedness of the global solution in bounded and unbounded region [3,6]. Cao and Sun [1] derived a finite difference scheme by the method of reduction of order for the numerical solution of strongly coupled reaction-diffusion system with Neumann boundary values conditions. They proved the solvability and convergence by using the energy method. Several researchers focused on analytical solutions of nonlinear equations by using approximate analytical methods. Examples include Ghoreishi et al [7] who obtained the analytical solution for a strongly coupled reaction-diffusion system by using the Homotopy Analysis Method. The solution of the system was calculated in the form of an infinite series with easily computed components. This method cannot always guarantee the convergence of approximate series [8,9] and the method depends on choosing the proper linear operator and initial guesses. The study of spline functions is a key element in computer aided geometric design [10] and also several other applications. It has also attracted attention in the literature for the numerical solution of linear and non-linear system of second-order boundary value problems [11,12] that arise in science and engineering. Some researchers have considered spline collocation method for diffusion problems [15][16][17][18][19]. Advection-diffusion equation arises frequently in the study of mass, heat, energy and vorticity transfer in engineering. Bickley [13] introduced the idea of using a chain of low-order approximation (cubic splines) rather than a global highorder approximation to obtain better accuracy for a linear ordinary differential equation. Fyfe [14] used the method proposed by Bickley [13] and conducted an error analysis. It was concluded that the spline method is better than the usual finite difference scheme because it has the flexibility to obtain the solution at any point in the domain with greater accuracy. The numerical solution of some partial differential equations can be obtained using B-spline functions of various degrees which can provide simple algorithms. As an example, a combination of finite difference approach and cubic B-spline method was used to solve the Burgers' equation [15], heat and wave equation [16,17], advection-diffusion equation [18] and coupled viscous Burgers' equation [19].
Sahin [24] presented the B-spline methods in several degrees for the solution of following non linear reaction-diffusion system LU Lt~a The finite element method was employed for the solution of reaction-diffusion systems and the time discretization of the system was achieved by the using Crank-Nicolson formulae. The system was considered only reaction-diffusion but the problem we propose to investigate is a strongly reaction-diffusion system with an extra term U xx in the second equation of system (1).
In this paper we aim to apply the combination of finite difference approach and cubic B-spline method to solve the system (1). Some researchers have utilized the B-spline collocation methods to solve systems of differential equations but so far as we are aware not the system (1). A usual finite difference scheme is used to discretize the time derivative. Cubic B-spline is applied as an interpolation function in the space dimension. The unconditional stability property of the method is proved by von Neumann method. The feasibility of the method is shown by a test problem and the approximated solutions are found to be in good agreement with the exact solution.
This paper is structured as follows: Firstly, we discuss a numerical method incorporating a finite difference approach with cubic B-spline. Secondly, the von Neumann approach is used to prove the stability of method and compare the numerical solution with exact solution of system (1). Thirdly, a test problem is considered to show the feasibility of the proposed method. Finally, the conclusion of this study is given.

Description of Cubic B-Spline Collocation Method
In this section, we discuss the cubic B-spline collocation method for solving numerically the strongly coupled reaction-diffusion system (1). The solution domain aƒxƒb is equally divided by knots x i into n subintervals ½x i ,x iz1 , i~0,1,2,:::,n{1 where a~x 0 vx 1 v:::vx n~b . Our approach for strongly coupled reaction-diffusion system using collocation method with cubic Bspline is to seek an approximate solution as [20] where C i (t) and D i (t) are (time dependent) quantities which are to be determined for the approximated solutions U i (x,t) and V i (x,t) to the exact solutions u(x,t) and v(x,t) respectively, at the point (x i ,t j ) and B 3,i (x) are cubic B-spline basis functions which are defined by [21] where h~(b{a)=n. The approximations U j i and V j i at the point (x i ,t j ) over subinterval ½x i ,x iz1 can be defined as where i~0,1,2,:::,n. So as to obtain the approximations to the solutions, the values of B 3,i (x) and its derivatives at nodal points are required and these derivatives are tabulated in Table 1.
Using approximate functions (6) and (7), the values at the knots of U j i and V j i and their derivatives up to second order are determined in the terms of time parameters C j k and D j k as The approximations of the solutions of the system (1) at t th jz1 time level can be given by as [22] where a i ,i~1,2,3,4 are constants and the subscripts j and jz1 are successive time levels, j~0,1,2,::::. Discretizing the time derivatives in the usual finite difference way and rearranging the equations, we obtain Table 1. Values of B 3,i (x) and its derivatives.
where Dt is the time step. It is noted that the system becomes an explicit scheme when h~0, a fully implicit scheme when h~1, and a Crank-Nicolson scheme when h~0:5 [15] with time stepping process being half explicit and half implicit. In this paper, we use the Crank-Nicolson approach. Hence, (10) becomes The system thus obtained on simplifying (11) after using (8) consists at the time level t~t jz1 . So as to obtain a unique solution to the resulting system, four additional linear equations are required. Thus, equation (7) is applied to the boundary conditions (3) to obtain From equations (11) and (12), the system becomes a matrix system of dimension 2(nz3)|2(nz3) which is a bi-tridiagonal system that can be solved by the Thomas Algorithm [23]. From equation (10) and (12), the system can be written in the matrix vector form as: where b~½a 1 (t jz1 ),0,0,:::,b 1 (t jz1 ),a 2 (t jz1 ),0,0,:::,b 2 (t jz1 ) T ,j~0,1,2,:::: and M is an 2(nz3)|2(nz3) dimensional matrix given by   Also the entries in sub-matrices A i ,i~1,2,3,:::,8 have the following form

Initial state
After the initial vectors C 0 and D 0 have been computed from the initial conditions, the approximate solutions U jz1 i and V jz1 i at a particular time level can be calculated repeatedly by solving the recurrence relation [19].
Now on inserting the trial solutions (one Fourier mode out of the full solution) at a given point x m C j m~A j j exp(imyh) and D j m~B j j exp(imyh) into equation (16) and rearranging the equations, where A and B are the harmonics amplitudes, y is the mode number, h is the element size and i~ffi ffiffiffiffiffiffi ffi {1 p we get ½X 2 ziY j jz1~½ X 1 ziY j j ,

Y~0:
On direct calculation of equation (17) we obtain j~X 1 ziY X 2 ziY , For stability, the maximum modulus of the eigen-values of the matrix has to be less than or equal to one. As h~0:5 is used in the proposed scheme, we thus substitute the value of h into equation (18) and after some algebraic calculation, it can be noticed that Thus, from (19), the proposed scheme for strongly coupled reaction-diffusion equations is unconditionally stable since the modulus of the eigen-values must be less than one. This means that there are no constraints on grid size h and step size in time level Dt but we should prefer those values of h and Dt for which we obtain the best accuracy of the scheme.

Results and Discussion
To test the accuracy of present method, one example is given in this section with L ? and relative L 2 error norms are calculated by and in similar way for the numerical solution V (x,t).
The numerical order of convergence R for both U(x,t) and V (x,t) of present method is obtained by using the formula [19].
where L ? (n) and L ? (2n) are the errors at number of partitions n and 2n respectively. We compare the numerical solution obtained by cubic B-spline collocation method for strongly coupled reaction-diffusion system (1) with known exact solution.

Example 1
In this example, we present a strongly coupled reactiondiffusion system (1) with numerical solution to show the capability and efficiency of cubic B-spline collocation scheme.
We consider the following problem with initial conditions and boundary conditions as follows: It is straightforward to verify that the following are the exact solutions [1] u(x,t)~e {t sin 2 px, v(x,t)~e {t cos 2 px, We use the cubic B-spline collocation method (11)- (12) and (14)- (15) to compute the numerical solution of (21)-(23). Table 2 and  Table 3 show the acceptable comparison between numerical and exact solutions at some grid points at t~1:0 with different number of partitions. This problem is tested using different values of h and Dt to show the capability of the proposed method for solving the system (21)- (23). The final time is taken as T~1. The maximum absolute errors of the method at some grid points are comparable with finite difference scheme in [1]. The numerical errors of the proposed method are presented in Table 4 and Table 5 and are also depicted graphically in Figure 1 and Figure 2. The absolute errors of the numerical solution using finite difference scheme in [1] are shown in Table 6 and Table 7 The solutions are also tabulated in Table 8 and Table 9 with different number of partition and different time levels. Results are presented graphically for U(x,t) and V(x,t) in Figure 3 and Figure 4 for 0vtƒ1 respectively.
The order of convergence of the present example is calculated by the use of the formula given in (20) and which is tabulated in Table 10 and Table 11 for U(x,t) and V (x,t) respectively. An examination of these tables indicates the method has a nearly second order of convergence.

Concluding Remarks
In this paper, a numerical method which incorporates a usual finite difference scheme with cubic B-spline is presented for solving the strongly coupled reaction diffusion system. A finite difference approach is used to discretize the time derivatives and cubic Bspline is used to interpolate the solutions at each time level. It is noted that sometimes the accuracy of solution reduces as time increases due to the time truncation errors of time derivative term [19]. However the cubic B-spline method used in this work is simple and straight forward to apply. The computed results show that the cubic B-spline gives reasonable solutions which are comparable with finite difference scheme with smaller space steps. The obtained solution to the reaction diffusion system for various time levels have been compared with the exact solution by finding the L ? and L 2 errors. An advantage of using the cubic B-spline method outlined in this paper is that it can give accurate solutions at any intermediate point in the space direction.