Galerkin-Bernstein Approximations of the System of Time Dependent Nonlinear Parabolic PDEs

The purpose of the research is to find the numerical solutions to the system of time dependent nonlinear parabolic partial differential equations (PDEs) utilizing the Modified Galerkin Weighted Residual Method (MGWRM) with the help of modified Bernstein polynomials. An approximate solution of the system has been assumed in accordance with the modified Bernstein polynomials. Thereafter, the modified Galerkin method has been applied to the system of nonlinear parabolic PDEs and has transformed the model into a time dependent ordinary differential equations system. Then the system has been converted into the recurrence equations by employing backward difference approximation. However, the iterative calculation is performed by using the Picard Iterative method. A few renowned problems are then solved to test the applicability and efficiency of our proposed scheme. The numerical solutions at different time levels are then displayed numerically in tabular form and graphically by figures. The comparative study is presented along with L 2 norm, and L ∞ norm.


Introduction
Reaction-diffusion systems have been extensively studied during the 20 th century.The study of the reaction-diffusion system reveals that different species have interactions with one another and that after these interactions, new species are created via chemical reactions.The solution of the reaction-diffusion system shows the chemical reaction's underlying mechanism and the various spatial patterns of the chemicals involved.Animal coats and skin coloration have been linked to reaction-diffusion processes, which have been considered to constitute a fundamental basis for processes associated with morphogenesis in biology.
There are numerous notable examples of coupled reaction-diffusion systems such as the Brusselator model, Glycolysis model, Schnackenberg model, Gray-Scott model, etc.With the help of the system size expansion, a stochastic Brusselator model has been suggested and investigated in the study cited in [1].The reaction-diffusion Brusselator model has been addressed by Wazwaz et al. through the decomposition technique [2].Because of its potential to provide a close analytical solution, the fractional-order Brusselator model was studied by Faiz et al [3].The Brusselator system stability of a reaction-diffusion cell as well as the Hopf bifurcation analysis of the system have been detailed by Alfifi [4].Qamar has analyzed the dynamics of the discrete-time Brusselator model with the help of the Euler forward and nonstandard difference schemes [5].The research article cited in [6] has been prepared by investigating the numerical analysis of the Glycolysis model using a well-known finite difference scheme.
Adel et al [7] have examined the synchronization problem of the Glycolysis reaction-diffusion model and designed a novel convenient control law.David et al [8] have analyzed the stability of turing patterns of the Schnackenberg model.Liu et al [9] have developed the bifurcation analysis of the aforementioned model.Khan et al. [10] have established a scheme for the solution of the fractional order Schnackenberg reaction-diffusion system.Numerical explorations have been applied to analyze the pattern formations of the model in the research article cited in [11].Gray and Scott [12] were the first to introduce the Gray-Scott model.They have proposed this model as an alternative to the autocatalytic model of Glycolysis [13].For this model, Pearson [14] has employed experimental studies to depict several sophisticated spot-type structures.Mazin et al. [15] have conducted an experiment using a computer simulation to investigate a range of far-from-equilibrium occurrences that emerge in a bistable Gray-Scott model.Many renowned authors [16,17] have evaluated the preceding model in which self-replicating structures have been noticed.McGough et al. [18] have conducted research on the bifurcation analysis of the patterns that are depicted in the model.In the research cited in [19], the linear stability and periodic stationary solutions of this model have been investigated.Some analytical results of this model have also been explored [20].Several prominent authors have studied the spatiotemporal chaos of the model in the research studies cited in [21] and [22].Furthermore, Wei [23] has analyzed the pattern formation of the two-dimensional Gray-Scott model.The model has also been explored by Kai et al. [24] using an innovative technique known as the second-order explicit implicit methodology.In recent years, the nonlinear Galerkin finite element approach has become increasingly prevalent as a means to investigate the model [25,26].Mach [27] has performed an in-depth examination of the quantitative evaluation of the model's numerical solution.
In references [28] and [29], the Gray-Scott reaction-diffusion system has been the subject of extensive wave modeling studies by eminent scholars.The simulation of the coupled model has been carried out by Owolabi et al. [30] using the higher-order Runge-Kutta method.The well-known Gray-Scott model's numerical findings have been calculated using the help of the hyperbolic B-spline [31].In order to analyze the ionic version of the model while it is being affected by an electric field, the Galerkin method has been deployed [32].With the use of the hybrid-asymptotic numerical method, Chen et al. [33] have investigated the model's dynamic behavior and stability.In the research study cited in [34], special polynomials have been employed to numerically solve the Gray-Scott model.Han et al. [35] have conducted an exhaustive investigation on the three-dimensional Gray-Scott model.In the process of assessing the model, the cubic B-spline has proven to be of considerable use by Mittal et al [36].
In the disciplines of engineering and mathematical physics, the Weighted Residual Method is an approximation method that can be leveraged to resolve problems.Analysis of structures, thermal expansion, stream of fluids, movement of masses, and the electromagnetic potential, etc. are examples of prominent problem fields of concern.Several distinct Weighted Residual Method variations are within our reach.The Galerkin Weighted Residual Method (also known as GWRM) has been put into practice for centuries, long before the invention of computers.It is generally agreed that this strategy is one of the best and most often used approaches available.Lewis and Ward have provided a comprehensive overview of the process in the article that is referenced in [37].This methodology has been effectively implemented in the well-known Black-Scholes model by Hossan et al. [38].Shirin et al. [39] have employed the Galerkin method in conjunction with other special polynomials to analyze the Fredholm equations.In the research referred to in [40], the approach was utilized to solve boundary value problems.In addition, this method has been used to perform a numerical calculation of the eigenvalues associated with the Sturm-Liouville problem [41].There have been several successful uses of this method for problems involving metal beams and polygonal ducts with rounded edges [42,43].
The objective of this study is to employ the modified Galerkin Weighted Residual Method in con-junction with the appropriate special polynomials to numerically evaluate the one-dimensional reactiondiffusion systems.Based on our best information, this study is presently unavailable.In addition to that, the study has provided the validation necessary to use the approach in one-dimensional reaction-diffusion systems.The main merit and advantage of the study are that by solving this type of system of equations, we will be able to analyze the behavior of the ecological system and forecast its future.The article is split up into four sections.Section 2 provides a detailed explanation of the formulation of our proposed method to solve the system of nonlinear parabolic partial differential equations.In the third section, the approach's implications are shown while analyzing the aforementioned system.Numerical and graphical representations are included here as well.The fourth section contains some concluding remarks and a general discussion.

