An Advanced Galerkin Approach to Solve the Nonlinear Reaction-Diffusion Equations With Different Boundary Conditions

This study proposed a scheme originated from the Galerkin finite element method (GFEM) for solving nonlinear parabolic partial differential equations (PDEs) numerically with initial and different types of boundary conditions. The scheme is applied generally handling the nonlinear terms in a simple way and throwing over restrictive assumptions. The convergence and stability analysis of the method are derived. The error of the method is estimated. In the series, eminent problems are solved, such as Fisher’s equation, Newell-Whitehead-Segel equation, Burger’s equation, and Burgers-Huxley equation to demonstrate the validity, efficiency, accuracy, simplicity and applicability of this scheme. In each example, the comparison results are presented both numerically and graphically.


Introduction
In this paper, we consider a standard form of nonlinear parabolic partial differential equations (PDEs) in the following way defined on the spatio-temporal domain, R = Λ × (0, T ], T > 0 ∂A(x, t) ∂t = η ∂ 2 A(x, t) ∂x 2 + Ψ x, t, A(x, t), The initial condition is defined as subject to the boundary conditions where η > 0 is the diffusion coefficient, x is the space, and t, the time. The variable A is the function of time, t and space, x. The initial condition A(x, 0) = Θ describes the function of the spatial variable, x. The notation Λ is a region bounded in R j , typically j = {1, 2, 3}, with boundary ∂Λ ∈ C 2+δ , 0 < δ < 1. The function Ψ presents the reactionadvection map (either reaction or advection or both of them simultaneously). The reaction function is smooth enough, i.e. Ψ is continuous and at least first partial derivative exist. The advection term is defined as ∂A ∂x in (1). The notations α i ∈ R, β i ∈ R, γ i ∈ R (i = 1, 2) are constants and either α 1 = α 2 = 0 or, β 1 = β 2 = 0. It is noted that α 1 and β 1 both are not zero at the same time. Similarly α 2 and β 2 are not zero simultaneously. The left and right boundary points are 1 and 2 respectively and both 1 , 2 ∈ R. Also, in this study Ψ and Θ describe the nonlinear function of A(x, t) and x, respectively. It is also remarked that in many cases, Θ(x) can be linear steady state. The boundary conditions (3) cover a wide class depending on the parameters, for example: 3. Robin boundary conditions for all non-zero parameters, α 1 , α 2 , β 1 and β 2 .
The mixed boundary conditions are the next possibility based on proper parameters selection.
These types of nonlinear parabolic partial differential equations (PDEs) (1) frequently arise in the field of engineering and science. They are generally used to describe a wide variety of time-dependent phenomena. Some instantaneous applications are drawn in the following series (Mittal & Jain, 2012): • particle diffusion and the modeling of turbulence, • pricing of derivative investment instruments, • filtration of liquids, gas dynamics and elasticity, • modeling of biological species (organisms) with harvesting and Allee effect, and • chemical reactions, environmental pollution etc.
So the solutions of those types of equations have great significance in science and engineering, see, (Ozis, Aksan & Ozdes, 2003;Caldwell, Wanless & Cook, 1981;Reddy, 2014;Rao, 2017;Lewis & Ward, 1991;Chen, 2005;Akter et al., 2020) and references therein. The analytic solutions of second or higher order nonlinear PDEs are not easy and impossible in many cases except some simple cases. Many difficulties encountered, such as the nature of the governing equation, to construct an appropriate mesh, complex geometrical shapes, ill conditions, and singularity problems. So researchers devoted their attention to developing new and efficient solution methodologies to solve those problems numerically.
Instead of analytic solutions, the scientists introduced robust numerical techniques to find an approximation of the solutions, and in this study, we were interested in numerical solutions. Numerous authors had proposed different types of numerical methods for solving nonlinear parabolic PDEs (Hassan & Rashidi, 2014;Sarwar & Rashidi, 2016). Jacques solved this type of equations by predictor-corrector method (Jacques, 1983). Mittal and Jain solved nonlinear parabolic PDEs by Cubic B-splines collocation method (Mittal & Jain, 2012) which is limited only for Neumann boundary conditions. As a consequence, Caldwell et al. applied a finite element approach to solve Burger's but linear equation (Caldwell, Wanless & Cook, 1981). After transforming the nonlinear parabolic PDE into linear PDE by Hopf-Cole transformation, Ozis et al. found the solutions of the Burger's equation by finite element method with homogeneous Neumann boundary conditions (Ozis, Aksan & Ozdes, 2003). During the solution process, they used the kinematics viscosity within [0.004 − 0.01]. Using the classical Lie approach with additional generating conditions method, Cherniha and Dutka established the analytic solution of the Fisher's equation (Cherniha & Dutka, 2001). Ali et al. used quadratic element method for solving fractional-order differential equation (Ali, Kamrujjaman & Shirin, 2021). At the initial time, they compared the exact and the numerical solutions found by the finite element method with local perturbations of the analytic solution. In this method, they replaced the infinite interval into a finite interval so that the method should satisfy the homogeneous Neumann boundary conditions with sufficient exactness in this minimal interval. Ali et al. used Galerkin finite element method for solving nonlinear ordinary differential equations (Ali & Islam, 2017) . The FitzHugh-Nagumo equation and several other nonlinear reaction-diffusion were solved using the standard linear shape functions for various parameter values (Ali, Kamrujjaman & Islam, 2020;Lima et al., 2020Lima et al., & 2021. In this study, we have solved the nonlinear parabolic PDEs with various types of boundary conditions by our proposed scheme. In this method, different types of shape functions (linear, quadratic, cubic etc.) can be used instead of trial functions which are challenging to find that satisfy the boundary conditions; linear shape function is used in this paper with multiple nodes. In the literature, the current approach is not available as far as we know. So the main motive in this paper is to find the numerical solution of the nonlinear parabolic PDEs with different types of boundary conditions by our proposed scheme. This paper comprises five sections. In Section , a proper derivation of our proposed scheme for nonlinear parabolic PDEs is formulated. The convergence and stability analysis along with the iterative procedures are discussed in Section . The implication of the method is presented in Section considering four problems with strong nonlinearity. Error graphs and numerical results are also presented in this Section. General discussion and concluding remarks can be found in Section .
In the following section, we are moving to discuss the key idea of our scheme that is based on Galerkin finite element method to solve nonlinear parabolic PDEs (1).

