Numerical investigation of nonlinear extended Fisher-Kolmogorov equation via quintic trigonometric B-spline collocation technique

1.   Introduction

Nonlinear partial differential equations (PDEs) are extremely significant in real-world applications, which is why they play a significant part in the modeling of natural events. Most real-world problems can be characterized by nonlinear PDEs. However, analytical methods are unable to resolve some of these problems. Several numerical methods have been applied to solve these equations.

The Fisher Kolmogorov (F-K) models have been extensively utilized in the study of physical, material, and biological systems ^[1]. Some of these applications include, for instance, pattern formation ^[2], spatiotemperal chaos ^[3], traveling waves phenomena ^[4], liquid crystals domain walls propagation ^[5], reaction model for Alzheimer's disease ^[6], and tumor growth dynamics ^[7].

A cubic B-spline differential quadrature technique (DQM) was applied to simulate the fourth-order extended F-K equation in ^[8]. Error norms were used to assess the accuracy of the method. DQM and discretization were used to solve PDEs in ^[9]. Authors of ^[10,11] utilized differential quadrature and collocation methods based on quintic B-spline (QB) functions, i.e., QBDQM and QBCM, respectively, to simulate a numerically extended F-K (EFK) model. An analysis was conducted on stability, rate of convergence, and error norms. A numerical approach using trigonometric cubic B-spline functions was introduced to solve the time fractional gas dynamics equation was introduced in ^[12]. The generalized nonlinear time-fractional Klein-Gordon equation (TFKGE) was solved numerically using extended cubic B-spline (ECBS) functions in ^[13].

In ^[14], a trigonometric quintic B-spline method is employed to solve singularly perturbed boundary value problems, proving its ease and cost-effectiveness. A recent study used a trigonometric quintic B-spline method to solve the Korteweg-de Vries (KdV)-Kawahar equation, which has time-dependent parameters ^[15]. The numerical results obtained from this method were shown to be both effective and efficient. The trigonometric cubic and quintic B-splines are utilized to simulate the Burgers equation ^[16]. The accuracy of the outcome is evaluated by error norms and compared to the available approximate solutions. Results showed that the approach is exceptionally effective in resolving coupled Burgers' equation (cBE). A messless method (finite pointset) was applied in ^[17], to numerically solve the EFK equation. In ^[18], mesh-free numerical techniques were used based on DQM and radial basis functions (RBF) to simulate the F-K extended model. In ^[19], the authors used the Fourier-based meshless approach to solve a (3+1)-dimensional PDE. In ^[20], a Gaussian-cubic backward substitution method was utilized to solve nonlinear and linear problems in two dimensions and three dimensions with an irregular domain. The method's computational accuracy was effectively demonstrated.

The authors of ^[21] adopted element-free Galerkin interpolation to solve EFK numerically. In ^[22], the study on stochastic EFK was expanded by employing the numerical technique of Euler-Maruyama. In ^[23], He's variational method was adopted to study the solitons and period wave solution for EFK. This method simplifies calculation by reducing the order of the equation compared to previous methods. The authors of ^[24] proposed an approach to solve the EFK problem utilizing 2nd order splitting and orthogonal cubic spline collocation. Applying an $H^{1}$-Galerkin mixed finite element (H1-GMFE) method to the extended F-K equation via the splitting method, ^[25] derived optimal order error estimates. The authors of ^[26] utilized the C1-conforming finite element method to estimate error optimally for two-dimensional EFK equations (semi and wholly discrete). The authors of ^[27] analyzed the solution to the EFK equation in terms of its long-term behavior. In ^[28], the author presented a Crank-Nicolson finite difference scheme (CNFDS) to solve the EFK problem in two space dimensions with Dirichlet boundary conditions. In ^[29], the stability of the CNFDS in $L_{\infty}$ discrete norm for the EFK equation was investigated. The numerial solution to the EFK equation was derived in ^[30] using quartic B-spline based DQM (QAB-DQM). A novel meshless local collocation approach for numerically solving the 3D extended EFK equation was proposed in ^[31]. To discretize the time and spatial derivatives of the EFK equation, the second-order Crank-Nicolson scheme and meshless generalized finite difference method (GFDM) were used. Furthermore, in ^[32], a hybrid numerical method was introduced to address 2D nonlinear transient heat conduction concerns with temperature-dependent conductivity. They used Krylov deferred correction (KDC) for temporal discretization and GFDM for spatial discretization. In light of the aforementioned study, we propose a quintic trigonometric B-spline collocation technique to numerically simulate the extended F-K equation, which is expressed as:

subject to initial condition (IC)

and boundary condition (BC)

