Numerical Discrete-Domain Integral Formulations for Generalized Burger-Fisher Equation

Okey Oseloka Onyejekwe, Beruk Minale, Fikru Habtamu, Tesfaye Amha, Getenet Tamiru, Bethelhem Mengistu, Yohaness Demiss, Nahom Alemseged, Computational Science and Dynamics Systems Group

Computational Science Program, Addis Ababa University, Arat-Kilo Campus, Addis Ababa, Ethiopia.

DOI: 10.4236/am.2020.113012                

Abstract

In this study we use a boundary integral element-based numerical technique to solve the generalized Burger-Fisher equation. The essential feature of this method is the fundamental integral representation of the solution inside the problem domain by means of both the boundary and domain values. The occurrences of domain integrals within the problem arising from nonlinearity as well as the temporal derivative are not avoided or transferred to the boundary. However, unlike the classical boundary element approach, they are resolved within a finite-element-type discrete domain. The utility and correctness of this formulation are proved by comparing the results obtained herein with closed form solutions.

Keywords

Burger-Fisher Equation, Domain Integrals, Boundary Element, Finite-Element-Type Discrete Domain, Boundary Integral, Numerical, Nonlinearity

1. Introduction

Burger-Fisher equation describes the interaction between linear diffusion, nonlinear convection and nonlinear logistic-reaction and has wide applications in several areas of engineering and applied science. Several numerical methods have been proposed for its solution. A few of them include cubic B-spline interpolation quasi-interpolation technique (Zhu and Wang [1] ), variational iteration method (Abdou and Soliman [2], Zhao et al. [3] ) non-standard finite difference scheme (Mickens and Gumel [4] ), spectral collocation approach (Javidi [5] ), homotopic mapping methods (MO [6] ), direct discontinuous Galerkin method (Zhang and Zhang [7] ) and other methods (Wazwaz and Gorgius [8], Wazwaz [9] ), Wazwaz [10] ).

A major thrust of this paper is to further display the utility of one of the emerging families of domain-based integro-differential numerical techniques (Onyejekwe [11] [12], and [13] ) for the solution of transient nonlinear partial differential equations. The approach is based on discretizing the governing differential equation into a local or element based system of equations which finally yield a sparse coefficient matrix thereby making it competitive with FEM for problems of similar rigor.

For the problem set up, all the nodal points are inter-connected just like in the finite element method (FEM). In case of a linear element, only two nodal points are required for discretization in each sub-domain. Both the dependent variable and its spatial derivative constitute the two degrees of freedom at each node. Within each element, the problem constants and parameters represent those of the element under consideration. Thus, the method is completely domain-discretized and element based. It should however be realized that the implementation of such an idea, requires a deeper numerical insight into the properties of the corresponding nonlinear integral operators and their application to elements resulting from a discretized problem domain. The paper is aimed at further strengthening this approach by applying these ideas to a highly nonlinear, transient one-dimensional problem which more often than not, poses a considerable numerical challenge to classical boundary element method (BEM) application.

In this paper, we consider the generalized Burger-Fisher equation described as:

$D\frac{{\partial }^2\theta }{\partial x^2}=\frac{\partial \theta }{\partial t}+\alpha {\theta }^d\frac{\partial \theta }{\partial x}+\mu \theta \left({\theta }^d-1\right),0\le x\le 1,t>0$ (1)

where, D is the diffusion coefficient, and $\alpha ,\mu$ and d are the equation parameters. The exact solution is

$\theta \left(x,t\right)={\left[\frac{1}{2}+\frac{1}{2}\tanh \left[{\pi }_1\left(x-{\pi }_{2}t\right)\right]\right]}^{1/d},t\ge 0$ (2)

with initial and boundary conditions:

$\theta \left(x,0\right)={\left[\frac{1}{2}+\frac{1}{2}\tanh \left({\pi }_{1}x\right)\right]}^{1/d}$ (3)

$\theta \left(0,t\right)={\left[\frac{1}{2}+\frac{1}{2}\tanh \left({\pi }_1{\pi }_{2}t\right)\right]}^{1/d},t\ge 0$ (4)

$\theta \left(1,t\right)={\left[\frac{1}{2}+\frac{1}{2}\tanh \left[{\pi }_1\left(1-{\pi }_{2}t\right)\right]\right]}^{1/d},t\ge 0$ (5)