Mathematical Formulation
First of all, we span the domain R into n sub-domains. Let S l be the l−th sub-domain of Λ with length h (= ∆x), which is the maximum length of all sub-domains and k (= ∆t) is the step length of t. An illustrative case as such is in Figure 1. If N is the total node number over the whole domain, then the total number of degrees of freedom over the domain is, N = (λ − 1) × n + 1, where n is the total number of elements and λ is the local node number per element. The elements are numbered similarly with parenthesis [e] (Ali & Islam, 2017) .
To find the solution, let us assume the trial solutioñ From the governing equation (1) which yields [e] ∂ã ∂t Where ψ i (i = 1, 2, 3, · · · , λ) are the shape functions. After substitutingã(x, t), ∂ã(x,t) ∂x from equation (4), we obtain The last equation is the matrix form of the system of ordinary differential equations with respect to t, where Equation (6) can be rewritten in a convenient matrix form as Here K and C are called the stiffness and forced matrices and is called the load vector. Since C is symmetric, it is also referred to as damping matrix.
To transform the ordinary differential equation (10) into a recurrent equation, lots of different methods are available in the literature. For good accuracy, it is imperative to choose an effective but generalized iterative method. So we choose the α-family of approximation instead of traditional finite difference method that comprises forward difference, backward difference, or central difference approaches.
In α-family of approximations, for different values of α, we find the following well-known numerical integration schemes (Reddy, 2014): By applying α-family of approximation, we form a system of algebraic equations to take a sequence of time steps of length h from time level j to j + 1.
Which can be written as Where the triangular matrices H & [H] are N × N and the vectors F j is an N × 1 in sizes, whose entries are as follows: f However, to treat the nonlinearity, we convert the nonlinear term as a multiple of A(x, t) to express in matrix form. Which we transform as a function of a j by using equation (4). Applying the initial condition into the recurrent equation (12), we iterate to find a better approximation for a particular time t. Putting back the values of a j for a particular element in equation (4), we find a piece-wise function of space and time for that element. Using the previously obtained function, we can find values of A(x, t) at any point on the element for a particular time t. Next Section contains one important branch of this study, the convergence and stability analysis which ensures that our method is convergent and our iterative procedure is stable for small time steps.

