# Analyzing Similarity Solution of Modified Fisher Equation

Authors: Esen Hanaç Duruk; Mehmet Emir Koksal; Ram Jiwari
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this paper, we first examine the type of structure of the solutions to the modified form of a nonlinear Fisher’s reaction-diffusion equation. The existence of the traveling wave solution to the equation in the long term is observed by using dynamical system theory and exhibiting a phase space analysis of its stable points. In parallel, we represent radial basis functions (RBFs)-based differential quadrature methods (DQMs) to close the solution of the equation. The stability analysis of the recommended method is demonstrated. Some initial-boundary value problems are considered test problems. The numerical results indicate extremely exact and stable initial and boundary conditions in the same domain with dissimilar time ranges.


---

## Body

## 1. Introduction

As one of the main fields of applied mathematics, differential equations arise in many areas of science and technology such as electromagnetic, fluid dynamics, mixing problems, population problems, wave motion, and in many other branches of engineering. Especially, partial differential equations can be used to describe a wide variety of phenomena in nature such as acoustics, electrodynamics, fluid flow, heat, and sound. Nonlinear partial differential equations are mostly renowned for describing the underlying behavior of nonlinear phenomena related to the nature of the real world [1–8].It is well-known that obtaining analytical solutions of nonlinear partial differential equations has a significant role in defining physical phenomena that are rising in several areas such as physics, biology, chemistry, and engineering, and there is a tremendous amount of work on the theory and techniques for solving partial differential equations both analytically [6–8] and numerically [9–14].Numerous nonlinear equations have numerous applications in fluid mechanics, chemical and plasma physics, and other fields. One of those equations is the Fisher equation:(1)ut=εuxx−θu1−u,which is a reaction-diffusion equation initially introduced by Fisher in the structure of population dynamics [15]. Especially, the Fisher equation is used as a powerful tool for modelling, analyzing, and solving real engineering problems and for discussing the results turned up at the end of analysis for the resolution of natural problems. It arises in varied physical phenomena, modelling an enormous diversity of problems, for instance, in ecology [16], combustion [17], nuclear reactor theory [18], and so forth. Many different methods, i.e., extended modified cubic B-spline algorithm [9], cubic trigonometric B-spline differential quadrature method [10], exponential modified cubic B-spline differential quadrature method [12], finite difference methods [19], B-spline methods [20], and extended cubic B-spline finite element method [21] have been applied to find the analytical and numerical solutions to the Fisher equation for many years.In Fisher equations, wave speeds determine the solutions. Fisher established that the solution of (1) evolves into a traveling wave with a wave speed greater than and equal to two whereas the solution is oscillatory when the wave speed is less than two which may cause nonphysical formation in some applications. Due to this situation, the desired wave movement may not be captured. However, finding the traveling wave for the smallest wave speed as possible is of great importance in physical and biological systems. For example, it is of vital importance to predict how long it will take for a cancer cell to spread to a certain area in the body, according to its rate of spread. So, that is the reason why we consider and demonstrate the existence of traveling wave solutions for a modified and complex form of equation (1), and this aspect constitutes the originality of the article.In this paper, the following dynamic system of a modified form of equation (1) is(2)ut=uxx−u1−u,−∞<x<∞,t>0,with initial condition (IC)(3)ux,0=fx,−∞<x<∞,and boundary conditions (BCs)(4)u−∞,t=u∞,t=0,have been considered. We specify analytical model (2) and analyze their properties using the phase plane approach.We suggest an organized process depending upon the presence of three different methods for equation (2) in order to obtain approximate traveling wave solutions. For this purpose, in Section 2, we first examine the dynamic system of the modified form of the Fisher equation by using phase space analysis to find a traveling wave solution occurrence that corresponds to a heteroclinic orbit in phase space. In Section 3, we demonstrate an outline introduction to the radial basis functions-based DQMs and so determine our numeric results depending on the radial basis functions-based DQMs using multiquadric in spatial direction for the modified reaction-diffusion equation (2). In Section 4, we present a numerical algorithm for the radial basis functions based on DQMs for the modified reaction-diffusion equation (2). In Section 5, we observe some test examples with different initial and boundary conditions to exhibit the accomplishment of the radial basis functions based on DQMs. In the last section, the determination of our work is given.


