*Article* **A Cubic Spline Collocation Method to Solve a Nonlinear Space-Fractional Fisher's Equation and Its Stability Examination**

**Adel R. Hadhoud 1, Faisal E. Abd Alaal 2, Ayman A. Abdelaziz <sup>2</sup> and Taha Radwan 3,4,\***


**Abstract:** This article seeks to show a general framework of the cubic polynomial spline functions for developing a computational technique to solve the space-fractional Fisher's equation. The presented approach is demonstrated to be conditionally stable using the von Neumann technique. A numerical illustration is given to demonstrate the proposed algorithm's effectiveness. The novelty of the present work lies in the fact that the results suggest that the presented technique is accurate and convenient in solving such problems.

**Keywords:** Caputo sense; space-fractional Fisher's equation; cubic polynomial spline; von Neumann stability

### **1. Introduction**

Fractional calculus is a generalization of classical calculus that deals with fractional order differentiation and integration operations. Everyone who has studied elementary calculus is familiar with the differentiation operator, *D* = *d*/*dx*. Additionally, the *n* th derivative of a suitable function, *D<sup>n</sup> f*(*x*) = *d<sup>n</sup> f*(*x*)/*dnx*, is properly defined as long as *n* is a positive integer. In 1695, L'Hopital enquired to Leibniz about the meaning of *D<sup>n</sup> f* if *n* was fractional. Many prominent mathematicians, including Euler, Laplace, Fourier, Abel, Liouville, Riemann, and Laurent, have since studied time-fractional calculus. However, it was not until 1884 that the theory of generalized operators progressed to the point where it could be used as a starting point for modern mathematicians. Fractional calculus has been utilized in thermodynamics, viscoelasticity, bioengineering, aerodynamics, control theory, electromagnetics, finance, chemistry, and signal processing in recent decades [1–3]. Fractional calculus theory dates back to Leibniz in the sixteenth century, and many different types of fractional operators have been created and developed virtually as far as the classical theory has evolved. The Riemann–Liouville, Caputo, conformable, and Riesz operator approaches are among these [4–9]. Fractional differential equations (FDEs) appear in engineering processes, a number of physical phenomena, financial products, and biological systems, such as non-exponential relaxation patterns and anomalous diffusion [10]. Several numerical methods have been developed to obtain an approximate solution of FDEs. These methods include the Adomian decomposition method [11], finite difference method [12], variational iteration method [13], spectral methods [14], and homotopy perturbation method [15,16].

Fisher [17] introduced the nonlinear fractional Fisher equation as a mathematical model to describe the kinetic advancing rate of an advantageous gene. The fractional

**Citation:** Hadhoud, A.R.; Abd Alaal, F.E.; Abdelaziz, A.A.; Radwan, T. A Cubic Spline Collocation Method to Solve a Nonlinear Space-Fractional Fisher's Equation and Its Stability Examination. *Fractal Fract.* **2022**, *6*, 470. https://doi.org/10.3390/ fractalfract6090470

Academic Editors: Angelo B. Mingarelli, Leila Gholizadeh Zivlaei and Mohammad Dehghan

Received: 20 July 2022 Accepted: 22 August 2022 Published: 26 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Fisher's equation is also used in autocatalytic chemical reactions, nuclear reactor theory, neurophysiology, flame propagation, and branching Brownian motion processes [18–20]. In the open literature, many works describe solving the nonlinear fractional Fisher equation, see [20–23]. In our article, we propose a cubic polynomial spline-based technique for obtaining approximation solutions for the space-fractional Fisher's equation in the following form [24]:

$$\frac{\partial^a u(\mathbf{x},t)}{\partial \mathbf{x}^a} = \frac{\partial u(\mathbf{x},t)}{\partial t} - u(\mathbf{x},t)(1 - u(\mathbf{x},t)) + f(\mathbf{x},t), \mathbf{x} > 0,\tag{1}$$

subject to the conditions:

$$
\mu(a, t) = \mathbb{Z}\_1(t), \mu(b, t) = \mathbb{Z}\_2(t), 0 < t \le \mathcal{T}, \tag{2}
$$

and

$$
\mu(\mathbf{x},0) = \mathbf{g}(\mathbf{x}), a \le \mathbf{x} \le b,\tag{3}
$$

where *α* is a parameter in the Caputo sense that describes the fractional order of the space derivative, and 1 < *α* ≤ 2; if we take *α* = 2, Equation (1) reduces to the classical nonlinear Fisher's equation. Several scholars have studied and solved the space-fractional nonlinear Fisher's equation using various numerical techniques, for example, the generalized differential transform method (GDTM) and the variational iteration method (VIM) [25], the radial basis functions method [26], the quadratic polynomial spline-based method (QPSM) [24], and the Legendre spectral collocation method [27]. The above-mentioned articles show how to solve the space-fractional Fisher's equation in which the exact solution is unknown, but, as far as we are aware, no papers are using cubic polynomial spline for space-fractional Fisher's equation; hence, that is the motivation for our study.

**Definition 1.** *The Caputo space-fractional derivative operator of order α* > 0 *is given by [25]*:

$$\mathbb{E}\_{a}^{\mathbb{C}}D\_{\mathbf{x}}^{a}u(\mathbf{x},t) = \frac{\partial^{a}u(\mathbf{x},t)}{\partial \mathbf{x}^{a}} = \begin{cases} \frac{1}{\Gamma(n-a)} \int\_{a}^{\mathbf{x}} \frac{\partial^{n}u(\mathbf{r},t)}{\partial \mathbf{r}^{a}} (\mathbf{x}-\mathbf{r})^{n-a-1} d\mathbf{r}, & n-1 < a < n\_{r} \\\ \frac{\partial^{n}u(\mathbf{x},t)}{\partial \mathbf{x}^{a}}, & n = n \in \mathbb{N}\_{+} \end{cases}$$

*where* Γ( ) *is the gamma function. The right and left Caputo fractional derivatives of order α are given, respectively, by [28]*:

$$\begin{aligned} \, \_a^C D\_b^a u(\mathfrak{x}, t) &= \frac{1}{\Gamma(n - a)} \int\_{\mathfrak{x}}^b (-1)^n \frac{\partial^n u(\mathfrak{x}, t)}{\partial \mathfrak{x}^n} (\mathfrak{x} - \mathfrak{x})^{n - a - 1} d\mathfrak{x}, n - 1 < a \le n \, \_a \\\\ \_a^C D\_a^a u(\mathfrak{x}, t) &= \frac{1}{\Gamma(n - a)} \int\_{\mathfrak{x}}^{\mathfrak{x}} \frac{\partial^n u(\mathfrak{x}, t)}{\partial \mathfrak{x}^n} (\mathfrak{x} - \mathfrak{x})^{n - a - 1} d\mathfrak{x}, n - 1 < a \le n \, . \end{aligned}$$

The numerical solutions of fractional differential equations using numerical methods with spline functions and their use is an important area of research in numerical analysis, see [29–34]. Such problems appear very often in physics, engineering, medicine, biology, and other sciences. Spline functions can be integrated and differentiated due to being piecewise polynomials and, since they have a basis with small support, many of the integrals that occur in numerical methods are zero. Numerical methods with spline functions in obtaining the numerical solution of the differential equations lead to band matrices with a low computational cost. Cubic splines were used instead of using linear or quadratic splines, because they produce smooth approximations while maintaining a low computational cost.

#### **Definition 2.** *Polynomial Spline Function [35,36]:*

*Let a* = *x*<sup>0</sup> < ··· < *xn*−<sup>1</sup> < *xn* = *b be a subdivision of the interval* [*a*, *b*] *and m* ∈ *N. A function P* : [*a*, *b*] → *R is called a spline of degree m for this subdivision if P is* (*m* − 1) *times* *continuously differentiable on* [*a*, *b*] *and if the restriction of P to each subinterval* [*xi*, *xi*+1] *for i* = 0, 1, ··· , *n* − 1 *reduces to a polynomial of degree at most m. Thus, the polynomial spline function has the form:*

$$P(\boldsymbol{x}) = \begin{cases} \begin{array}{c} P\_0(\boldsymbol{x})\_\prime \\ P\_1(\boldsymbol{x})\_\prime \end{array} & \text{if } \boldsymbol{x} \in [\boldsymbol{x}\_{0\prime}, \boldsymbol{x}\_1], \\\ \begin{array}{c} \vdots \\ \vdots \\ \vdots \\ P\_{n-1}(\boldsymbol{x})\_\prime \end{array} & \begin{array}{c} \vdots \\ \vdots \\ \end{array} \\\ \begin{array}{c} \vdots \\ \end{array} \\\ \text{if } \boldsymbol{x} \in [\boldsymbol{x}\_{n-1\prime}, \boldsymbol{x}\_n]. \end{array} \end{cases}$$

*such that each polynomial spline segment is given by:*

$$P\_{\mathbf{i}}(\mathbf{x}) = \sum\_{r=0}^{m} a\_{\mathbf{i}r} (\mathbf{x} - \mathbf{x}\_r)^r,\\ \mathbf{i} = 0, \ 1, \ \cdot, \ \cdot, \ n - 1, \ \mathbf{x} \in \left[\mathbf{x}\_{\mathbf{i}}, \mathbf{x}\_{\mathbf{i}+1}\right].$$

The remainder of this article is structured as follows. In Section 2, a method is proposed that is based on the cubic spline polynomial method. In Section 3, the von Neumann method is used to discuss stability theoretically. In Section 4, we use a numerical example to demonstrate the suggested method's efficiency and accuracy. The conclusion is presented in Section 5.

#### **2. Derivation of the Method**

To approximate *u*(*x*, *t*) by a cubic polynomial spline-based method over finite elements, let the region *M* = [*a*, *b*] × [0, *T*] be discretized by a set of points *Mij*, which are the vertices of a grid of points  *xi*, *tj* where *xi* = *a* + *ih*, with *h* = Δ*x* = *<sup>b</sup>*−*<sup>a</sup> <sup>N</sup>* for each *i* = 0, 1, . . . , *N* and *tj* = *jk*, *k* = Δ*t* for *j* = 0, 1, . . . , *R*.

Let *W<sup>j</sup> <sup>i</sup>* be an approximation to *u xi*, *tj* obtained by the segment *Pi xi*, *tj* of the spline function passing through the points *xi*, *<sup>W</sup><sup>j</sup> i* and *xi*+1, *<sup>W</sup><sup>j</sup> i*+1 . Each segment has the form [37,38]:

$$P\_i(\mathbf{x}, t\_j) = a\_i(t\_j) \left(\mathbf{x} - \mathbf{x}\_i\right)^3 + b\_i(t\_j) \left(\mathbf{x} - \mathbf{x}\_i\right)^2 + c\_i(t\_j) \left(\mathbf{x} - \mathbf{x}\_i\right) + d\_i(t\_j), \tag{4}$$

for each *i* = 0, 1, ..., *N* − 1. To obtain expressions for the coefficients of Equation (4) in terms of *W<sup>j</sup> <sup>i</sup>* , *<sup>W</sup><sup>j</sup> <sup>i</sup>*+1, *<sup>S</sup><sup>j</sup> i* , and *S<sup>j</sup> <sup>i</sup>*+1, we first define:

$$P\_i(\mathbf{x}\_{i\prime} \, t\_j) = \mathcal{W}\_{i\prime}^j \, \tag{5}$$

$$P\_i(\mathbf{x}\_{i+1}, t\_j) = W\_{i+1'}^{j} \tag{6}$$

$$P\_{\mathbf{i}}^{(a)}\left(\mathbf{x}\_{i\prime},t\_{\mathbf{j}}\right) = \frac{\partial^{a}}{\partial \mathbf{x}^{a}} P\_{\mathbf{i}}\left(\mathbf{x}\_{i\prime},t\_{\mathbf{j}}\right) = S\_{\mathbf{i}\prime}^{\mathbf{j}} \tag{7}$$

$$P\_i^{(a)}\left(\mathbf{x}\_{i+1}, t\_j\right) = \frac{\partial^a}{\partial \mathbf{x}^a} P\_i\left(\mathbf{x}\_{i+1}, t\_j\right) = S\_{i+1'}^{j} \tag{8}$$

where 1 < *α* ≤ 2. By using (4), (5), and (6), we have:

$$d\_{\vec{i}} = \mathcal{W}\_{\vec{i}}^{\vec{j}}\,. \tag{9}$$

$$h^3 a\_i + h^2 b\_i + h c\_i + d\_i = W^{\bar{j}}\_{i+1'} \tag{10}$$

where *ai* ≡ *ai tj* , *bi* ≡ *bi tj* , *ci* ≡ *ci tj* , and *di* ≡ *di tj* . Using the right Caputo fractional derivative and (4) in (7), we obtain:

$$\begin{split} \frac{\partial^{a}}{\partial x^{\prime}} P\_{i}(\boldsymbol{x}\_{i}, \boldsymbol{t}\_{j}) &= \frac{1}{\Gamma(2-a)} \int\_{\substack{\boldsymbol{x}\_{i} \\ \boldsymbol{x}\_{i} \\ \Gamma(4-a)}}^{\boldsymbol{x}\_{i}+1} (\boldsymbol{\tau} - \boldsymbol{x}\_{i})^{1-a} P\_{i}^{(2)}(\boldsymbol{\tau}, \boldsymbol{t}\_{j}) d\boldsymbol{\tau} \\ &= \frac{6(2-a)\hbar^{3-a}}{\Gamma(4-a)} a\_{i} + \frac{2\hbar^{2-a}}{\Gamma(3-a)} b\_{i} = S\_{i'}^{\boldsymbol{j}} \end{split} \tag{11}$$

and from using the left Caputo fractional derivative and (4) in (8), we obtain:

$$\begin{split} \frac{\partial^{\boldsymbol{a}}}{\partial \mathbf{x}^{\boldsymbol{a}}} P\_{\boldsymbol{i}}(\mathbf{x}\_{i+1}, \boldsymbol{t}\_{\boldsymbol{j}}) &= \frac{1}{\Gamma(2-\boldsymbol{a})} \int\_{\begin{subarray}{c} \mathbf{x}\_{i} \\ \mathbf{x}\_{i} \end{subarray}}^{\mathbf{x}\_{i+1}} (\mathbf{x}\_{i+1} - \boldsymbol{\tau})^{1-\boldsymbol{a}} P\_{\boldsymbol{i}}^{(2)}(\boldsymbol{\tau}, \boldsymbol{t}\_{\boldsymbol{j}}) d\boldsymbol{\tau} \\ &= \frac{6h^{3-\boldsymbol{a}}}{\Gamma(4-\boldsymbol{a})} a\_{\boldsymbol{i}} + \frac{2h^{2-\boldsymbol{a}}}{\Gamma(3-\boldsymbol{a})} b\_{\boldsymbol{i}} = S\_{\boldsymbol{i}+1'}^{\boldsymbol{j}} \end{split} \tag{12}$$

We obtain the following expressions by solving (9)–(12):

$$a\_{i} = \frac{\Gamma(4-a)}{6(1-a)h^{3-a}} \left( S\_{i}^{j} - S\_{i+1}^{j} \right),$$

$$b\_{i} = \frac{\Gamma(3-a)}{2(1-a)h^{2-a}} \left( (2-a)S\_{i+1}^{j} - S\_{i}^{j} \right),$$

$$c\_{i} = \frac{1}{h} \left( W\_{i+1}^{j} - W\_{i}^{j} \right) + \frac{\Gamma(4-a)}{6(1-a)h^{1-a}} \left( S\_{i+1}^{j} - S\_{i}^{j} \right) + \frac{\Gamma(3-a)}{2(1-a)h^{1-a}} \left( S\_{i}^{j} - (2-a)S\_{i+1}^{j} \right),$$

$$d\_{i} = W\_{i}^{j}. \tag{13}$$

Using the first derivative's continuity condition at *<sup>x</sup>* <sup>=</sup> *xi*, that is *<sup>P</sup>*(1) *i xi*, *tj* = *P*(1) *i*−1 *xi*, *tj* , we obtain:

$$c\_i - c\_{i-1} = 3h^2 a\_{i-1} + 2hb\_{i-1} \tag{14}$$

Substituting the expressions *ai*, *bi* , and *ci* from Equation (13) into Equation (14), we obtain:

$$\begin{array}{l} \mathcal{W}\_{i+1}^{j} - 2\mathcal{W}\_{i}^{j} + \mathcal{W}\_{i-1}^{j} + \frac{\Gamma(4-a)}{6(1-a)h^{-a}} \left( S\_{i+1}^{j} - S\_{i}^{j} \right) + \frac{\Gamma(3-a)}{2(1-a)h^{-a}} \left( S\_{i}^{j} - (2-a)S\_{i+1}^{j} \right) \\ \quad + \frac{\Gamma(4-a)}{6(1-a)h^{-a}} \left( S\_{i-1}^{j} - S\_{i}^{j} \right) + \frac{\Gamma(3-a)}{2(1-a)h^{-a}} \left( (2-a)S\_{i}^{j} - S\_{i-1}^{j} \right) \\ \quad = \frac{\Gamma(4-a)}{2(1-a)h^{-a}} \left( S\_{i-1}^{j} - S\_{i}^{j} \right) + \frac{\Gamma(3-a)}{(1-a)h^{-a}} \left( (2-a)S\_{i}^{j} - S\_{i-1}^{j} \right). \end{array} \tag{15}$$

Equation (15) can be rewritten in the form:

$$=\begin{pmatrix} \mathcal{W}\_{i+1}^{j} - 2\mathcal{W}\_{i}^{j} + \mathcal{W}\_{i-1}^{j} \\ \frac{\Gamma(4-a)}{2(1-a)h^{-a}} - \frac{\Gamma(3-a)}{(1-a)h^{-a}} - \frac{\Gamma(4-a)}{6(1-a)h^{-a}} + \frac{\Gamma(3-a)}{2(1-a)h^{-a}} \end{pmatrix} \mathcal{S}\_{i-1}^{j} \\ \quad + \left( \frac{2\Gamma(4-a)}{6(1-a)h^{-a}} - \frac{\Gamma(3-a)}{2(1-a)h^{-a}} - \frac{(2-a)\Gamma(3-a)}{2(1-a)h^{-a}} \right) \\ \quad - \frac{\Gamma(4-a)}{2(1-a)h^{-a}} + \frac{(2-a)\Gamma(3-a)}{(1-a)h^{-a}} \right) \mathcal{S}\_{i}^{j} \\ \quad + \left( \frac{(2-a)\Gamma(3-a)}{2(1-a)h^{-a}} - \frac{\Gamma(4-a)}{6(1-a)h^{-a}} \right) \mathcal{S}\_{i+1}^{j} \end{pmatrix} \tag{16}$$

After slight rearrangements, Equation (16) becomes:

$$\mathcal{W}\_{i-1}^{j} - 2\mathcal{W}\_{i}^{j} + \mathcal{W}\_{i+1}^{j} = \beta \mathcal{S}\_{i-1}^{j} + \gamma \mathcal{S}\_{i}^{j} + \beta \mathcal{S}\_{i+1}^{j}, i = 1, \dots, N - 1,\tag{17}$$

where *<sup>β</sup>* <sup>=</sup> (2*α*−3)*hα*Γ(3−*α*) <sup>6</sup>(*α*−1) and *<sup>γ</sup>* <sup>=</sup> *<sup>α</sup>hα*Γ(3−*α*) <sup>3</sup>(*α*−1) . As *<sup>α</sup>* <sup>→</sup> 2, system (17) reduces to:

$$\mathcal{W}\_{i-1}^{j} - 2\mathcal{W}\_{i}^{j} + \mathcal{W}\_{i+1}^{j} = \frac{1}{6} \left( S\_{i-1}^{j} + 4S\_{i}^{j} + S\_{i+1}^{j} \right), i = 1, \ldots, N - 1. \tag{18}$$

**Remark 1:** *For Equation (17), the local truncation error is:*

$$T\_i^{\*j} = \beta \left( u\_{i-1}^j + u\_{i+1}^j \right) - 2\gamma u\_i^j - \beta \left( D\_x^2 u\_{i-1}^j + D\_x^2 u\_{i+1}^j \right) - \gamma D\_x^2 u\_{i\prime}^j \tag{19}$$

and can be obtained by expanding this equation in terms of *u xi*, *tj* and its derivatives in the Taylor series, as follows:

$$T\_i^{\*j} = \left(h^2 - (\gamma + 2\beta)\right) D\_x^2 u\_i^j + \left(\frac{h^4}{12} - \beta h^2\right) D\_x^4 u\_i^j + \left(\frac{h^6}{360} - \frac{\beta h^4}{12}\right) D\_x^6 u\_i^j + \dots \tag{20}$$

Then, the last expression *T*∗*<sup>j</sup> <sup>i</sup>* can be written as:

$$T\_i^{\*j} = h^a \left( h^{2-a} - \theta \right) D\_x^2 u\_i^j + h^{2+a} \left( \frac{h^{2-a}}{12} - \delta \right) D\_x^4 u\_i^j + h^{4+a} \left( \frac{h^{2-a}}{360} - \frac{\delta}{12} \right) D\_x^6 u\_i^j + \dots \tag{21}$$

where *<sup>θ</sup>* <sup>=</sup> <sup>Γ</sup>(<sup>3</sup> <sup>−</sup> *<sup>α</sup>*) and *<sup>δ</sup>* <sup>=</sup> (2*α*−3)Γ(3−*α*) <sup>6</sup>(*α*−1) . From the expression (21) of the local truncation error, our scheme (17) is reduced to *<sup>O</sup>*(*hα*), 1 <sup>&</sup>lt; *<sup>α</sup>* <sup>≤</sup> 2.

The fractional Fisher's Equation (1) can now be defined as follows in terms of *S<sup>j</sup> i* :

$$S\_i^j = \frac{\partial^\alpha W\_i^j}{\partial x^\alpha} = \frac{\partial W\_i^j}{\partial t} - \mathcal{W}\_i^j \left(1 - \mathcal{W}\_i^j\right) + f\left(\mathbf{x}\_i, t\_j\right).$$

Using the finite difference method, we obtain:

$$\frac{\partial \mathcal{W}\_i^j}{\partial t} \approx \frac{\mathcal{W}\_i^j - \mathcal{W}\_i^{j-1}}{k}.$$

This allows us to express *S<sup>j</sup> <sup>i</sup>*−1, *<sup>S</sup><sup>j</sup> i* , and *S<sup>j</sup> <sup>i</sup>*+<sup>1</sup> as follows:

$$\begin{array}{lcr} S\_{i-1}^{j} = \frac{W\_{i-1}^{j} - W\_{i-1}^{j-1}}{k} - \sigma\_{i-1}^{j} W\_{i-1}^{j} + f\_{i-1}^{j} \\ S\_{i}^{j} = \frac{W\_{i}^{j} - W\_{i}^{j-1}}{k} - \sigma\_{i}^{j} W\_{i}^{j} + f\_{i}^{j} \\ S\_{i+1}^{j} = \frac{W\_{i+1}^{j} - W\_{i+1}^{j-1}}{k} - \sigma\_{i+1}^{j} W\_{i+1}^{j} + f\_{i+1}^{j} \end{array} \tag{22}$$

where *σ<sup>j</sup> <sup>i</sup>* = <sup>1</sup> <sup>−</sup> *<sup>W</sup><sup>j</sup> i* . Substituting (22) into (17) gives us the following system:

$$\begin{aligned} A\_i \mathcal{W}\_{i-1}^j + B\_i \mathcal{W}\_i^j + \mathbb{C}\_i \mathcal{W}\_{i+1}^j &= A\_i^\* \mathcal{W}\_{i-1}^{j-1} + B\_i^\* \mathcal{W}\_i^j + \mathbb{C}\_i^\* \mathcal{W}\_{i+1}^{j-1} + \rho\_i^j, \\ \quad i = 1, \dots, N-1, \ j = 1, \dots, R, \end{aligned} \tag{23}$$

where:

$$A\_i = k - \beta + \beta k \sigma\_{i-1'}^{j} A\_i^\* = -\beta,$$

$$B\_i = -2k - \gamma + \gamma k \sigma\_{i'}^{j} B\_i^\* = -\gamma,$$

$$\mathbf{C}\_i = k - \beta + \beta k \sigma\_{i+1'}^{j} \mathbf{C}\_i^\* = -\beta,$$

$$i = 1, \ldots, i \qquad \mathbf{j} \quad \mathbf{j} \quad \mathbf{j} \quad \mathbf{j} \quad \mathbf{j}.$$

and:

$$
\rho\_i^j = \beta k f\_{i-1}^j + \gamma k f\_i^j + \beta k f\_{i+1}^j.
$$

$$
\text{satisfies of } N - 1 \text{ equations in the unknowns } W\_i, i
$$

System (23) consists of *N* − 1 equations in the unknowns *Wi*, *i* = 0, 1, ... , *N*. Two additional equations are required to solve this system. The boundary conditions in (2) are used to derive these equations, which can be represented as:

$$\mathcal{W}\_0^j = \mathcal{Q}\_1(\mathbf{t}\_j), \mathcal{W}\_N^j = \mathcal{Q}\_2(\mathbf{t}\_j), j = 0, 1, \dots, R. \tag{24}$$

Writing (23) and (24) in matrix form gives:

$$
\Phi \mathsf{W}^{j} = \Phi^{\*} \mathsf{W}^{j-1} + r^{j}, \tag{25}
$$

where:

$$W^j = \left(W\_0^j, W\_1^j, W\_2^j, \dots, \dots, \dots, \dots, W\_{N-1}^j, W\_N^j\right)^T,$$

$$\Phi = \begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0\\ A\_1 & B\_1 & C\_1 & 0 & \cdots & 0 & 0 & 0\\ 0 & A\_2 & B\_2 & C\_2 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & A\_{N-1} & B\_{N-1} & C\_{N-1}\\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 \end{bmatrix},$$

$$\Phi^\* = \begin{bmatrix} 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0\\ A\_1^\* & B\_1^\* & C\_1^\* & 0 & \cdots & 0 & 0 & 0\\ 0 & A\_2^\* & B\_2^\* & C\_2^\* & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & \cdots & A\_{N-1}^\* & B\_{N-1}^\* & C\_{N-1}^\*\\ 0 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \end{bmatrix}'$$

and,

$$\boldsymbol{r}^{\boldsymbol{j}} = \left(\boldsymbol{\mathcal{Q}}\_1(\boldsymbol{t}\_{\boldsymbol{j}}), \boldsymbol{\rho}\_{1'}^{\boldsymbol{j}}, \boldsymbol{\rho}\_{2'}^{\boldsymbol{j}}, \dots, \dots, \dots, \boldsymbol{\rho}\_{N-1'}^{\boldsymbol{j}} \otimes\_2(\boldsymbol{t}\_{\boldsymbol{j}})\right)^T$$

where <sup>Φ</sup><sup>∗</sup> and <sup>Φ</sup> are (*<sup>N</sup>* <sup>+</sup> <sup>1</sup>) <sup>×</sup> (*<sup>N</sup>* <sup>+</sup> <sup>1</sup>) matrices and *<sup>W</sup><sup>j</sup>* and *<sup>r</sup>*<sup>j</sup> are (*<sup>N</sup>* <sup>+</sup> <sup>1</sup>) vectors.

Now, the initial condition, *<sup>u</sup>*(*x*, *<sup>t</sup>*0) <sup>=</sup> *<sup>g</sup>*(*x*) for <sup>a</sup> <sup>≤</sup> *<sup>x</sup>* <sup>≤</sup> b, implies that *<sup>W</sup>*<sup>0</sup> *<sup>i</sup>* for *i* = 0, 1, . . . , *N*. These values can be utilized in (25) to obtain the value of *W*<sup>1</sup> *<sup>i</sup>* for*i* = 0, 1, . . . , *N*. If the technique is reapplied once all the approximations *W*<sup>1</sup> *<sup>i</sup>* are determined, the values of *W*<sup>2</sup> *<sup>i</sup>* , *<sup>W</sup>*<sup>3</sup> *<sup>i</sup>* ,... , *<sup>W</sup><sup>R</sup> <sup>i</sup>* can be derived in the same way.

The nonlinear term in the system is linearized using the technique in (25). For example, if *j* = 1, we approximate *σ*<sup>1</sup> *<sup>i</sup>* by *<sup>σ</sup>*1<sup>∗</sup> *<sup>i</sup>* computed from *<sup>w</sup>*<sup>0</sup> *<sup>i</sup>* , which obtains an approximation to *w*<sup>1</sup> *<sup>i</sup>* . For *<sup>j</sup>* = *<sup>m</sup>*, we approximate *<sup>σ</sup><sup>m</sup> <sup>i</sup>* by *<sup>σ</sup>m*<sup>∗</sup> *<sup>i</sup>* computed from *<sup>w</sup>m*−<sup>1</sup> *<sup>i</sup>* and obtain an approximation to *w<sup>m</sup> i* .

#### **3. Stability Analysis**

The von Neumann technique will be used to investigate our system's stability (23). To do this, we must linearize the nonlinear term *u*(*x*, *t*)(1 − *u*(*x*, *t*)) of Fisher's Equation (1), by assuming that the corresponding quantities *σ<sup>j</sup> <sup>i</sup>*−1, *<sup>σ</sup><sup>j</sup> <sup>i</sup>* , and *<sup>σ</sup><sup>j</sup> <sup>i</sup>*+<sup>1</sup> are equal to a local constant *λ*∗ in (23). The key part of the von Neumann analysis is to assume a solution of the form [39]:

$$\mathcal{W}\_i^j = \xi^j \exp(\varphi qil\hbar),\tag{26}$$

where *<sup>q</sup>* <sup>=</sup> √−1, *<sup>h</sup>* is the element size, *<sup>ϕ</sup>* is the mode number, and *<sup>ξ</sup><sup>j</sup>* is the amplification factor at time level *j*. As *j* increases (more time steps are computed), adding the latter *W<sup>j</sup> i* expression to the system (23) gives the following characteristic equation:

$$\begin{array}{l} \xi^{j}(A\_{i}\exp(\varrho q(i-1)h) + B\_{i}\exp(\varrho qih) + \mathbb{C}\_{i}\exp(\varrho q(i+1)h)) \\ = \xi^{j-1}(A\_{i}^{\*}\exp(\varrho q(i-1)h) + B\_{i}^{\*}\exp(\varrho qih) + \mathbb{C}\_{i}^{\*}\exp(\varrho q(i+1)h)), \end{array} \tag{27}$$

where:

$$A\_i = k - \beta + \beta k \lambda^\*, A\_i^\* = -\beta,$$

$$B\_i = -2k - \gamma + \gamma k \lambda^\*, B\_i^\* = -\gamma,$$

$$C\_i = k - \beta + \beta k \lambda^\*, C\_i^\* = -\beta,$$

After simple basic calculations, (27) becomes:

$$\zeta^{\ddagger} = \frac{A\_i^\* \exp(-q\phi) + B\_i^\* + C\_i^\* \exp(q\phi)}{A\_i \exp(-q\phi) + B\_i + C\_i \exp(q\phi)},\tag{28}$$

where φ = *ϕh*. After using Euler's formula exp(*q*φ) = cos φ + *q* sinφ, the last equation (28) can be simplified as:

$$\zeta = \frac{\left(A\_i^\* + C\_i^\*\right)\cos\phi + B\_i^\* + q\left(C\_i^\* - A\_i^\*\right)\sin\phi}{\left(A\_i + C\_i\right)\cos\phi + B\_i + q\left(C\_i - A\_i\right)\sin\phi}.\tag{29}$$

Equation (29) can be rewritten in the forms:

$$\xi = \frac{-2\beta\cos\Phi - \gamma}{(2k - 2\beta + 2\beta k\lambda^\*)\cos\Phi + (-2k - \gamma + \gamma k\lambda^\*)}.\tag{30}$$

or

$$\xi = \frac{(2\beta\cos\Phi + \gamma)}{(2\beta\cos\Phi + \gamma) + 2k(1 - \cos\Phi) - \lambda^\*(2\beta k \cos\Phi + \gamma k)}. \tag{31}$$

The quantity (1 − cos φ) is positive or equal to zero but 2*β* cos φ + *γ* is surely positive if we choose *γ* > 0 and *β* > 0 such that *γ* > 2*β*. If we choose *k*, *β*, and *γ* small enough to make *λ*∗(2*βk* cos(φ) + *γk*) → 0, then last (31) is close to:

$$\zeta = \frac{(2\beta\cos\phi + \gamma)}{(2\beta\cos\phi + \gamma) + 2k(1 - \cos\phi)}.\tag{32}$$

For stability, we must have |*ξ*| ≤ 1 (otherwise *ξ<sup>j</sup>* in (26) would expand in an unbounded manner). This condition is satisfied for *γ* > 0 and *β* > 0 such that *γ* > 2*β*. Finally, we can say that our numerical scheme is conditionally stable for *γ* > 0, *β* > 0, and *γ* > 2*β*, if *k*, *β*, and *γ* values chosen are small enough.

#### **4. Numerical Results**

In this part, we employ the suggested method for solving the nonlinear space-fractional Fisher's equation [17,18] with *<sup>α</sup>* <sup>=</sup> 1.5, *<sup>f</sup>*(*x*, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup>*x*2, and *<sup>u</sup>*(*x*, 0) <sup>=</sup> *<sup>x</sup>* and subject to conditions:

```
u(0.0125, t) ≈ 0.0125(1 + t) +0.00609375t
                                         2 −0.082176t
                                                      3 −0.0210541t
                                                                    4 −7.16634×10−6 t
                                                                                        5,
```

```
u(1.0125, t) ≈ 1.0125(1 + t) − 0.518906t
                                          2 − 0.921366 t
                                                         3 + 0.310529t
                                                                       4 + 0.0845434 t
                                                                                       5.
```
The solution of the space-fractional Fisher's equation, like that of most fractional partial differential equations, is unknown, so we will compare it with the known computational approaches such as GDTM and VIM. We used *k* = 0.0005 and *N* = 80 in all our numerical results. Tables 1–3 compare our method, which was improved in Section 2, to other existing methods.

**Table 1.** The comparison between the proposed method and other existing methods for *t* = 0.1 and various *x* values.



**Table 2.** The comparison between the proposed method and other existing methods for *t* = 0.2 and various *x* values.

**Table 3.** The comparison between the proposed method and other existing methods for *t* = 0.4 and various *x* values.


Tables 1–3 show that for all values of *x* and *t*, our approximate solutions are in good agreement with the approximate solutions using other existing methods. This is also shown in Figures 1 and 2.

**Figure 1.** The behaviors of the numerical solution *W*(*x*, *t*) for various time levels.

**Figure 2.** Graphs of the approximate solutions *<sup>W</sup>*(*x*, *<sup>t</sup>*) for *<sup>t</sup>* <sup>≤</sup> 0.5 (*h*<sup>=</sup> 0.0125, *<sup>k</sup>* <sup>=</sup> 0.0005).

#### **5. Conclusions**

The nonlinear space-fractional Fisher's equation is numerically treated in this article using a collocation method with cubic polynomial spline functions. The proposed method is demonstrated to be conditionally stable using the von Neumann stability technique. Compared to other useful techniques, such as VIM, GDTM, and QPSM, the numerical results illustrate that our proposed method retains efficiency, capability, and good versatility. The findings of our study are hoped to pique the attention of academics interested in using the cubic polynomial spline to numerically solve nonlinear fractional partial differential equations of the same type. It is proven that our scheme's local truncation error is to *O*(*hα*), 1 < *α* ≤ 2. Furthermore, the novelty of the present work lies in the fact that the results illustrate that the proposed technique is convenient, accurate, and very efficient in solving such problems.

**Author Contributions:** A.R.H., F.E.A.A., A.A.A. and T.R.; methodology, A.R.H., F.E.A.A., A.A.A. and T.R.; software, A.R.H., F.E.A.A., A.A.A. and T.R.; validation, A.R.H., F.E.A.A., A.A.A. and T.R.; formal analysis, A.R.H., F.E.A.A., A.A.A. and T.R.; investigation, A.R.H., F.E.A.A., A.A.A. and T.R.; resources, A.R.H., F.E.A.A., A.A.A. and T.R.; data curation, A.R.H., F.E.A.A., A.A.A. and T.R.; writing—original draft preparation, A.R.H., F.E.A.A., A.A.A. and T.R.; writing—review and editing, A.R.H., F.E.A.A., A.A.A. and T.R.; visualization, A.R.H., F.E.A.A., A.A.A. and T.R.; supervision, A.R.H.; project administration, T.R.; funding acquisition, All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No new data were created or analyzed in this study. Data sharing is not applicable to this article.

**Acknowledgments:** The researchers would like to thank the Deanship of Scientific Research, Qassim University, for funding the publication of this project.

**Conflicts of Interest:** All the authors declare that they have no conflict of interest.

#### **References**