where, ${\pi }_1=\frac{-\alpha d}{2D\left(1+d\right)},{\pi }_2=\frac{\alpha }{\left(1+d\right)}+\frac{\left(D*\mu \left(1+d\right)\right)}{\alpha }$

2. Mathematical Formulation

The numerical approach adopted herein is initiated by implementing the integral replication of Equation (1). This is achieved by recasting the governing differential equation into its Poisson form and applying the Green’s second identity to obtain:

${\displaystyle \underset{x_x}{\overset{x_L}{\int }}\left[\theta \delta \left(x-x_i\right)-\frac{G}{D}\left(\alpha {\theta }^d\frac{\partial \theta }{\partial x}+\frac{\partial \theta }{\partial t}+\mu \theta \left({\theta }^d-1\right)\right)\right]\text{d}x}=\left[\theta \frac{\text{d}G}{\text{d}x}-G\frac{\text{d}\theta }{\text{d}x}\right]$ (6)

where, $G\left(x,x_i\right)=\left[\frac{\left|x-x_i\right|+K}{2}\right]$ represents the free-space Green’s function and

K is an arbitrary constant. We may mention in passing that when the Poisson equation is solved using the boundary integral technique, a domain integral exists. The way and manner of handling the domain integral constitutes the theme of several boundary element method (BEM) solution techniques. Equation (6) is solved on each element of the problem domain and then assembled in a finite-element sense.

Major aspects of the scheme are as follows: The formulation is motivated by two important characteristics with some inherent properties of BEM formulation. The first is the “elementization” or discretization of the problem domain. This finite-element-like formulation guarantees local compact support. This attributes guarantees that a slim and banded coefficient matrix which can be stored efficiently even for large field problems and is easily amenable to numerical solution. The second characteristic is BEM’s straightforward formulation and accuracy

The problem domain is discretized into elements and within each element, the dependent variables and their functions are approximated by interpolating functions in space to guarantee the continuity of the dependent variable and its derivative throughout the discrete elements.

$\begin{array}{l}\frac{\partial \theta \left(x,t\right)}{\partial x}=\phi \left(x,t\right)\approx {\Omega }_1{\phi }_1+{\Omega }_2{\phi }_2,\theta \left(x,t\right)\approx {\Omega }_1{\theta }_1+{\Omega }_2{\theta }_2\\ \frac{\partial \theta \left(x,t\right)}{\partial t}\approx {\Omega }_1{\left\{\frac{\partial \theta \left(x,t\right)}{\partial t}\right\}}_1+{\Omega }_2{\left\{\frac{\partial \theta \left(x,t\right)}{\partial t}\right\}}_2\end{array}$ (7)

Substitute Equation (7) into Equation (6) to give:

$\begin{array}{l}-\lambda \theta \left(x_i,t\right)+{\left[\theta \frac{\text{d}G}{\text{d}x}-G\frac{\partial \theta }{\partial x}\right]}_{x_0}^{x_L}+{\displaystyle \underset{x_0}{\overset{x_L}{\int }}\frac{G}{D}(\alpha \left\{{\Omega }_1{\theta }_1^d+{\Omega }_2{\theta }_2^d\right\}\left\{{\Omega }_1{\phi }_1+{\Omega }_2{\phi }_2\right\}\begin{array}{c}\\ \end{array}}\\ +\left\{{\Omega }_1\frac{\text{d}{\theta }_1}{\text{d}t}+{\Omega }_2\frac{\text{d}{\theta }_2}{\text{d}t}\right\}+\left\{{\Omega }_1\left(\mu {\theta }_1\left[{\theta }_1^d-1\right]\right)+{\Omega }_2\left(\mu {\theta }_2\left[{\theta }_2^d-1\right]\right)\right\})\text{d}x=0\end{array}$ (8)

In order to make Equation (8) amenable to element-by-element computation, we covert the spatial linear interpolating functions into local coordinate that has its origin at node 1 of any generic element. viz: ${\Omega }_1\left(\xi \right)=1-\xi$, ${\Omega }_2\left(\xi \right)=\xi$, $\xi =\left(x-x_i\right)/l$, $0\le \xi \le 1$. The spatial derivative of the free space Green function is $\text{d}G\left(x,x_i\right)/\text{d}x=0.5\left[H\left(x-x_i\right)-H\left(x_i-x\right)\right]$, where H(), is the Heaviside function. A general representation of the numerical solution within the problem domain in discrete form is obtained by representing Equation (8) on each of the elemental nodes.

For node 1 Equation (8) becomes:

$\begin{array}{l}\frac{1}{2}\left\{\upsilon \left(-{\theta }_1+{\theta }_2+{\phi }_1-\left(l+l_m\right){\phi }_2\right)\right\}\\ +\frac{1}{2}\{l{\displaystyle \underset{0}{\overset{1}{\int }}\left(l\xi +1\right)[a\left({\Omega }_1{\theta }_1^d+{\Omega }_2{\theta }_2^d\right)\left({\Omega }_1{\phi }_1+{\Omega }_2{\phi }_2\right)\begin{array}{c}\\ \end{array}}\\ +\left({\Omega }_1\frac{\text{d}{\theta }_1}{\text{d}t}+{\Omega }_2\frac{\text{d}{\theta }_2}{\text{d}t}\right)+\left({\Omega }_1\left(\mu {\theta }_1\left[{\theta }_1^d-1\right]\right)+{\Omega }_2\left(\mu {\theta }_2\left[{\theta }_2^d-1\right]\right)\right)]\text{d}\xi \}=0\end{array}$ (9)

Similarly at node 2 we obtain

$\begin{array}{l}\frac{1}{2}\left\{\left({\theta }_1-{\theta }_2+\left(l+l_m\right){\phi }_1-{\phi }_2\right)\right\}\\ +\frac{1}{2}\{l{\displaystyle \underset{0}{\overset{1}{\int }}\left[l\left(1-\xi \right)+1\right][a\left({\Omega }_1{\theta }_1^d+{\Omega }_2{\theta }_2^d\right)\left({\Omega }_1{\phi }_1+{\Omega }_2{\phi }_2\right)\begin{array}{c}\\ \end{array}}\\ +\left({\Omega }_1\frac{\text{d}{\theta }_1}{\text{d}t}+{\Omega }_2\frac{\text{d}{\theta }_2}{\text{d}t}\right)+\left({\Omega }_1\left(\mu {\theta }_1\left[{\theta }_1^d-1\right]\right)+{\Omega }_2\left(\mu {\theta }_2\left[{\theta }_2^d-1\right]\right)\right)]\text{d}\xi \}=0\end{array}$ (10)

Equation (9) and Equation (10) are put in compact matrix form to yield:

$\left(R_{ij}{\theta }_j+L_{ij}{\phi }_j\right)+\alpha {\Lambda }_{ijk}{\theta }_j{\phi }_k+\frac{T_{ij}}{D}\left(\frac{\text{d}{\theta }_j}{\text{d}t}+\mu {\varphi }_j\left\{{\varphi }_j^d-1\right\}\right)=0\text{where}i,j,k=1,2$ (11)

Equation (11) is system of first order nonlinear differential equations in time which is solved for $\theta$ and $\text{d}\theta /\text{d}x$ at the nodes. The temporal term in Equation (11) still needs to be discretized. There are several ways of achieving this but for the time being, we use a two-level time scheme. This approximates the time scheme, i.e., $t=t_{m+\alpha }=t_m+z\Delta t$. As a result, Equation (11) becomes:

$\begin{array}{l}z\left[\left(R_{ij}{\theta }_j^{\left(m+1\right)}+L_{ij}{\phi }_j^{\left(m+1\right)}\right)+a{V}_{ijl}{\theta }_j^{\left(m+1\right)}{\phi }_j^{\left(m+1\right)}\right]\\ +\left(1-z\right)\left[\left(R_{ij}{\theta }_j^{\left(m\right)}+L_{ij}{\phi }_j^{\left(m\right)}\right)+a{V}_{ijl}{\theta }_j^{\left(m\right)}{\phi }_j^{\left(m\right)}\right]\\ +T_{ij}\left(\frac{1}{\Delta t}\left\{{\theta }_j^{\left(m+1\right)}-{\theta }_j^{\left(m\right)}\right\}+\left(z\right)\mu {\theta }_j^{\left(m+1\right)}\left\{{\theta }_j^{d\left(m+1\right)}-1\right\}+\left(1-z\right)\mu {\theta }_j^{\left(m\right)}\left\{{\theta }_j^{d\left(m\right)}-1\right\}\right)\\ \equiv s_i=0\end{array}$ (12)

where z is the weighting factor and subscripts m and m + 1 denote the previous and current time levels and the time step is defined as $\Delta t=t_{m+1}-t_m$. We elect to use the Newton Raphson’s technique to resolve the nonlinearity in Equation (12). The computational procedure starts with a known estimate (guessed estimate) of an unknown quantity say:

${\left\{{\theta }_j^{\left(m+1\right)},{\phi }_j^{\left(m+1\right)}\right\}}^{\text{T}}$, and update it to ${\left\{{\theta }_j^{\left(m+1,k\right)}+\Delta {\theta }_j^{\left(m+1,k+1\right)},{\phi }_j^{\left(m+1,k\right)}+\Delta {\phi }_j^{\left(m+1,k+1\right)}\right\}}^{\text{T}}$, eventually, the increment is accounted for by the solution of the matrix equation:

$J_{ij}^{\left(m+1,k\right)}\left\{\begin{array}{c}\Delta {\theta }_j^{\left(m+1,k+1\right)}\\ \Delta {\phi }_j^{\left(m+1,k+1\right)}\end{array}\right\}=-s_i^{\left(m+1,k\right)}$ (13)

where $J_{ij}^{\left(m+1,k\right)}$ is the coefficient or the Jacobian matrix. The key item here is the determination of the Jacobian matrix. The solution of Equation (13) yields the increments used in subsequent iterations.

The Jacobian for Equation (12) is:

$\begin{array}{l}\frac{\partial s_i}{\partial {\theta }_j^{\left(m+1\right)}}=\left[\alpha \left\{R_{ij}+a{V}_{ijl}{\phi }_j^{\left(m+1\right)}+\left(2d\right)\mu \left(2{\theta }^{\left(2d-1\right)}-1\right)\right\}+\frac{T_{il}}{\Delta t}\right],\\ \frac{\partial s_i}{\partial {\upsilon }_j^{\left(m+1\right)}}=\alpha \left[L_{ij}+a{V}_{ijl}{\theta }_j^{\left(m+1\right)}\right],\text{where}i,j,l=1,2\end{array}$ (14)

The main advantage of these systems is that they have a direct construction, in addition they are local, in the sense that the dependent variables depend only on the values in the neighborhood of any specific node under consideration.

3. Numerical Experiments and Discussion

In order to test the accuracy and utility of the current formulation, numerical solutions of some test problems are compared with analytical results for different problem parameters. The following parameters were found to provide stable and accurate results for all the problems tested herein. $\Delta t=0.001$, $d=1$, $N\left(\text{total}\text{grid}\text{points}\right)=21$. Relative and absolute errors $abs\left(\left({\theta }_{num}-{\theta }_{exact}\right)/{\theta }_{num}\right)$, $abs\left({\theta }_{num}-{\theta }_{exact}\right)$ are used to compute the error magnitude.

Table 1 displays the characteristics of the scalar profile for equal weighting of the diffusion and nonlinear inertia terms. The ability of the integral formulation to deal with such cases is demonstrated by the magnitudes of the absolute and relative errors. It can be observed that at each computational point, the error magnitudes display a slight increase as time progresses. This is probably due to the influence of the reaction term.

On the other hand, when the diffusion term is dominant, as shown in Table 2, the damping effect on the scalar profile is shown by the smaller magnitude of the relative and absolute errors.

Table 1. Comparison of numerical and analytic results for $0\le x\le 1,\alpha =1.0$ and $\mu =1.0,D=1.0$.

Table 3 and Table 4 are for cases where the diffusion term again plays a more prominent role that than the nonlinear inertia term. For both cases there is a damping effect extends to the negative nonlinear reaction term as well. The error terms for both cases are significantly less than the two previous cases.

The ability of the formulation to handle relatively steep gradients as the value of the reaction term is significantly increased is shown in Table 5. Convergent results are noticeable as can be seen by the homoclinic decrease of error magnitudes as time increases.

For Table 6 and Table 7, the diffusion coefficient was set at 0.5; $\alpha =1.0$ in Table 6. This results in the dominance of the inertia term over diffusion and the equation adopts a more hyperbolic form. The reverse is the case in Table 7 where $\alpha =0.001$. The previous case results in a more computational challenge, because of the overarching influence of the nonlinear inertia term. This effect is displayed by the magnitudes of the errors displayed in Table 6 and Table 7. As can be seen greater errors are displayed in Table 6 than the Table 7. The influence of the reaction terms remain the same for both cases.

Table 2. Comparison of numerical and analytic results for $0\le x\le 1,\alpha =1.0$ and $\mu =0.0,D=1.0$.

Table 3. Comparison of numerical and analytic results for $0\le x\le 1,\alpha =0.1$ and $\mu =-0.0025,D=1.0$.