## 2. Stability of Travelling Wave Theory

In this section, we examine the traveling wave solutions of the modified Fisher equation (2) by using the transformation, namely,(5)s=x−ct,u=Us,where c is the wave speed. By substituting the transformation (5) into equation (2), then we obtain that(6)−cUs=Uss−U1−U.On writingUs=V, we obtain the dynamical system as follows:(7)Us=V,Vs=−cV+U1−U.The dynamical system (7) has been examined by a number of authors, including Murray [22]. The dynamical system (7) has two equilibrium points at P:0,0 and R:1,0. We require a monotone solution in 0≤U≤1 with Uss≤0. We next classify the equilibrium points by linearization. We first consider the equilibrium point P:0,0. The associated linear system is given by(8)Us=VVs=−cV−U,A+=01−1−c.Eigenvalues ofA+ and associated eigenvectors are given by(9)pw±=−c±c2−42,v±=1ρw±.In any condition ofc≥2, the point P:0,0 establishes a stable point situation in phase space. We next consider the equilibrium point R:1,0. The associated linear system is given by(10)Us=VVs=−cV−U,A=01−1/2−c.Eigenvalues ofA and associated eigenvectors are given by(11)p±=−c±c2−1/42,v±=1ρ±.Now, sincec≥1/2, the point R:1,0 is a stable point. As a result of the linearization theorem, point R = (1.0) is a stable point for a nonlinear system. We are looking for a heteroclinic occurrence for a specific wave speed c. Therefore, when we take c in different cases, namely.(i)
1/2≤c<2,(ii)
c≥2.In the first case while1/2≤c<2, the point P:0,0 indicates the spiral node while R:1,0 indicates a stable point. So, at this specific value of c, there is a heteroclinic orbit running from point P: (0,0) to R: (1,0). In the last case (iii), both points indicate stable nodes. We now demonstrate that the trajectories of system (7) are featured in the physical plane for c≥2 and c≥1/2 in time by using MATLAB ode45 package implementation in Figures 1 and 2, respectively.Figure 1
Trajectories of system (7) for a variety of c≥2. Each path symbolizes the trajectory for a variety of speeds c like 2 (blue), 3 (red), 4 (yellow), and 5 (purple).Figure 2
Trajectories of system (7) for a variety of c≥0.5. Each path symbolizes the trajectory for a variety of speeds c like 0.5 (blue), 2 (red), 3 (yellow), and 4 (purple).We note from above Figures that all paths for both casesc in the physical space are connected from P:0,0 to R:1,0 as time t in 0,50. Specifically, in Figure 2, when wave speed is chosen as c=0.5, even if there is a tiny oscillation, the blue line merges with other smooth lines while t tends to infinity. Hence, Figures 1 and 2 prove there is a connection from the point P:0,0 to R:1,0, which is called the “heteroclinic orbit” that represents waves.


## 3. RBFs-Based Differential Quadrature Methods

This section includes global and local RBFs-based DQMs. The methodologies of these methods are as follows:


### 3.1. Concise Description of DQM

The main task of DQM is to approximate the derivatives that occur in problems similar to conventional integral quadrature. Suppose we have a smooth functionux. Then, by DQM, the m-th order derivative of the function is approximated by the following formula:(12)∂mu∂xmxi≅∑j=1NWijmuxj,i=1,2,…,N.Here,N is the total number of nodes and Wijm are the unknown weighting coefficients (WCs) with respect to m-th order derivative of the function. The next step is to compute the WCs. In the literature, there are many schemes [23–28] available for computing the WCs. Due to connectivity between the nodes, these schemes are not easy to apply to the complex-shaped types of domains.To overcome the lack of these schemes, RBFs are used as test functions to compute the WCs. To complete this work, we used multiquadric (MQ) RBFs with constant shape parameters. Some well-known RBFs are listed in Table1.