If $\gamma = 0$ in (1.1), we obtain the standard F-K equation, and its generalization is given as follows:

When approaching Lifshitz points, ^[33] phase transitions require the inclusion of the fourth-order derivative due to the necessity of including higher-order gradient factors in the energy-free functional. The EFK equation appears at the phase transition Lifshitz point, and the gradient systems evolution equation is given as follows:

where $F$ is the double-well potential and is given as follows:

If we choose $\beta = 1$, then we obtain the EFK equation.

For the structure of this paper: Section 2 includes the preliminaries about the quadratic trigonometric B-spline (QTB-spline) functions. In Section 3, we present the discretization of the extended F-K equation. Section 4 shows the stability analysis of the discretized system of the F-K equation using the von Neumann method. Section 5 shows the simulation of the F-K equation through numerical examples and compares the accuracy of the results. The conclusions are given in Section 6.

2.   QTB-spline basis functions

The problem domain $x\in[a, b]$ is uniformly partitioned into a mesh of length $h = \Delta x = \frac{{b - a}}{M}$, using the knots $x_i = a + ih$, $i = 0, 1, ..., M$, ensuring that $a = x_0 < x_1 < x_2 < \ldots < x_M = b$. Now, we specify the QTB-spline functions $\hat Q_i(x)$ for $i = -2, -1, 0, ..., M+$ as follows:

where $\rho(x_i) = sin (\frac{x-x_i}{2})$, $i = 0, 1, .., M$, and $\theta = sin(2.5h)sin(2h)sin(1.5h)sin(h)sin(0.5h)$.

The QTB-spline functions $\left\{ {QT_{-2}, QT_{-1}}, ..., QT_{M+1}, QT_{M+2} \right\}$ form a basis over the problem domain. Table 1 displays QTB-spline values and derivatives at knots.

Where

Now, we undertake that the estimation $u_M(x, t)$ to the function $u(x, t)$ at $(x, t_j)$ is represented by:

where, $C_i(t_j)$ are unknowns and need to be found using the initial and BCs, as well as the collocation technique. Each QTB-spline consists of six components, and each component is included in six QTB-splines. The function $u_M(x, t)$ can be represented as the variation over the component and stated as:

Using Eq (2.3), $u$, $u_x$, $u_{xx}$, $u_{xxx}$, and $u_{xxxx}$ the knots can be expressed as:

where $C_i^j = C_i(t_j)$.

3.   The problem's discretization

Currently, the time derivative of the problem (1.1) is discretized by the usual FDM, while $\theta$-weighted scheme is used for spatial derivatives as follows:

The Rubin-Graves ^[34] approach linearizes the term $u^3$ as follows:

Taking $\theta = \frac{1}{2}$ and using it in the above Eqs (3.1) and (3.2), we get

Simplifying the above equation and manipulating terms, we have

Let

Then, the above equation becomes

Now, using approximated $u$, $u_{xx}$, and $u_{xxxx}$ via the QTB-spline collocation technique, we get

We suppose that $\hat{A}^j_i = \tau_1 A^{j}_i - \frac{\Delta t}{2} \tau_6 + \frac{\gamma}{2} \Delta t \tau_{11}$, $\hat{B}^j_i = \tau_2 A^{j}_i - \frac{\Delta t}{2} \tau_7 + \frac{\gamma}{2} \Delta t \tau_{12}$,

$\hat{D}^j_i = \tau_3 A^{j}_i - \frac{\Delta t}{2} \tau_8 + \frac{\gamma}{2} \Delta t \tau_{13}$, $\hat{E}^j_i = \tau_1 B^{j}_i +\frac{\Delta t}{2} \tau_6 - \frac{\gamma}{2} \Delta t \tau_{11}$,

$\hat{F}^j_i = \tau_2 B^{j}_i + \frac{\Delta t}{2} \tau_7 - \frac{\gamma}{2} \Delta t \tau_{12}$, and $\hat{G}^j_i = \tau_3 B^{j}_i + \frac{\Delta t}{2} \tau_8 - \frac{\gamma}{2} \Delta t \tau_{13}$.

Then, the above equation becomes

The discretization of the BCs is as follows:

and

Solving Eqs (3.8)–(3.11), we have

where $\hat{\tau_m} = \tau_2 \tau_6 - \tau_1 \tau_7$, $\hat{\tau_p} = \tau_3 \tau_6 - \tau_1 \tau_8$, and $\hat{\tau_q} = \tau_3 \tau_7 - \tau_2 \tau_8$.

For $i = 0$, using Eqs (3.12) and (3.13) in (3.7) and manipulating terms, we get

For $i = 1$, using Eq (3.12) in (3.7), we get

For $i = M-1$, using Eq (3.13) in (3.7), we get