Mathematical Formulation
Let us commence with the following system over the domain The boundary and initial conditions are as follows: and Let us assume the approximate solutions of System (2.1) be of the form where B j 's are the modified Bernstein polynomials and c j and d j are the coefficients dependent on time.
The first terms of the approximate solutions (2.4) have come from the boundary conditions of the system.The modified Bernstein polynomials are defined as follows: where U & L are the upper and lower limits of x.The last terms of Solution (2.4) will vanish at the boundary points.Therefore, the residual functions are Now we form the residual equations as: From the first residual equation, we can write Now we apply integration by parts in the above equation Then we substitute solution (2.4) in Equation (2.9).Therefore, the equation becomes, The first terms on both sides, and third terms on the left-hand side Equation (2.10) become zero because of boundary conditions.Therefore, the equation reduces to, The derivative and non-derivative terms of Equation (2.11) can be summarized via standard matrix notation as follows: where Here, K 1 and K 2 are n × n matrices, C 1 is n × n matrix, and F 1 is n × 1 matrix.The first two matrices K 1 and K 2 are called stiffness matrices.The other two matrices C 1 and F 1 are called forced matrix, and load vector respectively.Therefore, we apply the backward difference method on the first term of Equation (2.12) and rearrange the resulting terms as follows: Or, The second residual equation can be written as, After employing integration by parts and then substitution of (2.4) reduces the above equation, Since the first, and third terms on the left-hand side and the second term on the right-hand side of Equation (2.14) become zero, the equation reduces to, The derivative and non-derivative terms of Equation (2.15) can be summarized via standard matrix notation as follows: where Here, K 3 and K 4 are n × n matrices, C 2 is n × n matrix, and F 2 is n × 1 matrix.They are called stiffness matrices, forced matrices, and load vectors respectively.The application of the backward difference method on the first term of Equation (2.16) results in the following equation, By assembling Equations (2.13) and (2.17), we get the following recurrent system, To calculate the initial values of c j and d j , the initial conditions are set in Galerkin sense as follows, This process will help us to evaluate the numerical solutions of the nonlinear reaction-diffusion systems.