### 3.2. Global RBFs-Based DQM

Herein, to compute the WCsWijm of the DQM mentioned in Section 3.2, we will determine MQ RBFs as test functions globally. That is why the method is called Global RBF-DQM. Except MQ RBFs, we can use other RBFs listed in Table 1. The MQ RBFs in 1D can be defined as follows:(13)φkr=r2+c2wherer=x−xk.Table 1
Definitions of some well-known RBFs.
RBFφr,r=x−xj,c>0MQc2+r2Inverse multiquadric (IMQ)1/c2+r2Gaussiane−cr2Thin plate spliner2logrMonomialr2k−1To explain the method, let us considerN the number of nodes in the computational domain a,b, i.e. a=x1<x2<⋯<xN=b. Now, according to equation (12) and the total number nodes, we have unknown N×N number of WCs. Now, put the MQ RBFs φkr,k=1,2,…,N from equation. (13) into equation (12); then, we have(14)φkmxi−xk=∑j=1NWijmφkxk−xj,i=1,2,…,N.Equation (14) is a system of N×N linear equations for each fixed i. The system can be written in matrix form:(15)BimWim=Gim,where(16)Bim=φ1x1−x1φ1x1−x2⋯φ1x1−xNφ2x2−x1φ2x2−x2⋯φ2x2−xN⋮⋮⋮⋮φNxN−x1φNxN−x2⋯φNxN−xN,Wim=Wi1mWi2m⋯WiNmT,Gim=φ1mxi−x1φ2mxi−x2⋯φNmxi−xNT.The above system is easy to solve. The solution of the above system will give us the WCs of first-order, second-order, and higher-order derivatives. The formula is as follows:(17)Wim=Bim−1Gim.For the large number of node points, the matrixBim becomes highly ill-conditioned and the inverse of the matrix becomes difficult. To overcome this situation, we used local domains from the global domain with a small set of nodes.


### 3.3. Local RBFs-Based DQM

In this method, we consider local support domains in place of the global domain. In pursuance of this consideration, equation (12) became as follows:(18)φkmxi−xk=∑j=1niWijmφkxk−xk,k=i1,i2,…,ini,i=1,2,…,N.Here,ni is the total number of nodes chosen in the local domain of the i-th node of the global domain and Wijm is the corresponding WCs of the m-th order derivative. The solution of equation (18) in matrix form is as follows:(19)Wnim=Bnim−1Gnim,where(20)Bnim=φi1xi1−xi1φi1xi1−xi2⋯φi1xi1−xiniφi2xi2−xi1φi2xi2−xi2⋯φi2xi2−xini⋮⋮⋮⋮φinixini−xi1φinixini−xi2⋯φinixini−xini,Wnim=wi1mwi2m...winimT,Gnim=φi1mxi−xi1φi2mxi−xi2...φinimxi−xiniT.For computational purpose, we have takenni=3 in the local domain.Theorem 1.
We consider the three uniform nodesxi−h,xi,xi+h, i.e. ni=3 in local RBF-based DQM. Then, the absolute errors in the approximation of first- and second-order derivatives are given as [29] follows:(i)
(21a)εx1≡h26u‴x1+h22c2u′x1.(ii)
(21b)ε^x1≡h212uIVx1+h2c2u‴x1−3h24c2u′x1.whereεx1≡u^′x1−u′x1 and ε^x1≡u^″x1−u″x1.


### 3.4. Shape Parameter (SP)

The SPc used in RBFs-based DQM has the utmost importance from the stability and accuracy point of view. Choosing an optimal SP c is a big challenge for the numerical simulation community and, till date, it is an open problem. There are various techniques available in the literature for choosing the optimal SP c. In the present work, we choose the trigonometry variable (TV) [30] technique to compute the SP c, which is given by(22)cj=smin+smax−smin×sinj,where j=1,2,…,N, Smin=1/N and Smax=3/N.