Convergence and Stability Analysis
Let the shape functions are ψ j (x) ∈ H 1 2 , j = 1, 2, · · · , λ where H 1 2 is the Hilbert space. Since η is a constant and Ψ x, t, A(x, t), ∂A(x,t) ∂x is continuous functions, so the solution of equation (1) uniquely exist (Babuska & Rheinboldt, 1978). Now substitute the trial function from equation (4) into (1), then we find equation (10) as Equation (16) is an initial value problem The integrating factor of equation (17) is e ζ j dt = e ζ j t .
Multiplying on both sides of the first equation of (17) by integrating factor and then integrating over (0, T ], we have and So in energy norm for a fixed t, we have Therefore, the series in equation (4) converges for a specific value of t.
Let D p,q be the space of (p − 1) times continuously differentiable functions onΛ, the closure of the domain, for which the restriction to Λ j , j = 0, 1, 2, · · · , n − 1 is a polynomial of degree at most λ − 1, where λ is the number of local nodes in each element. LetÃ ∈ D p,q be the finite element solution for a particular value of t, then for all υ ∈ D p,q where (., .) 0 is the inner product and B(., .) is the bilinear transformation defined in (Babuska & Rheinboldt, 1978;Babuska, 1972). Equation (20) provides an initial value problem instead of the system of ordinary differential equations. It is well known that for decreasing the increments in x, this result will converge to the solution of equation (10) for a particular t.
Now the iterative schemes can be applied in a relaxed way, if it is numerically stable. So write equation (12) as where is the amplification matrix. Here {a} j+1 , {a} j are the solution vector at time (t + 1) and t, respectively.
The solution {a} j+1 at time (t + 1) relies on the solution {a} j at time t. So error can grow with iteration. An iterative method is said to be stable if the error does not grow in an unbounded way with iteration. The necessary and sufficient conditions to bound the error within a borderline, the eigenvalue χ max of the amplification matrix [A] must be less or equal unity such that Equation (23) is an eigenvalue problem which is unconditionally stable if χ max is less or equal to unity for any time steps k (Reddy,2014). If χ max depends on the step size of time to be less or equal unity, then the procedure will be called conditionally stable. Thus α-family of approximation for α < 1 2 , all numerical schemes are stable for the time increment which satisfies the condition that is written as And for α ≥ 1 2 , the largest eigenvalue of the amplification matrix satisfies the following inequality which reveals that α-family of approximation is unconditionally stable if α ≥ 1 2 . In the following Section , we are going to present the solution of several nonlinear parabolic PDEs to demonstrate the efficiency of this method. Also, we will study the role played among different parametric values, exact-approximate solution, number of linear elements, and error terms while time varies.

Numerical Examples and Applications
In this section, the solutions of some re-known nonlinear parabolic equations with different boundary conditions are considered. All results are presented graphically and numerically, along with the exact solution.

Burger's Equation
Nonlinear Burger's equation is one of the most famous equations that frequently raised in fluid dynamics, magnetohydrodynamics, etc. Here we solve this by our proposed scheme with homogeneous Dirichlet boundary conditions (Wood, 2006) and the problem is where Re > 0 is the Reynolds number and σ > 1 is a parameter. The analytic solution of the problem (26) is Re cos(πx) .
For computation, we use only 40 linear elements with Re = 1, σ = 2, k = 0.01, h = 0.025, and α = 0.5. The numerical results are presented numerically in Table 3 (Appendix A). The graphical parallelism between the exact and approximate solutions are depicted in Figure 2 and in Figure 4. To clarify the understanding, let us consider the three dimensional surface plot, Figure 3 of the Burger's equation (26). Compared with left and right diagrams in Figure 3, it is clear that the method is more suitable to solve such nonlinear parabolic PDEs without any complexity.

