The Differential Quadrature Solution of Reaction-Diffusion Equation Using Explicit and Implicit Numerical Schemes

In this paper, two different numerical schemes, namely the Runge-Kutta fourth order method and the implicit Euler method with perturbation method of the second degree, are applied to solve the nonlinear thermal wave in one and two dimensions using the differential quadrature method. The aim of this paper is to make comparison between previous numerical schemes and detect which is more efficient and more accurate by comparing the obtained results with the available analytical ones and computing the computational time.


Introduction
Thermal wave is reaction-diffusion equation that plays an ever-increasing role in the study of material parameters.It has been employed in optical investigations of solids, liquids and gases with photo-acoustic and thermal lens spectroscopy.Thermal waves have also been used to analyze the thermal and thermodynamic properties of materials and image thermal and material features within a solid sample [1].
In the past several decades, there has been greeting activity in developing numerical and analytical methods for the thermal wave equation.Due to the nonlinearity and complexity of such problems, only limited cases can be analytically solved [2][3][4][5].Yan applied the projective Riccati equation method to solve Schrodinger equation in nonlinear optical fibers [2].Then Mei, Zhang and Jiang employed the same method to get the exact solutions for some reaction-diffusion problems [3].Abdusalam applied a factorization technique to find exact traveling wave solutions [4].Chowdhury and Hashim obtained analytical solution for Cauchy reaction-diffusion problems using homotopy perturbation method [5].Literature on the numerical solution of reaction-diffusion equations is sparse, and singular perturbation method has been applied to solve reaction-diffusion equations by Puri et al. in [6].David, Curtis and John introduced time integration methods to solve thermal wave propagation [7].Marcus applied finite difference method to study the dynamics of predator-prey interactions [8].As well as, Chen et al. employed the finite element method to solve adjective reaction-diffusion equations [9].Then Christos et al. also applied the same method to solve the problem with boundary layers [10].Meral and Sezgin used this method and finite difference method with a relaxation parameter to solve nonlinear reaction-diffusion equation in one and two dimensions [11].Recently, differential quadrature method has been efficiently employed in a variety of engineering problems [12].Wu and Liu had introduced the generalization of the differential quadrature method to solve linear and nonlinear differential equations [13].Kajal applied differential quadrature and Runge-Kutta method to solve thermal wave, a blow-up and a Brusselator chemical dynamics system [14].Kajal achieved high accuracy.But there are some difficulties in the previous method which are explicit schemes used to update the solution using very small step size due to the limitation of stability condition that leads to more computational cost and lower efficiency.Therefore, Meral applied differential quadrature method and implicit Euler method with Newton method to solve one dimensional density dependent nonlinear reaction-diffusion equation [15].Meral obtained stable solutions, and larger time steps could be used.
In this research, the thermal wave propagation model is solved by using two numerical methods to make comparison between them.In the first method, we used the hybrid technique method of Runge-Kutta fourth order method (RK4) and differential quadrature method (DQM).In the second method, we used the combined algorithm of DQM, Perturbation method of second degree and implicit Euler method.Perturbation method is used to avoid the nonlinear term.The obtained results are compared with the previous analytical ones to complete the comparison between previous different numerical schemes.

Numerical Procedure of Thermal Wave
Propagation of thermal waves through a rectangular plate, is governed by [14]: where: U is a temperature, α and β are diffusion parameters in direction of x and y, respectively, γ is reaction parameter, a and b are plate dimensions in direction of x and y, respectively, U max is a maximum temperature of the system.Along the external boundaries, the temperatures can be described as: , , , 0, , , , , , ,

U a U a y t b f y t a y t x
, 0, , , 0, , , , , , i i i a b f i = are known functions.Then initial temperature may be described as: where ( ) , g x y is a known function.

Numerical Procedure Using First Method (RK4)
The main strategy is to employ DQM to reduce the problem to a system of ordinary differential equations then to apply RK4 to solve the reduced system as follows: 1) Discretize the spatial domain using Chebyshev-Gauss-Lobatto grid points [12], such as: ( ) ( ) 2) Apply the method of differential quadrature in terms of nodal temperature, such that: , , , , , , , , ) are the first and second order weighting coefficients with respect to p [12]. 3) On sustainable substitution from Equations ( 5) into (1), one can reduce the problem to the following system of ordinary differential equations as: where .