For $i = M$, using Eq (3.13) and (3.14) in (3.7), we get

At the time $t^j$, $j = 0, 1, ..., N$, and Eqs (3.16), (3.17), (3.7), (3.18), and (3.19) form a linear system of $(M+1)\times (M+1)$ order. We must establish the initial vectors $\left({C_0^0, C_1^0, ..., C_{M - 1}^0, C_M^0 } \right)$ from the IC to solve the system with the $M+1$ equation and $M+3$ unknowns. To remove the $C_{-1}^0$ and $C_{M+1}^0$, we use the IC $u(x, 0) = \phi(x)$ and its first and second derivatives at boundaries as follows:

Solving Eqs (3.20)–(3.24), we get

where $\tau_5 \tau_6 + \tau_4 \tau_7 = \eta_3$ and $\tau_5 \tau_6 - \tau_4 \tau_7 = \eta_4$.

For $i = 0$, using (3.25) and (3.26) in (3.20), we get

For $i = 1$, using (3.25) in (3.20), we get

For $i = M-1$, using (3.27) in (3.20), we get

For $i = M$, using (3.27) and (3.28) in (3.20), we get

Equations (3.25), (3.26), (3.20), (3.27), and (3.28) form a system $(M+1) \times (M+1)$ order at the time $t^0$.

4.   Analysis of stability

This section goes into stability analysis for the discretized system of the extended F-K equation via the von Neumann method ^[35]. According to the Duhamels' principle ^[36], the stability analysis of an inhomogeneous problem follows promptly from the homogeneous one. Thus, it seems necessary to examine the stability of the discretized system for the extended F-K equation with the force function $f = 0$. Taking $\theta = \frac{1}{2}$ and linearizing the nonlinear term $u^3$ by taking $u^2 = \hat{k}_1^2$ as a locally constant, the Eq (3.1) can be written as

The above equation can be written as

where

Now, using approximated values of $u$, $u_{xx}$, and $u_{xxxx}$ by the QTB-spline collocation technique in Eq (4.2), we get

where

We suppose one Fourier mode from the full solution $C_i^j = \delta ^j e^{\kappa i \vartheta}$ is used as trial solutions at $x_i$, where $\vartheta = \ \varepsilon h$. The $h$ is the element size, $\ \varepsilon$ is the mode number, and $\kappa = \sqrt{-1}$. Inverting this solution in Eq (4.3), we have

Simplifying and manipulating some terms, we have

where,

Inverting values of coefficients, we observe that $E^{*} \le \ A^{*}$, $F^{*} \le \ B^{*}$, $G^{*} \le \ D^{*}$, and, so, $\left| {\delta} \right| \le 1$. Therefore, the extended F-K equation discretized system is unconditionally stable.

5.   Computational results

Example 1. Consider extended F-K Eq (1.1) in $x \in[-4, 4]$ with the IC as

and the BCs as

Figures 1(a)–(c) illustrate the computations for $h = 0.1$, $\Delta t = 0.001$ at $t =$ 0, 0.05, 0.1, 0.15, and 2.0 for $\gamma = 0.0001$ and $\gamma = 0.001$, respectively. The figures indicate that solutions are identical for $\gamma = 0.0001$ and $\gamma = 0.001$, whereas solutions for $\gamma = 0.1$ rapidly decrease toward 0, confirming the extended F-K equation's stabilizing property. The solutions that are acquired efficiently reproduce the satisfactory qualitative properties of the extended F-K equation. Figures 2(a)–(c) depict the three-dimensional visualization of numerical solutions for different values of $\gamma$ ($\gamma = 0.0001$, $\gamma = 0.001$, and $\gamma = 0.1$) at $t = 0.2$, with $h = 0.1$ and $\Delta t = 0.001$. It can be noted from these figures that the outcome of solutions is nearly identical for very small values of $\gamma$. However, the deterioration of solutions to zero is very immediate in the case of $\gamma = 0.1$, which proves the stabilizing nature of the extended F-K equation. Table 2 presents a comparison between the current approach and existing methods in terms of $L_2$ and $L_\infty$ error norms. The comparison is made with a value of $\Delta t = 0.001$, $\gamma = 0.1$, and $t = 0.2$. At $M = 20$, our technique surpasses quartic B-spline differential quadrature method (QAB-DQM)^[30], QBCM^[11], and modified cubic B-spline based differential quadrature method (MCB-DQM)^[8] in terms of outcomes. Furthermore, our study shows that our results exhibit greater performance compared to the findings in QBDQM^[10] in terms of $L_2$ error norms. Therefore, we are able to conclude that the approach yields better results compared to certain methods described in the literature for a small grid size.