Fisher's Equation
In this illustration, we consider the renowned Fisher's equation with non-homogeneous Dirichlet boundary conditions (Wazwaz & Gorguis, 2004;Kamrujjaman, Ahmed & Alam, 2019;. This equation arises in many applications of heat and mass transfer, biology, and ecology. Let us now consider the Fisher's equation The exact solution of equation (27) is given by We compute the numerical solution of this problem at different nodes of x and at different time levels by using only 40 linear elements with α = 1, k = 0.01, and h = 0.025. We reported the parallelism between exact and approximate solutions in Table 1. The graphical representation of exact and approximate solutions for different time levels are displayed in Figure 5 and in Figure 6. From Figure 6, we can observe that the graphical representation between the exact and approximate solution of (27).
Since the approximate solution is very close to the exact solution, so from these diagrams, it is not straightforward and easy to identify the difference between them. Both solutions are exhibited by two transparent graphs concerning time and space to overcome this difficult situation.
To clarify this, let us consider the absolute error map of (27) over space using various time levels. Eventually, it is apparent that this method is more applicable for solving such types of nonlinear parabolic PDEs with a perfect consent between the approximate and exact solution of (27); Figure 6. Now, we are devoted to demonstrating the famous Burgers-Huxley equation with homogeneous Neumann boundary conditions.

Burgers-Huxley Equation
Another renown nonlinear parabolic partial differential equation is the Burgers-Huxley equation with no-flux boundary conditions (Mittal & Jain, 2012;Zhou & Cheng, 2011) Here ϕ, , and τ are real parameters. It is easily verified that the exact solution of Burgers-Huxley equations (28) is where ϕ = −1, = 1, and τ = 2 is being used to get the numerical solution while time varies. We compute the numerical solution of the above system (28) Table 2. Graphical visualizations are shown in Figure 7. The three-dimensional solutions of (28) are presented in Figure 8 where there is a good correlation between the exact and imminent solutions, which ensured the acceptance of the numerical illustrations obtained by GFEM.
Numerical simulations depicted in Figure 7 indicate that the visible exact and approximate solutions have an excellent agreement. The error curves are presented in Figure 8 and are showing the higher accuracy, which ensures the stability of this algorithm. Finally, we illustrate the Newell-Whitehead-Segel equation with non-homogeneous Neumann boundary conditions.

Newell-Whitehead-Segel Equation
Let us consider the Newell-Whitehead-Segel equation that models the interaction of the impact of the diffusion term with the nonlinear effect of the reaction term (Salman Nourazar, Soori & Nazari-Golshan, 2015).
The exact solution of (29) is defined as  Table 4 (Appendix A). Additionally, the graphical visualization between the exact and approximate results is depicted in Figure 9. Visibly it is seen that the approximate and exact solutions coincide. For better understanding, we also affix the 3D solution maps (Figure 10) of (29) over the independent time and space variables, which depicts the acceptability of this method by displaying high accuracy with the exact solution.
Journal of Mathematics Research Vol. 14, No. 1; 2022 Figure 11. Absolute error of (29) at different time levels over the domain x ∈ [0, 1] The error curves, the absolute difference between the exact and numerical solutions of (29) are presented in Figure 11. The absolute error map provided a minimal error which is acceptable for numerically illustrated results.

Conclusion
In this study, an efficient scheme based on the Galerkin finite element method has been introduced to find the numerical solutions of a general class of nonlinear parabolic partial differential equations with different boundary conditions. The method was applied directly to find the approximate solutions of the equations without any limitations; it excluded linearization technique, no further assumption on trial functions, and no need to transform into lower-order PDE. To verify the method, we recalled several real-life problems, such as Burger's equation, Newell-Whitehead-Segel equation, and additional two renowned PDEs where all the equations carried out various boundary conditions. The approximate results were compared with the exact solutions and presented graphically (2D & 3D) and in tabular form. By the way of conclusion, it is apparent about the efficiency, accuracy, and acceptance of the proposed scheme because of its low cost, easily implementable, high accuracy with low compuatation time compared to the well-published literature.
This section contains some additional figures and tables to observe the accuracy of the solution methodology; the Galerkin Finite Element Method.