Numerical Procedure Using Second Method (Implicit Euler)
The main strategy is to apply perturbation method of second order [17,18] then applying DQ discretization to reduce the problem to a system of ordinary differential equations then applying implicit Euler method to trans-form the previous system to a system of linear algebraic equations as follows: 1) We can solve ( ) subjected to the prescribed to boundary and initial conditions in Equations ( 2) and (3), assuming where 0 1 2 , , U U U are unknowns functions and ε is a perturbation parameter.
The following condition is tested to ensure the convergence condition [19] in previous series in Equation (11). 1 where 0,1, , 2) On sustainable substitution from Equation ( 11) into (10), one can reduce the problem to the following equation.
3) Applying zero order perturbation method such that, Subjected to boundary and initial conditions in Equations ( 2) and ( 3), where differential quadrature method and implicit method are used to reduce Equation (14) to a system of linear algebraic equations such that, By substitution of Equation ( 5) into (14) result that, ( ( ) ( ) 4) First order perturbation method is applied such that, ( ) Subjected to the same boundary and initial conditions in Equations ( 2) and (3), reduced to the following algebraic system in equations 5) Also second order perturbation method is applied such that, ( ) Subjected to the same boundary and initial conditions in Equations ( 2) and ( 3), reduced to the following algebraic system Finally, the series solution can be written as We carry on previous procedure until the specified time is reached.

Results and Discussions for One Dimension Analysis
To ensure the accuracy of the proposed numerical techniques, the thermal wave propagating model is solved using presented methods and compared with the available analytical solution [14,20].Consider a one-dimensional problem of thermal wave propagation along x-direction as: ( ) The exact solution for such problem can be obtained as [20]: To validate the accuracy of numerical results, the following errors [16] are computed, ( ) Root mean square of errors . ..of errors Max. absolute error max , , , ,

Numerical Results of First Method (RK4)
For the numerical computation, the time domain is limited to 0 20 t ≤ ≤ and N = 7.The efficiency of presented techniques is tested by CPU time required when the computation reaches to t = 20 s.Two time step sizes of 0.001 and 0.005 are used for layer marching in the time direction.The numerical results of these two cases are listed respectively in Tables 1 and 2. From these two tables, it can be seen that RK4 method can achieve high accuracy at very small step size at Δt = 0.001 with .Moreover, as shown in Figure 1, as Δt increased slightly to 0.00515 the stability condition will not achieved and the oscillation will occurs in the period 6 25 t ≤ ≤ .On the other hand the efficiency is very small as the CPU time required to reach t = 20 s is much larger.

Numerical Results of Second Method (Implicit Euler)
In the obtained results the advantage of using an implicit scheme has been observed.Stability problems are not encountered due to the use of implicit time integration step and larger time increments can be used, e.g. for t = 30, 3.0 t ∆ = can be taken.Table 3 shows the maximum absolute errors and root mean square of errors for a fixed time (t = 30.0)for various numbers of grid points.The accuracies by using N = 8, 11 are almost the same and there is a drop for N = 15.From the table, DQM is observed to give very good accuracy with a small number of grid points.For N = 15, the drop of accuracy is due to the ill-conditioned Vandermonde-system obtained after the DQM discretization, which is the known nature of DQM for large N [15].Tables 4-6 give the comparison of the DQM solution with the exact solution in terms of maximum absolute error and root mean square of errors for small time levels and for the times tending to steady-state, respectively.The computations are carried out with N = 11 and it is seen to be enough to obtain the solution with five digits accuracy at steady-state.Moreover, Figure 2 shows the absolute error at different times and locations.Also convergence condition in Equation ( 12) is tested achieving higher accuracy at second order perturbation method as shown in Figure 3.

Results and Discussions for Two Dimensions Analysis
Consider also a simple two dimensional problem with    which can be solved exactly as [21]: , , e sin π sin π t U x y t x y − = .The design of the numerical scheme is extended to two dimensions.
Table 7 shows that for 0.01, 5, 4 t N M ∆ = = = , the obtained results agree with the analytical ones [21] in both methods.Also Figure 4 shows that absolute error for the hybrid method absolute error <0.01, and for implicit Euler, the absolute error 3 0.6 10 − × .

Conclusion
Throughout this study, thermal wave propagation model which is the type of reaction-diffusion equations is solved by using DQM for space discretization and two different time-integration schemes.Moreover, one can use a small number of discretization points, which lead to higher accuracy.Also for the nonlinear wave equation, the use of DQM with non-uniform grid discretization increases the accuracy and stability of solution.The resulting system of ordinary differential equations is solved by using two different time integration schemes in order  to make comparison between two methods and detect which of them is better.The numerical results obtained in this paper ensure that the problems have small desired time to reach it.Thus they have very small step size which is preferred and use RK4 to solve the system of ordinary differential equations in order to decrease the computational time.On the other hand, the problems which have high desired time to reach it, thus have large incremental time (stiff problems) which are preferred and use implicit Euler with perturbation method to solve the system of ordinary differential equations.

Figure 1 .
Figure 1.Absolute errors at different times and locations.

Figure 2 .Figure 3 .
Figure 2. Absolute errors at different times and locations.

Table 7 .
Root mean square of errors for two dimensional thermal wave propagation at Δt = 0.01, N = 5, M = 4.

Figure 4 .
Figure 4. Absolute error at t = 0.1 s and different locations of x, y.