Numerical Examples and Applications
In this section, the previously described approach has been implemented into practice by solving a few examples of practical issues.Our methodology has been shown to be valid after being applied to the first test problem.The aforementioned procedure is then used, with a variety of parameters, to assess the subsequent test problems.The L 2 norm and L ∞ norm has been determined by the following expression, Where ∆t is the time increment and M ∆t is the approximate solution obtained using time increment ∆t.
Test Problem 1.Let us consider the system of the parabolic equations from the study of Manaa et.al. [44] where f(M, N) = M 2 N and x ∈ [a, b], t ⩾ 0. The boundary conditions and the initial conditions are considered as: Here, to obtain the numerical approximation, the effect of boundary conditions is insignificant because all terms of B j (x) are zero at the boundary points.We have employed the modified Galerkin method to the system of nonlinear partial differential equations (3.1) and therefore obtained the system of ordinary differential equations with respect to t.In this stage, we have used the α family of approximation in order to convert the system into recurrent relations and then we applied Picard iterative procedure.To find the initial guess of the given system, we have applied the weighted residual procedure on the initial conditions (3.2).
Tables (1) and ( 2) provide the numerical results of concentrations M(x, t) and N(x, t) for various values of x.For computation, we have taken ∆t = 0.1.The numerical approximations are derived at time levels t = 1 and t = 2. Throughout these tables, we have compared the results which we have obtained with the numerical approximations that have already been published in other well-known literature.The table demonstrates that our outcomes are reasonably comparable to those that have been published.It validates the accuracy of our approach to approximating the reaction-diffusion system numerically.The approximate results M(x, t) and N(x, t) of Equation (3.1) are presented in the following figure (1).In Figure (1) we have employed a three-dimensional graphical depiction of approximate solutions of M(x, t) and N(x, t) at different time levels for better understanding.The graphical representations agree with the results that we have obtained in the tables.Eventually, it makes sense clearly that the method is more applicable to solving such nonlinear parabolic PDE systems.In Figure (2), we have presented the error graph of M(x, t) and N(x, t) at time t = 10, where the absolute errors are computed between two different time increments, say ∆t = 0.2, ∆t = 0.4 and ∆t = 0.1, ∆t = 0.2.The L 2 norm and L ∞ norm, are presented in Table (3), which shows that the comparative errors are reduced significantly according to the reduction of the size of the time increments.Test Problem 2. The Gray-Scott Model is one of the most important models whose wave formations are similar to many waves formed in real life such as butterfly wings, gesticulation, damping, turning patterns, embryos, multiple spots, and so on [29,45].Let us consider the following model, where f(M, N) = MN 2 .The boundary conditions and the initial conditions are considered as follows: The domain of the model is [−50, 50].The values of the parameters are taken as ε 1 = 1, ε 2 = 0.01, p = 0.01, q = 0.12.Here for computational purposes, we have used 7 modified Bernstein polynomials.By applying the modified Galerkin method, we have used the backward difference method to transform the system of ordinary differential equations into the recurrent relations which is therefore solved by Picard iterative procedure.Numerical data of M(x, t) and N(x, t) of (3.1) are also presented in tabulated form in the following table at different time steps.The L 2 norm, and L ∞ norm, are presented in table (5), which shows that the comparative errors are reduced significantly according to the reduction of the size of the time increments.However, the order of convergences increased noticeably along with the reduction of the time length.

Conclusion
This research study has provided numerical approximations of nonlinear reaction-diffusion systems with specified boundary and initial conditions through the employment of the modified Galerkin method.To generate the trial solution, modified Bernstein Polynomials have been used.The simplification of the weighted residual leads to a system of ordinary differential equations which is then transformed into the recurrent relation by applying the backward difference formula.At this stage, we have used Picard's iterative procedure to approximate the trial solution.After successful derivation, we applied our proposed method to several models in order to test their applicability and effectiveness.We have solved and displayed the results both numerically and graphically.From those figures and numerical results, it is indisputable that our proposed method is an unconditionally stable, efficient, highly modular, and easily expandable method that can be applied to any type of system of nonlinear parabolic partial differential equations regardless of the type of the boundary conditions, type of non-linearity of the functions, coefficients are constants or function of independent variables.

Figure ( 3 )
Figure (3) is deployed to provide pictorial representations of the numerical concentrations M and N at different time levels.The results that are obtained in the table are shown graphically.The graphs are obtained for different time levels.The graphical presentation shows that the changes in concentrations are sufficiently small for different time levels.The L 2 norm, and L ∞ norm, are presented in table(5), which shows that the comparative errors are reduced significantly according to the reduction of the size of the time increments.However, the order of

Table 1 :
Numerical results of concentrations M(x, t) at different time levels with ∆t = 0.1 and first 7 modified Bernstein polynomials.

Table 2 :
Numerical results of concentrations N(x, t) at different time levels with ∆t = 0.1 and first 7 modified Bernstein polynomials.

Table 3 :
The L 2 and L ∞ norm at t = 10 for M(x, t) of equation (3.1).

Table 4 :
Numerical results of concentrations M(x, t) and N(x, t) ate different time levels with ∆t = 0.1 and first 7 modified Bernstein polynomials.

Table 5 :
The L 2 and L ∞ norm at t = 10 for M(x, t) of equation(3.4).