## 3.1. Concise Description of DQM

The main task of DQM is to approximate the derivatives that occur in problems similar to conventional integral quadrature. Suppose we have a smooth functionux. Then, by DQM, the m-th order derivative of the function is approximated by the following formula:(12)∂mu∂xmxi≅∑j=1NWijmuxj,i=1,2,…,N.Here,N is the total number of nodes and Wijm are the unknown weighting coefficients (WCs) with respect to m-th order derivative of the function. The next step is to compute the WCs. In the literature, there are many schemes [23–28] available for computing the WCs. Due to connectivity between the nodes, these schemes are not easy to apply to the complex-shaped types of domains.To overcome the lack of these schemes, RBFs are used as test functions to compute the WCs. To complete this work, we used multiquadric (MQ) RBFs with constant shape parameters. Some well-known RBFs are listed in Table1.


## 3.2. Global RBFs-Based DQM

Herein, to compute the WCsWijm of the DQM mentioned in Section 3.2, we will determine MQ RBFs as test functions globally. That is why the method is called Global RBF-DQM. Except MQ RBFs, we can use other RBFs listed in Table 1. The MQ RBFs in 1D can be defined as follows:(13)φkr=r2+c2wherer=x−xk.Table 1
Definitions of some well-known RBFs.
RBFφr,r=x−xj,c>0MQc2+r2Inverse multiquadric (IMQ)1/c2+r2Gaussiane−cr2Thin plate spliner2logrMonomialr2k−1To explain the method, let us considerN the number of nodes in the computational domain a,b, i.e. a=x1<x2<⋯<xN=b. Now, according to equation (12) and the total number nodes, we have unknown N×N number of WCs. Now, put the MQ RBFs φkr,k=1,2,…,N from equation. (13) into equation (12); then, we have(14)φkmxi−xk=∑j=1NWijmφkxk−xj,i=1,2,…,N.Equation (14) is a system of N×N linear equations for each fixed i. The system can be written in matrix form:(15)BimWim=Gim,where(16)Bim=φ1x1−x1φ1x1−x2⋯φ1x1−xNφ2x2−x1φ2x2−x2⋯φ2x2−xN⋮⋮⋮⋮φNxN−x1φNxN−x2⋯φNxN−xN,Wim=Wi1mWi2m⋯WiNmT,Gim=φ1mxi−x1φ2mxi−x2⋯φNmxi−xNT.The above system is easy to solve. The solution of the above system will give us the WCs of first-order, second-order, and higher-order derivatives. The formula is as follows:(17)Wim=Bim−1Gim.For the large number of node points, the matrixBim becomes highly ill-conditioned and the inverse of the matrix becomes difficult. To overcome this situation, we used local domains from the global domain with a small set of nodes.


## 3.3. Local RBFs-Based DQM

In this method, we consider local support domains in place of the global domain. In pursuance of this consideration, equation (12) became as follows:(18)φkmxi−xk=∑j=1niWijmφkxk−xk,k=i1,i2,…,ini,i=1,2,…,N.Here,ni is the total number of nodes chosen in the local domain of the i-th node of the global domain and Wijm is the corresponding WCs of the m-th order derivative. The solution of equation (18) in matrix form is as follows:(19)Wnim=Bnim−1Gnim,where(20)Bnim=φi1xi1−xi1φi1xi1−xi2⋯φi1xi1−xiniφi2xi2−xi1φi2xi2−xi2⋯φi2xi2−xini⋮⋮⋮⋮φinixini−xi1φinixini−xi2⋯φinixini−xini,Wnim=wi1mwi2m...winimT,Gnim=φi1mxi−xi1φi2mxi−xi2...φinimxi−xiniT.For computational purpose, we have takenni=3 in the local domain.Theorem 1.
We consider the three uniform nodesxi−h,xi,xi+h, i.e. ni=3 in local RBF-based DQM. Then, the absolute errors in the approximation of first- and second-order derivatives are given as [29] follows:(i)
(21a)εx1≡h26u‴x1+h22c2u′x1.(ii)
(21b)ε^x1≡h212uIVx1+h2c2u‴x1−3h24c2u′x1.whereεx1≡u^′x1−u′x1 and ε^x1≡u^″x1−u″x1.


