An Hermitian Boundary Integral Hybrid Formulation for Nonlinear Fisher-Type Equations
Okey Oseloka Onyejekwe
Computational Science Program, Addis Ababa University, Arat Kilo campus, Addis Ababa, Ethiopia
Email address:
To cite this article:
Okey Oseloka Onyejekwe. An Hermitian Boundary Integral Hybrid Formulation for Nonlinear Fisher-Type Equations. Applied and Computational Mathematics. Vol. 4, No. 3, 2015, pp. 83-99. doi: 10.11648/j.acm.20150403.11
Abstract: This paper explores the application of an Hermitian hybrid boundary integral formulation for handling Fisher-type equations. The Hermite system incorporates the problem unknowns with their space derivatives and as a consequence produces a relatively larger coefficient matrix than the corresponding linear approximation. However by adopting a finite element-like integral numerical procedure, the modified boundary integral formulation otherwise known as the Green element method (GEM) produces slender and sparse coefficient matrices which enhance an efficient solution algorithm. The resulting equations appear in the form of local elemental integral equations whose contributions add up to the coefficient matrix. This process is amply simplified in the Green element method due to the presence of the source point inside an element thereby encouraging integration to be carried out locally and accurately. This so called ‘divide and conquer’ approach is significantly much better than working with the entire matrix especially for nonlinear problems where an encounter with the problem domain can not be totally avoided. Numerical tests are carried out to illustrate the utility of this technique by comparing results obtained from both the Hermite and non-Hermite discretizations. It is observed that for each of the problems tested, not only do the results agree with those from literature, it took the Hermitian approximation fewer number of elements to achieve the same level of accuracy than its non-Hermitian version. However, application of same technique to multi-dimensional problems may not be as straightforward due to the construction and storing of the Hermite system matrix which will not only involve non-trivial operations in terms of a high computational cost but also a compromise in the quality of the numerical solution arising from significant round-off errors.
Keywords: Fisher-Type Differential Equations, Boundary Element Method, Nonlinearity, Hybrid Boundary Integral Method, Green Element Method, Hermite Interpolation
1. Introduction
Cubic Hermite approximation possesses a unique feature namely; both the primary variable and its space derivative are continuous across inter-nodal boundaries of the problem domain. The utility of this approximation rests in the fact that not only is a differentiable interpolation of the dependent variable achieved but also a higher order more accurate scheme using the same nodes as those employed for the simple linear element is guaranteed. It goes without saying that there is also a price to be paid namely; extra computational cost, introduction of round-off errors, as well as storage requirements needed for the solution of additional equations. Nevertheless some of these numerical challenges can be dealt with by GEM’s element-based numerical procedure which produces a banded and sparse coefficient matrix.
Transient diffusion coupled with nonlinear source or reaction terms are typical of Fisher-type equation, a prototypical reaction-diffusion differential equation. Such problems arise in several areas of engineering physics, epidemiology, applied mathematics, biological systems and nuclear reactor physics. The ubiquitous nature of these equations has provided the motivation for numerous theoretical and numerical work in several areas of application and has subsequently provided rich information for numerous problems targeted.
The numerical solutions of Fisher’s equation often exhibit propagating fronts that can be very steep for large values of reaction rate coefficient. For such problems, linear interpolation is inadequate to capture the physics that accompany the rapid change of the scalar profile within a short spatial interval. Problems of this type,can therefore be dealt with by incorporating higher order interpolation into the formulation. Hence from a modeling perspective, the work presented herein offers a good canonical model of a cooperation between a hybrid boundary integral numerical formulation and the first- order cubic Hermitian polynomials that are guaranteed to ensure the continuity of both the primary variable and its flux across interzonal boundaries.
Fisher [1] was the first to study the movement of a viral mutant in a habitat by adopting the nonlinear reaction-diffusion equation. Fisher’s equation illustrates a case where the diffusion term is linear and the reaction term is quadratic [2,3]. However it is known that there are other versions of it where the diffusion term is nonlinear [4,5]. Due to its complexity, a closed form solution of the governing equation in general is not achievable and as a result, most of the efforts directed at solving this nonlinear reaction-diffusion equation is by perturbation theory or more usually by numerical methods [6,7]. Despite these attempts Fisher’s equation still remains a challenging problem from a numerical point of view. This has a lot to do with the high steep fronts created for cases where the reaction becomes more prominent than diffusion . As a consequence, solutions of Fisher equation have elicited a vast area of research in fields such as finite difference, finite element techniques, boundary element method, adaptive and non-adaptive algorithms and spectral techniques [8,9]. One of the earliest numerical solutions using a pseudo-spectral approach was given by Gazdag and Canosa [10]. This was later followed by Mittal and Kumar [11]. More rigorous mathematical approaches were also applied in this field. Prominent among them are the following: the non-standard finite difference methods [12,13], boundary integral schemes[14,15] and the Adomian method [16]. Further attempts include the B-spline Galerkin approach of Sahin et al.[17]. Meral and Tezgin [18] applied the dual reciprocity boundary element technique as well as the differential quadrature method. They found the dual reciprocity boundary element method quite useful in transforming domain integrals resulting from the nonlinearity of the governing differential equation into equivalent boundary integrals .
In the work reported herein we consider the numerical solutions of Fisher-type equations with a hybrid boundary integral technique known as the Green element method (GEM) first mentioned by Taigbenu [19,20] and further developed by Taigbenu and Onyejekwe[21]. GEM starts by transforming the governing partial differential equation into its integral analog via the Green’s second identity. The resulting singular boundary integral equations are then evaluated element-by-element. This is similar in implementation to the finite element method (FEM) yet it maintains BEM second order accuracy. Both the dependent variable and its flux are obtained at every node.
GEM like every other boundary element- based technique requires the fundamental solution of an auxiliary equation of the governing differential equation. But unlike other BEM -based technique, it does not convert to boundary integrals those other terms of the original PDE that are not considered when the fundamental solution is derived. Among the equations that come with domain integrals are those which have always posed considerable numerical challenges to boundary integral techniques and which do not always admit tidy solutions e.g. nonlinear and heterogeneous problems, transient problems, problems that involve body force and source/sink terms. These encompass a majority of those problems which form the bulk of real-life engineering computations. BEM could still be broken down into zones and compartments to deal with domain integrals and yet maintain its boundary-only thrust, but such an approach would seem to produce acceptable results only for problems involving weak nonlinear terms [22,23]. On the contrary, GEM deals with a large number of sub-regions or elements which at the extreme are called the ‘Green elements’ (Onyejekwe[24])
GEM facilitates the fundamental solution of the Laplace equation and accurately accounts for the remaining terms of the governing equation by adopting domain integration. Similar approach can be found in the BDIM of Hribersek and Skerget [ 25] for the solution of incompressible viscous fluid flow as well as the iterative solution schemes of DRM-MD technique of Portapila and Power [26] where they applied DRM-MD to iteratively solve the linear systems of equations arising from the dual reciprocity method in multi-domains.
2. Computational Approach
We consider the generalized Fisher equation given as:
(1)
where is a real valued function of the space variable and time
With initial and boundary conditions given as:
where are constants. The application of GEM to obtain a discrete analog of equation(1) involves the following procedure:
a) Solution of a prescribed auxiliary equation of (1), to obtain the Green function, also known as the fundamental solution or simply the free space Green function. Classical BEM approach requires that all the terms of the governing partial differential equation that
can not be included in the auxiliary equation be put in the domain integral. In GEM there is no question of reconverting those integrals back to the boundary. They are left in the domain and dealt with as such.
b) The Green’s second identity is invoked and appropriate substitutions are made from the auxiliary equation, the free space Green’s function and the governing differential equation offer the integral replication of the governing differential equation.
c) The point of departure from a classical boundary-based approach is effected in the discretization of the problem domain into suitable partitions or elements over which the scalar variables are prescribed.
d) The finite-element component of this hybrid formulation requires that equations be derived for each element and both the boundary and initial conditions be enforced in the system of element equations.
e) All the element equations are assembled into a global matrix, and the known scalar values put in the right hand-side vector with the boundary and continuity conditions enforced.
f) In addition to GEM retaining its boundary integral second order accuracy, finite element implantation guarantees that the coefficient matrix be sparse and banded and easily handled numerically to yield the required dependent variables. As a consequence, the so called ‘local support’ is guaranteed as well as handling problems which require domain integration.
With the above GEM formulation procedure outlined, the next is to set up a time-marching procedure to utilize known data at the beginning of a time step to predict new values at the end of the time step.
Linear interpolation so called zero order continuity across inter-element boundaries has been found to be inadequate for handling problems which exhibit fast changes in the solution profile. Among these are problems that display the formation of boundary layers, and the advection or reaction dominant types. Though the governing equation in such cases remains bounded as the scalar history evolves through time, its spatial derivative develops a steep profile which quickly becomes unbounded and not easy to handle numerically. Current boundary integral literature in this area points to an evolving active research aimed at efficiently handling both the problem domain as well as its boundary for those problems where domain consideration is a necessity [27 ,28, 29,30].
2.1. GEM Hermitian Cubic Interpolation
Let us assume that , then equation (1) takes a form typical of a Fisher-like nonlinear diffusion-reaction equation:
(2)
The singular integral representation of equation (2) for a typical element is formulated along the lines of our previous work [29], and some of the details are restated here for clarity.
(3)
where lambda takes the value of one half when the source node is placed at the end of the element and unity when the source node is placed within the element. The auxiliary differential equation: , and its solution, the so called fundamental solution is: , , the parameter k is an arbitrary constant, which for the purposes of stability is assigned the longest element length, the flux is specified as: , for a more general case, the diffusion coefficient admits nonlinearity and its reciprocal is given as: , with its logarithm specified as: ; while the derivative of the free space Greens function is: where H is the Heaviside function. In our previous work [31] we adopted linear shape functions with zero-order continuity to do the interpolation of these functional values. For the current work however, we concentrate on a higher level of interpolation functions known as Hermitian cubic interpolation which guarantees the continuity of both the primary variable and its derivative across inter-element boundaries.
Any variable is approximated as:
(4a)
The following cubic Hermitian approximations are needed for equation (3)
(4b)
(4c)
(4d)
(4e)
(4f)
where are Hermitian basis functions which are expressed in terms of a local coordinate as:
(5a)
(5b)
(5c)
(5d)
Substitution of all the Hermite approximated terms into equation(3) yields:
(6)
Equation (6) is a system of first order highly nonlinear differential equation. The next chore is to interpolate the temporal derivative terms for the flux and the primary variable terms. The following equation is obtained by adopting a two-level time scheme for both variables.
(7)
Linearization has to be adopted to solve equation (7). Because of the huge computational cost that will be involved in dealing with the derivatives of the dependent variables and computing the Jacobian matrix, the Newton-Raphson procedure is avoided. The Picard procedure is implemented instead to yield:
(8)
The Picard technique employs known estimates of the dependent variable in the coefficient matrix to calculate refined estimates. These are compared with previous values until the difference between the two are less that an a-priorily set convergence estimate. By so doing the chore of determining the derivatives as well as a careful choice of the first estimates are obviated.
2.2. Modification of the Fisher’s Equation
Apart from the quadratic nonlinearity that is inherent in the reaction term for the Fisher’s equation, other sources of nonlinearity can occur to account for cases involving nonlinear boundary conditions, temperature-dependent thermo-elastic properties of a material for both or either the conduction coefficient or heat capacity terms.
The governing equation for a temperature-dependent conductivity, specific heat and density is given as:
(9)
where is the heat capacity term, are the density and specific heat of the material respectively, is the two-dimensional (2D) gradient operator and and are unit vectors in the horizontal and vertical directions. Equation (1) can admit any of the following first (Dirichlet) and second kind (Neumann) boundary conditions:
(10a)
(10b)
Or a combination of Dirichlet and Neumann conditions where both the function value and the normal derivative are specified on the boundary of the domain. This is the Cauchy type boundary condition and is given as:
(10c)
where are specified on the boundary, stands for a unit outward normal vector from the boundary, are different sections of the boundary and s is a time parameter. Another requirement is that the scalar profile be established everywhere in the problem domain at an initial time.
Classical boundary element theory requires a complimentary differential equation be adopted for the solution process. This can be of the type for which is the coordinate of the field node, for the source node . The fundamental solution of this Laplacian operator is given as: The finite element component of this hybrid computational technique, requires the problem domain to be discretized into suitable polygonal elements over which equation (9) applies. This approach unlike those techniques that shun the problem domain facilitates the numerical handling of nonlinearity and heterogeneity because medium parameters are considered homogeneous within an element though they may vary from element to element (piece-wise homogeneous).
We now make use of both the fundamental solution and the Green’s second identity to obtain a singular integral analog of equation (9):
(11)
in which is equal to the twice the nodal angle at the source node. . The problem domain is discretized with rectangular polygonal elements over which functional variables are interpolated. Using linear functions in space and time to interpolate functional quantities yields a general equation of the type
(12)
where represents the time component; represents current and previous times respectively, stand for linear interpolation functions in both time and space and is the number of nodes in the element. When equation (12) is substituted into equation (11) for all the functional variables, we obtain a system of discrete equations for a generic element of the problem domain.
(13)
where the coefficient matrices are given as:
(14a)
(14b)
(14c)
(14d)
2.3. Piecewise Approach
Hybridization of the boundary element method by the introduction of a FEM-like domain integration procedure makes it possible for domain properties to be treated as piecewise homogeneous. As a result, equation (9) is simplified to read:
(15)
The integral formulation of equation (15) is given as:
(16)
Quantities with overhead bars are element centroidal values of the media properties.
We hasten to comment that equation (16) is still nonlinear because media properties are still functions of the dependent variable even though they now have elemental representative averaged-values. As a result equation (13) converts to:
(17)
where
The temporal derivative is approximated by a 2-level time discretization, and equation (17) becomes:
(18)
The global matrix equation is obtained by assembling the element equation (18) and can be represented by:
(19)
where is the right hand side vector of known values which takes care of the boundary and initial conditions as well as point and distributed sources. The banded global matrix stores the coefficients of the element equations, while the mixed vector of unknown values contains both the dependent variable and its flux. Equation (19) is linearized by the Picard’s algorithm.
(20)
3. Numerical Illustrations
The formulations developed in section 2 are applied to study the behavior of Fisher’s equation.
3.1. Example 1
We consider equation (1) with the following values of coefficients . The general solution has a travelling wave solution of the form [33]:
(21)
And satisfies initial and boundary conditions with the solution travelling at a constant analytical wave speed: . In equation (1), energy released by the non-linear reaction term is balanced by the energy consumed by the diffusion process and results in travelling wave fronts. These have important applications in tracking population and gene spread. Figure 1 shows the profile of solving equation (1) with a unit reaction rate coefficient. In order to compare the results obtained herein with those of [33], the reaction coefficient is given a value of unity; and the boundary conditions given as The problem domain was truncated to and equally spaced 451 grid points were used to correspond with the problem parameters given in [33]. Table 1 shows a comparison of the theoretical wave speed (c=2.041241452) with the numerical wave speed:
(22)
Numerical Method | t=0 | t=4 | t=8 | t=12 | t=16 | t=20.0 |
Explicit | 2.0412414 | 1.8582093 | 1.8164059 | 1.8055922 | 1.8014246 | 1.7994060 |
Explicit predict Corrector | 2.0412414 | 1.9872378 | 2.0000463 | 2.0111028 | 2.0175552 | 2.0215470 |
Implicit | 2.0412414 | 2.2988649 | 2.4175923 | 2.4619361 | 2.4827240 | 2.4946464 |
Implicit Predict-Corrector | 2.0412414 | 2.0501223 | 2.0670761 | 2.0765175 | 2.0815662 | 2.0845510 |
Crank-Nicolson | 2.0412414 | 2.0552063 | 2.0766924 | 2.0882345 | 2.0943385 | 2.0979418 |
Explicit Mthd. Of Lines | 2.0412414 | 2.0021823 | 2.0404156 | 1.9685848 | 1.9650300 | 1.9632689 |
Fourth- Order Accurate Mthd. Of Lines | 2.0412414 | 2.0412854 | 2.0404156 | 2.0393077 | 2.0385086 | 2.0379799 |
GEM | 2.0412414 | 1.8486051 | 1.8448975 | 1.8496224 | 1.8553683 | 1.8581186 |
GEM Hermitian | 2.0412414 | 2.0348765 | 2.0465783 | 2.0784722 | 2.0975893 | 2.0997537 |
Finite element | 2.0412414 | 2.0258639 | 1.9953227 | 1.9781165 | 1.9696573 | 1.9654903 |
Fig. 1. Numerical solutions for example 1 .
Obtained by different numerical techniques. Table1 illustrates directly the accuracy of different computational techniques in relation to the numerical technique adopted herein as well as the effects of the discretization errors on the computed results for different times. It can be seen that the explicit schemes, underestimate the theoretical wave speed as time increases. This can be attributed to the first order truncation errors of the explicit scheme. It should be noted also the GEM without the Hermitian interpolation yielded results that are of the same level of accuracy as the explicit schemes. On the other hand, results for the fourth-order method of lines as well as the Hermitian GEM and the finite element method are found to be closer to the theoretical wave speed at steady state. We hasten to report that half the total number of grid points were used for the Hermitian GEM technique
3.2. Example 2
Consider the following form of the 1-D nonlinear Fisher’s equation:
(23a)
A pulse profile initial condition given as:
(23b)
and the boundary conditions are:
(23c)
The following parameters apply to this problem:
(23d)
The problem domain should be large enough to accommodate wave propagation.
Figures 2,3,4,5 illustrate the profiles of the solution history from short to long times. The interaction between the reaction and diffusion components at short times, is illustrated by figures 2 and 3. Initially diffusion dominates the reaction term and as result the peak of the graph goes down rapidly and becomes flat. This can be explained by noting that both the reaction and diffusion terms, are very small in the beginning and since diffusion for this problem, the absolute values of the diffusion term are found to dominate those of the reaction terms. However this trend immediately changes once the peak of the contour gets to its lowest level, at this stage, the reaction term starts to dominate the diffusion term (Fig.3) . The profile then continues with an upward ascent and with an increasing amplitude until it gets to the top where the value of the dependent variable is unity. The solution profiles then assume a bell shaped curve and become flatter with the lateral sides becoming very steep (Figs 4 and 5).
Fig. 3. Emerging solution profiles for example 2 .
Fig. 4. Emerging solution profiles for example 2 .
Fig. 5. Emerging solution profiles for example 2 .
3.3. Example 3
The Fisher’s equation is given by Mittal and Jain [34] as well as Cattani and Kudreyko [35] as:
(24a)
Equation (35a) is used mainly in population genetics and determines the a wave front of a particular gene though a population. While the second term on the right hand side of the equation determines how fast the infected particles are diffusing, the third term measures the infection rate. The following problem parameters apply:
. The initial condition is given as:
(24b)
The exact solution of equation (35a) has been determined by [35] as:
(24c)
The numerical results for a=0.5 and b=c= 1.0 and analytic solution are shown in Table 2 and are found to be in good agreement with those of [34] and [35].
Results at t=5 | ||
x-coordinates | Present | Analytical |
-30.00 | 0.5000000e+00 | 0.4999784e+00 |
-27.00 | 0.4999540e+00 | 0.4999487e+00 |
-26.00 | 0.4999342e+00 | 0.4999315e+00 |
-25.00 | 0.4999096e+00 | 0.4999086e+00 |
-23.50 | 0.4998586e+00 | 0.4998591e+00 |
-22.00 | 0.4997818e+00 | 0.4997828e+00 |
-13.75 | 0.4976467e+00 | 0.4976566e+00 |
-12.00 | 0.4961095e+00 | 0.4961253e+00 |
0.00 | 0.3953775e+00 | 0.3954037e+00 |
1.00 | 0.3676782e+00 | 0.3676518e+00 |
3.00 | 0.2978695e+00 | 0.2976762e+00 |
5.00 | 0.2147389e+00 | 0.2143460e+00 |
10.00 | 0.4821771e-01 | 0.4783219e-01 |
15.00 | 0.4642077e-02 | 0.4572265e-02 |
Results at t=6 | ||
30.00 | 0.5000000e+00 | 0.4999858e+00 |
-27.00 | 0.4999700e+00 | 0.4999662e+00 |
-26.00 | 0.4999568e+00 | 0.4999548e+00 |
-25.00 | 0.4999404e+00 | 0.4999397e+00 |
-23.50 | 0.4999068e+00 | 0.4999071e+00 |
-22.00 | 0.4998558e+00 | 0.4998568e+00 |
-13.75 | 0.4984455e+00 | 0.4984533e+00 |
-12.00 | 0.4974281e+00 | 0.4974405e+00 |
0.00 | 0.4269779e+00 | 0.4270202e+00 |
1.00 | 0.4061343e+00 | 0.4061379e+00 |
3.00 | 0.3501897e+00 | 0.3500467e+00 |
5.00 | 0.2756906e+00 | 0.2753178e+00 |
10.00 | 0.8241704e-01 | 0.8181168e-01 |
15.00 | 0.9709230e-02 | 0.9553112e-02 |
3.4. Example 4
In order to admit rigor and test to the formulation developed herein, we address the numerically challenging strong reaction problem by considering the case where . The presence of a strong reaction force, makes the solution to evolve into a wave like pattern. The solution technique must be able to handle the rapidly changing solution profile. This problem has been solved by many researchers [36,37,38]. The computational parameters used for this test follow that of [36] and are: . Figs 5a shows the numerical solution profiles for times: t=0.0005,0.0010,0.0015 and t= 0.003. The steeper profiles shown by Fig. 6 demonstrates the physical effects of increasing the reaction rate coefficient. In addition ,the flux profile is illustrated in Fig. 7. Both are in complete agreement with those of [36].
3.5. Example 5
The Hermitian GEM is further tested on a Fisher- type equation complicated by a nonlinear diffusion and heat capacity term. The governing equation can be described as :
(25a)
where is the heat capacity term. Equation (25a) has a lot of application in heat transfer process as demonstrated by Segal and Praagman [39] where they considered nonlinear heat conducting two dimensional bar with all the walls insulated except at the left boundary. A constant heat supply of serves as an input at this end. Both the thermal capacity and thermal conductivity are equal and are functions of temperature specified as: . Before going much further it should be observed that the adiabatic type of the boundary conditions applied to the top and bottom of the two-dimensional domain of this problem makes conversion to an equivalent 1-D analog straightforward. However we apply a 2-D non Hermitian GEM approach ( equations (9) to equation (20)) to enhance comparison with the 2-D example solved in literature. Segal and Praagman [39] solved this problem without the source term and as such it was possible for them to test their formulation with a closed form solution expressed as:
(25b)
The same approach has also been adopted for this test. The following problem parameters are employed: . The 2-D finite element method (FEM) results of [39] as well as those from the 2D averaged GEM and the 1-D Hermitian were found to be almost identical and compared very well with the analytical solution for that reason only those of the 1-D Hermitian GEM are displayed in Fig.8. The next step was to find the influence of the nonlinear source term on the numerical results using the 1D Hermitian GEM approach. At initial time, the value of the nonlinear diffusion coefficient is unity but the reaction term is quite small. As the dependent variable evolves with time the diffusion increases further and still becomes lager than the reaction term. This is contrary to our observation in example (2) and results in a continuous domination of diffusion over the reaction term over time as depicted in Fig. 9. Fig. 10 is a 3-D illustration of the solution profiles for the time of computation and clearly illustrates the nonlinear response of the of the governing equation to the specified boundary and initial conditions.
4. Conclusion
A modified boundary integral approach has been applied to nonlinear Fisher-type transient differential equations. The method reduces the governing equations to their weak forms before applying an element-based numerical approach to seek their numerical solutions. The overlying nonlinearity was handled straightforwardly by the Picard’s algorithm. It is noted that the introduction of an element-based hybrid numerica technique resulted in the presence of the source node in each computational element. This has a concomitant effect of not only converting a boundary based numerical method into its hybrid elemental analog but also facilitates the hybridization procedure. The benefits of this are immense in terms of not only simplifying the rigors of computation but at the same time advancing straightforwardly from one to multidimensional applications. The numerical results obtained so far confirms the reliability of the numerical formulation.
Acknowledgement
The author gratefully acknowledges the facilities and conducive environment provided by the African Institute of Mathematical Sciences (AIMS) in South Africa during his visit as senior researcher.
References