Example 2. Consider the extended F-K Eq (1.1) in $x \in [-4, 4]$ with the IC as

and the BCs as

Figure 3 illustrates the computations for $h = 0.1$, $\Delta t = 0.001$ at $t =$ 0.25, 1, 1.75, 2.5, 3.5, and 4.5 for $\gamma = 0.0001$. This figure shows that solutions decay and reach a stable state approaching the value 1 as time increases, which is the instant replicate of the adequate qualitative performance of the extended F-K equation. Figure 4 illustrates the 3D view of the computations with $\gamma = 0.0001$, $h = 0.025$, and $\Delta t = 0.0001$ for $t\in[0.25, 5]$. It is also obvious from this figure that solutions start to decay and reach a stable state, approaching the value 1 as time increases. Table 3 presents the $L_2$ and $L_\infty$ error norms, in addition to the $R_c$, for the parameters $\gamma = 0.0001$ and $\Delta t = 0.0001$ at $t = 1$ and $t = 4.5$. The table shows a decrease in error norms as the mesh size increases. Additionally, it is noted that error norms are at their minimum at time t = 4.5, indicating that the extended F-K equation is in a stable state. Furthermore, the accuracy of the projected procedure in space variables is second order.

Example 3. Consider the extended F-K Eq (1.1) in $x \in [-4, 4]$ with the IC as

and the BCs

Figure 5 illustrates the computations for $h = 0.1$, $\Delta t = 0.001$ at $t =$ 0.25, 1, 1.75, 2.5, 3.5, and 4.5 for $\gamma = 0.0001$. Figure 6 depicts the 3D view of the numerical solutions with $\gamma = 0.0001$, $\Delta t = 0.0001$, and $h = 0.025$ for $t\in[0.25, 5]$. It is obvious from these figures that solutions start decaying and reach a stable state approaching the value -1 as time increases, which is the instant replication of the adequate qualitative performance of the extended F-K equation.

Example 4. Lastly, we look at the nonhomogenous extended F-K equation

with

and BCs

where $f(u) = u^3-u$ and $g(x, t) = e^{-t}sin(2 \pi x)(e^{-2t}sin^2(2 \pi x)+4\pi^2+16\pi^4-2)$.

Table 4 shows $L_2$ and $L_\infty$ error norms as well as numerical convergence $R_c$ for $\gamma = 1$, $\Delta t = 0.01$ at $t = 0.5$ and $t = 1$. The table clearly indicates that error norms are negligible and decrease as mesh sizes expand. Additionally, the proposed method converges to a space variable of the second order. Figure 7 depicts exact and computation $u(x, t)$ for $\gamma = 1$, and $\Delta t = 0.001$ for $h = 0.025$ with $t =$ 0.1, 0.3, 0.5, and 0.7, while Figure 8 depicts 3D plots of exact and computation $u(x, t)$ with $\gamma = 1$, $\Delta t = 0.0005$ for $h = 0.05$ and $t\in [0.02]$, and $h = 0.025$ and $t\in [0, 1]$, respectively. These figures show excellent agreement between extact and computation solutions. The absolute error norms are also publicized in Figure 8, which are approximately in $10^{-3}$.

6.   Conclusions

A collocation technique based on QTB-spline functions is reported for homogeneous as well as nonhomogeneous extended F-K equations. The nonlinear term is handled by the Rubin-Graves (R-G) type linearization process. To validate results and check efficiency, three examples of homogeneous and one example of nonhomogeneous extended F-K equations are considered. The projected method has been found to yield improved results in comparison to the methods described in ^[8,10,11,30]. The performance of the projected technique and a relative investigation are accomplished graphically. Figures 1 and 2 portray that the nature of solutions is nearly similar for $\gamma = 0.0001$ and $\gamma = 0.001$ while, for $\gamma = 0.1$, solutions decompose rapidly to 0, ensuring the stability of the extended F-K equation. Figures 3 and 5 portray that solutions decay and attain a stable state approaching the values 1 and -1, respectively, as time increases, which is the instant replicate of the adequate qualitative performance of the extended F-K equation. Figure 8 shows excellent agreement between exact and numerical solutions, with absolute error norms of approximately $10^{-3}$. We show that the projected technique is unconditionally stable for the discretized extended F-K equations. Additionally, the technique is determined to be accurate to a second order in space. The numerical analysis proves the projected technique is straightforward and yields very accurate results.

Author contributions

All authors of this article have been contributed equally. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The authors extend their appreciation to the Deputyship for Research and Innovation, Ministry of Education in Saudi Arabia, for funding this research work through the project number ISP-2024.

Conflict of interest

The authors declare no conflict of interest.