## 3.4. Shape Parameter (SP)

The SPc used in RBFs-based DQM has the utmost importance from the stability and accuracy point of view. Choosing an optimal SP c is a big challenge for the numerical simulation community and, till date, it is an open problem. There are various techniques available in the literature for choosing the optimal SP c. In the present work, we choose the trigonometry variable (TV) [30] technique to compute the SP c, which is given by(22)cj=smin+smax−smin×sinj,where j=1,2,…,N, Smin=1/N and Smax=3/N.


## 4. Numerical Scheme for the Modified Reaction-Diffusion Equation

The main task of this section is to perform the simulation of the Fisher equation with constant coefficients. We approximate the partial derivatives of equation (2) with local RBFs-based DQM, and then equation (2) becomes a system of ODEs as follows:(23)duxi,tdt=∑j=1NWij2uxj,t−uxi,t1−uxi,t,xi∈Ω,i=1,2,…,N,with IC: uxi,0=fxi,−∞<x<∞. After applying the BCs, the system of ODEs (23) with IC is solved by the RK4 method.


### 4.1. Stability Analysis of the Scheme

In order to do the stability analysis of the proposed scheme, we use the following theorem:Theorem 2.
consider a nonlinear system (NLS)dU/dt=ΛU+ΨU and the corresponding linear system (LS) dU/dt=ΛU. Let 0,0 be a simple critical point (CP) of the NS and limU⟶0ΨU/U=0. If the CP (0,0) of the LS is asymptotically stable, then the CP 0,0 of NLS is also asymptotically stable.Proof.
For more details, see [31].
After using the BCs on system (23), the system be of the form(24)dUdt=ΛU+ΨU,where(25)Λ=W212W222⋯W2N−22W312W322⋯W3N−22⋮⋮⋮⋮WN−1,12WN−1,22⋯WN−1,N−12N−2×N−2,ΨU=ψ2Uψ3U⋯ψN−1UT,U=u2,u3,…,uN−1T,ψiU=uxi,t1−uxi,t,i=2,3,…,N−1.
Now, we consider the corresponding LS of the NLS (24) as follows:(26)dUdt=ΛU,
Then, the stability of the corresponding LS (26) implies the stability of the NLS (24). For more details, see [27]. Now, it can say that the NLS (24) is stable if the corresponding LS (26) is stable. We use the result that LS (26) is stable if all the Eigenvalues of Λ have a nonpositive real part [31]. During simulation, the computed eigenvalues of the matrix Λ have been plotted in Figure 3. The figure shows that the Eigenvalues of the matrix satisfied the above conditions. Hence, our algorithms are stable.Figure 3
Eigen value profiles forN=50 (a) and N=120 (b).
(a)(b)


## 4.1. Stability Analysis of the Scheme

In order to do the stability analysis of the proposed scheme, we use the following theorem:Theorem 2.
consider a nonlinear system (NLS)dU/dt=ΛU+ΨU and the corresponding linear system (LS) dU/dt=ΛU. Let 0,0 be a simple critical point (CP) of the NS and limU⟶0ΨU/U=0. If the CP (0,0) of the LS is asymptotically stable, then the CP 0,0 of NLS is also asymptotically stable.Proof.
For more details, see [31].
After using the BCs on system (23), the system be of the form(24)dUdt=ΛU+ΨU,where(25)Λ=W212W222⋯W2N−22W312W322⋯W3N−22⋮⋮⋮⋮WN−1,12WN−1,22⋯WN−1,N−12N−2×N−2,ΨU=ψ2Uψ3U⋯ψN−1UT,U=u2,u3,…,uN−1T,ψiU=uxi,t1−uxi,t,i=2,3,…,N−1.
Now, we consider the corresponding LS of the NLS (24) as follows:(26)dUdt=ΛU,
Then, the stability of the corresponding LS (26) implies the stability of the NLS (24). For more details, see [27]. Now, it can say that the NLS (24) is stable if the corresponding LS (26) is stable. We use the result that LS (26) is stable if all the Eigenvalues of Λ have a nonpositive real part [31]. During simulation, the computed eigenvalues of the matrix Λ have been plotted in Figure 3. The figure shows that the Eigenvalues of the matrix satisfied the above conditions. Hence, our algorithms are stable.Figure 3
Eigen value profiles forN=50 (a) and N=120 (b).
(a)(b)