Table 4. Comparison of numerical and analytic results for $0\le x\le 1,\alpha =0.001$ and $\mu =-0.0025,D=1.0$.

Table 5. Comparison of numerical and analytic results for $0\le x\le 1,\alpha =1.0$ and $\mu =10.0,D=1.0$.

Table 6. Comparison of numerical and analytic results for $-0.5\le x\le 0.5,\alpha =1.0$ and $\mu =0.01,D=0.5$.

Table 7. Comparison of numerical and analytic results $-0.5\le x\le 0.5,\alpha =0.001$ and $\mu =0.01,D=0.5$.

4. Conclusion

In the work reported herein, a boundary-integral domain-discretized formulation with uniform grids and time steps is applied for the solution of the Burger-Fisher equation. The formulation is straightforward and amenable to handling nonlinearity. Encouraging results are obtained for different parameter values with relatively small number of elements. This numerical experience, further justifies the need to amend and adapt BEM to become a more competitive numerical tool.

Acknowledgements

A lot of thanks to the anonymous reviewers who significantly improved the quality of this paper

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zhu, C.-G. and Wang, R.-H. (2009) Numerical Solution of Burgers’ Equation by Cubic B-Spline Quasi-Interpolation. Applied Mathematics and Computation, 208, 260-272. https://doi.org/10.1016/j.amc.2008.11.045
[2]	Abdou, M.A. and Soliman, A.A. (2005) Numerical Solutions of Nonlinear Evolution Equations Using Variational Iteration Method. Journal of Computational and Applied Mathematics, 207, 111-120. https://doi.org/10.1016/j.cam.2006.07.016
[3]	Zhao, C.-Z, Yu, X.J. and Zhu, J. (2010) Direct Discontinuous Galerkin Method for the Generalized Burger-Fisher Equation. Chinese Physics B, 19, Article ID: 0702203.
[4]	Mickens, R.E. and Gumel, A.B. (2002) Construction and Analysis of a Nonstandard Finite Difference Scheme for the Burgers-Fisher Equation. Journal of Sound and Vibration, 257, 791-797. https://doi.org/10.1006/jsvi.2001.4240
[5]	Javidi, M. (2006) Spectral Collocation Method for the Solution of Generalized Burger-Fisher Equation. Applied Mathematics and Computation, 174, 45-52. https://doi.org/10.1016/j.amc.2005.04.084
[6]	Mo, J.-Q. (2009) A Variational Iteration Solving Method for a Class of Generalized Boussinesq Equations. Chinese Physics Letters, 26, Article ID: 060202. https://doi.org/10.1088/0256-307X/26/6/060202
[7]	Zhang, R.-P. and Zhang, L.-W. (2012) Direct Discontinuous Galerkin Method for Generalized Burgers-Fisher Equation. Chinese Physics B, 21, Article ID: 090206. https://doi.org/10.1088/1674-1056/21/9/090206
[8]	Wazwaz, A.M. and Gorguis, A. (2004) An Analytic Study of Fishers Equation by Using Adomian Decomposition Method. Applied Mathematics and Computation, 154, 609-620. https://doi.org/10.1016/S0096-3003(03)00738-0
[9]	Wazwaz, A.M. (2005) The Tanh Method for Generalized Forms of Nonlinear Heat Conduction and Burgers-Fisher Equation. Applied Mathematics and Computation, 169, 321-338. https://doi.org/10.1016/j.amc.2004.09.054
[10]	Wazwaz, A.M. (2008) Analytic Study of Burgers, Fisher, Huxley Equations and Two Combined Forms of These Equations. Applied Mathematics and Computation, 195, 754-761. https://doi.org/10.1016/j.amc.2007.05.020
[11]	Onyejekwe, O.O. (2005) An Efficient Green Element Algorithm for Radial Flow. Applied Mathematics and Computation, 165, 635-645. https://doi.org/10.1016/j.amc.2004.04.099
[12]	Onyejekwe, O.O. (2006) A Note on Green Element Method Discretization for Poisson Equation in Polar Coordinates. Applied Mathematics Letters, 19, 785-788. https://doi.org/10.1016/j.aml.2005.08.025
[13]	Onyejekwe, O.O. (2009) A Modified Boundary Integral Evolution Formulation for the Wave Equation. Advances in Engineering Software, 40, 519-526. https://doi.org/10.1016/j.advengsoft.2008.10.002