## 5. Numerical Simulation and Discussions

In this section, some test problems are considered to check the efficiency and accuracy of the proposed RBFs based on DQM.Example 1.
We consider the problem (2) with the following IC and BCs over the domain [−10, 10](27)ux,0=sech210x,u−10,t=0=u10,t,for constructing and analyzing solutions of dynamical systems by using the phase space method. Firstly, we rewrite equation (2) into a system of first-order differential equations in time by introducing additional variables. The original variables and the new variables form a vector in phase space. The solution then becomes a curve in phase space parameterized by time. The curve is usually referred to as a trajectory or orbit. The vector differential equation is reformulated as a geometric description of the curve, as a differential equation only in terms of the phase space variables, without the original time parameterization. Finally, a solution in phase space is transformed back to the original environment.
The problem is simulated withN=120,Δt=0.0001, and the wave profiles of the problem are depicted in Figures 4 and 5 for different values of time. In the figures, we see that the wave profiles increase up to x=0 and then decrease after x=0. As the time increases, the wave profiles decrease.Figure 4
Wave profile of problem 1 at different time.Figure 5
Wave profile of problem 1 at different time.Example 2.
We consider the problem (2) with the following IC and BCs over the domain [−10, 10](28)ux,0=11+ex2,u−10,t=1=u10,t.
The exact solution of the problem is given by(29)ux,t=11+ex−5t2.
The problem is simulated withN=120,Δt=0.0001 and the wave profiles of the problem are depicted in Figures 6–8 for different values of time.
In Figures6–8, we determine the numerical solution of the initial-boundary value problem in Example 2 against x at times t=0,0.25,0.5,0.75,1, t=1,2,3,4,5 and t=6,7,8,9,10, respectively. The solution approaches the traveling wave quickly as t⟶∞ in all figures. Therefore, a traveling wave profile develops as t⟶∞ in the solution of the initial-boundary value problem in Example 2.Figure 6
Wave profile of problem 2 at different time.Figure 7
Wave profile of problem 2 at different time.Figure 8
Wave profile of problem 2 at different time.


## 6. Conclusion

In this paper, we first consider the occurrence of the traveling wave solutions to the modified Fisher equation (2) using phase plane analysis. We obtain a heteroclinic orbit that represents traveling wave solutions for c≥0.5, for which there is even a small existence of oscillation (see Figure 2). Afterward, we examine the traveling wave solutions to the modified Fisher equation (2), applying RBFs to them based on DQMs for dissimilar initial and boundary conditions. In the case of the trigonometric initial condition problem, where the solution develops a shock-like wavefront, good results have been obtained for various time ranges within the same distance domain. For the case of the exponential IC problem, the existence of a wavefront has been observed and obtained good results in same distance domain for dissimilar values of time. Accurate results have been established by using a numerical algorithm. The gained results with the RBFs based on DQMs are extremely exact and stable. Eventually, the stability analysis of RBFs based on DQM for the modified Fisher equation (2) has been represented. As a result, the modified Fisher equation (2) generates the identical traveling wave profile for all three approaches.Furthermore, the occurrence of a traveling wave is practically possible if the wave speed is at least 0.5 and greater.



---

*Source: 6806906-2022-06-29.xml*