On the numerical solution of Fisher's equation by an efficient algorithm based on multiwavelets

1.   Introduction

The balance between linear diffusion and nonlinear reaction or multiplication was studied in the 1930s by Fisher ^[12]. The generalized Fisher's equation with boundary and initial conditions is given as

where $\mu$ is the (constant) reaction factor and $f$ is the nonlinear reaction term. This equation was first proposed to show a model for the propagation of a mutant gene, with $u$ denoting the density of an advantageous. This equation is encountered in population dynamics and chemical kinetics, which includes problems such as the neutron population in a nuclear reaction, the nonlinear evolution of a population in a one-dimensional habitat ^[14,19].

Many methods and computational techniques can be employed to deal with these topics. Mittal et al. ^[20] used an efficient B-spline scheme to solve the Fisher's equation. They proved the stability of the method and reduced the computational cost. The Sinc collocation method is proposed in ^[18] for solving this equation. Mittal et al. ^[21] have obtained numerical solutions of equation (1.1) using wavelet Galerkin method and they have shown that the present method can be computed for a large value of the linear growth rate. Olmos et al. ^[22] developed an efficient pseudospectral solution of Fisher's equation. The viability of applying moving mesh methods to simulate traveling wave solutions of Fisher's equation is investigated in ^[23]. Cattani et al. ^[9] proposed mutiscale analysis of the equation and this article is one of the first articles that has solved the problem numerically. For more related numerical results, we refer the interested readers to ^[1,11,13,16] and references therein.

Mixed methods including finite difference and Galerkin or collocation method have been used to solve the various PDEs. The use of this method can be observed in many studies. Here are some studies which we refer to them. Başan ^[4] applied a mixed algorithm based on Crank–Nicolson mixed by modified cubic B-spline DQM to solve the coupled KdV equation. In ^[7], mixed methods including quintic B-spline and Crank-Nicolson is utilized for solving the nonlinear Schrödinger equation. In ^[5], a numerical solution to the coupled Burgers' equation via contributions of the Crank-Nicolson and differential quadrature method. In ^[28], an algorithm is proposed based on $\theta$-weighted method and wavelet Galerkin method to solve the Benjamin-Bona-Mahony equation. For more related cases, we refer the reader to ^[6,8,25].

Wavelets and specially multiwavelets are found an interesting basis for solving a variety of equations ^[10,27]. Multiwavelets have some properties of wavelets, such as orthogonality, vanishing moments and compact support. Note that they can be both symmetric and orthogonal. In contrast to the wavelets and biorthogonal wavelets they can have high smoothness and high approximate order coupled with short support ^[17]. In addition, some multiwavelets such as Alpert's multiwavelets have the interpolating properties. Contrary to biorthogonal wavelets, multiwavelets can have the high vanishing moments without enlarging their support ^[15]. These bases are suited for high-order adaptive solvers of partial differential equations also the integro-differential equations have sparse representations in these bases. In general, the multiwavelets are a very powerful tool for expressing a variety of operators. At the present work, we apply the Alpert's multiwavelets constructed in ^[2,3].

This paper is organized as follows: In Section 2, we review the definition and properties of the Alpert multi-wavelets required for our subsequent development. In Section 3, we proceed to the main results where construction of a convergent method is given by multiwavelets Galerkin method and the proposed method is examined along with the analysis of the convergence and stability. Section 4, contains some numerical examples to illustrate the efficiency and accuracy of the scheme.

2.   Alpert's multi-wavelet

Assume that $\Omega: = [\alpha, \beta] = \cup_{b\in \mathcal{B}}X_{J, b}$ is the finite discretizations of $\Omega$ where $X_{J, b}: = [x_b, x_{b+1}], b\in\mathcal{B}: = \{\alpha 2^J, \ldots, (\beta-1)2^J+2^J-1\}$ with $J\in \mathbb{Z}^+\cup \{0\}$, are determined by the point $x_b: = b/(2^J)$. On this discretization, applying the dilation $\mathcal{D}_{2^j}$ and the translation $\mathcal{T}_b$ operators to primal scaling functions $\{\phi_{0, \alpha 2^J}^0, \cdots, \phi_{0, (\beta-1)2^J+2^J-1}^{r-1}\}$, one can introduce the subspaces

of scaling functions. Here $\mathcal{R} = \{0, 1, \cdots, r-1\}$ and the primal scaling functions are the Lagrange polynomials of degree less than $r$ that introduced in ^[2].

Every function $p\in L^2(\Omega)$ can be represented in the form

where $\langle., . \rangle$ denotes the $L^2$-inner product and $\mathcal{P}_J^r$ is the orthogonal projection that maps $L^2(\Omega)$ onto the subspace $V_J^r$. To find the coefficients $p_{J, b}^k$ that are determined by $\langle p, \phi_{J, b}^k\rangle = \int_{X_{J, b}}f(x)\phi_{J, b}^k(x)dx$, we shall compute these integrals. We apply the $r$-point Gauss-Legendre quadrature by a suitable choice of the weights $\omega_k$ and nodes ${\tau}_k$ for $k\in \mathcal{R}$ to avoid these integrals ^[2,26], via

Convergence analysis of the projection $\mathcal{P}_J^r(p)$ is investigated for the $r$-times continuously differentiable function $p\in \mathbb{C}^r(\Omega)$.

For the full proof of this approximation and further details, we refer the readers to ^[3]. Thus we can conclude that $\mathcal{P}_J^r(p)$ converges to $p$ with rate of convergence $O(2^{-Jr})$.

Let $\Phi_J^r$ be the vector function ${\Phi_J^r}: = [\Phi_{r, J, \alpha 2^J}, \cdots, \Phi_{r, J, (\beta-1)2^J+2^J-1}]^T$ and consists of vectors $\Phi_{r, J, b}: = [\phi_{J, b}^0, \cdots, \phi_{J, b}^{r-1}]$. The vector function $\Phi_J^r$ includes the scaling functions and called multi-scaling function. Furthermore, by definition of vector $P$ that includes entries $p_{J, b}^k$, we can rewrite Eq (2.2) as follows

where $P$ is an $N$-dimensional vector $(N: = (\beta-\alpha)r2^J)$. The building blocks of these bases construction can be applied to approximate a higher-dimensional function. To this end, one can introduce the two-dimensional subspace $V_J^{r, 2}: = V_J^r\times V_J^r\subset L^2(\Omega)^2$ that is spanned by

Thus by this assumption, to derive an approximation of the function $p\in L^2(\Omega)^2$ by the projection operator $\mathcal{P}_J^r$, we have

where components of the square matrix ${P}$ of order $N$ are obtained by

where $\hat{\tau}_k = (\tau_k+1)/2$. Consider the $2r$-th partial derivatives of $f:\Omega^2\rightarrow \mathbb{R}$ are continuous. Utilizing this assumption, the error of this approximation can be bounded as follows

where $\mathcal{M}_{\max}$ is a constant.

By reviewing the spaces $V_J^r$, it is obvious these bases are nested. Hence there exist complement spaces $W_J^r$ such that

where $\oplus$ denotes orthogonal sums. These subspaces are spanned by the multi-wavelet basis

According to (2.8), the space $V_{J}$ may be inductively decomposed to $V_J^r = V_0^r\oplus (\oplus_{j = 0}^{J-1}W_j^r)$. This called multi-scale decomposition and spanned by the multi-wavelet bases and single-scale bases. This leads us to introduce the multi-scale projection operator $\mathcal{M}_J^r$. Assume that the projection operator $\mathcal{Q}_j^r$ the maps $L^2(\Omega)$ onto $W_j^r$. Thus we obtain

and consequently, any function $p\in L^2(\Omega)$ can be approximated as a linear combination of multi-wavelet bases

where

Note that, we can compute the coefficients $p_{0, 0}^k$ by using (2.2). But multi-wavelet coefficients from zero up to higher-level $J-1$ in many cases must be evaluated numerically. To avoid this problem, we use multi-wavelet transform matrix $T_J$, introduced in ^[24,26]. This matrix connects multi-wavelet bases and multi-scaling functions, via,

where ${{{\Psi}_J^r}}: = [\Phi_{r, 0, b}, \Psi_{r, 0, b}, \Psi_{r, 1, b}, \cdots, \Psi_{r, J-1, b}]^T$ is a vector with the same dimension $\Phi_J^r$ (here $\Psi_{r, j, b}: = [\psi_{j, b}^0, \cdots, \psi_{j, b}^{r-1}]$). This representation helps to rewrite Eq (2.10) as to form

where we have the $N$-dimensional vector ${\tilde{P}_J}$ whose entries are $p_{0, 0}^{k}$ and $\tilde{p}_{j, b}^k$ and is given by employing the multi-wavelet transform matrix $T_J$ as ${\tilde{P}_J} = T_J{P_J}$. Note that according to the properties of $T_J$ we have $T_J^{-1} = T^T_J$.

The multi-wavelet coefficients (details) become small when the underlying function is smooth (locally) with increasing refinement levels. If the multi-wavelet bases have $N_{\psi}^r$ vanishing moment, then details decay at the rate of $2^{-JN_{\psi}^r}$ ^[15]. Because vanishing moment of Alpert's multi-wavelet is equal to $r$, one can obtain ${\tilde{p}_{J, b}^k}\approx O(2^{-Jr})$ consequently. This allows us to truncate the full wavelet transforms while preserving most of the necessary data. Thus we can set to zero all details that satisfy a certain constraint $\varepsilon$ using thresholding operator $\mathcal{C}_{\varepsilon}$

and the elements of $\bar{P}_J$ are determined by

where $D_\varepsilon: = \{(j, b, k):|\tilde{p}_{j, b}^k| > \varepsilon\}$. Now we can bound the approximation error after thresholding via

where $\mathcal{P}_{J, D_{\varepsilon}}^r (p)$ is the projection operator after thresholding with the threshold $\varepsilon$ and $C_{thr} > 0$ is constant independent of $J, \varepsilon$.

3.   The mixed finite difference and Galerkin method (MFDGM)

The building block for the time discretization is the $\theta$-weighted scheme applied to the generalized Fisher's equation. For this purpose we introduce the time discretization $t_{n+1}: = t_n+\delta t$, where the time step $\delta t$ is assumed to be constant. To derive a stable method, we will see in next section that $\delta t$ must satisfy to certain conditions. Here we consider generalized Fisher's equation with initial and boundary conditions

where $\kappa$ and $\mu$ are the constants. We suppose that the initial data $g(x)$ is several times differentiable and the nonlinear term $uu^{\kappa}$ satisfies a Lipschitz condition with Lipschitz constant $\mathcal{L}$.

Let $u^n: = u(x, t_n)$ where $t_n = n\delta t$, $n\in \mathcal{N} = \{0, 1, \ldots, \frac{T}{\delta}\}$. Then the $\theta$-weighted scheme reads

where $0\leq \theta \leq 1$ and $\mathcal{R} < \delta t C$ is a small term for a positive constant $C$. Omitting the small term $\mathcal{R}$, and rearranging (3.2), we obtain from Crank–Nicolson method ($\theta = \frac{1}{2}$),

To derive a multiwavelets Galerkin method for the generalized Fisher's equation (3.1), assume that the approximate solution for each step can be written as

where $U_n$ is an $N$-dimensional vector whose elements must be found. The derivative operator $\frac{d^2}{dx^2}$, can be approximated by

where $D_{\psi}$ is the operational matrix of derivative for multiwavelets introduced in ^[8].

Inserting (3.4) and (3.5) into (3.3) and using multi-scale projection operator for $(u{u^\kappa})^n$, yields

where $(u{u^\kappa})^n\approx E_n^T\Psi_J^r$. To obtain the approximate solution of (3.6) using multiwavelets Galerkin method, we multiply (3.6) by $\Psi_J^T(x)$ and integrate over its support $\Omega$. Therefore using the orthonormality of this system of multiwavelets, we have

where $M_2$ is the right hand side of (3.6), and

The system (3.7) consists of $N$ equations. Since two of them are linearly dependent, we replace these dependent equations with two others that have derived from boundary conditions (3.1). A new system of equations is obtained by replacing the first and last columns of $M_1$ with $\Psi_J^r(\alpha)$ and $\Psi_J^r(\beta)$ and the first and last elements of $M_2$ with $h_1(t^{n+1})$ and $h_2(t^{n+1})$, i.e.,

To start the steps, we use the initial condition. The initial condition (3.1) can be approximated as

Putting $U_0 = \mathcal{G}$, a linear system of algebraic equations arise such that by solving this system, the unknown coefficients $U_{n+1}$ at every time steps can be found.

3.1.   Convergence analysis

Let $H^k(\Omega) = W^{k, 2}(\Omega)$ be the usual Sobolev space of order $k$ on $\Omega$. With this definition, the Sobolev spaces $H^2$ admit a natural norm

where $\langle., .\rangle$ is the $L^2(\Omega)$-inner product.

Theorem 3.1. Assume that the nonlinear term $p(u): = uu^{\kappa}$ satisfies Lipschitz condition with respect to $u$ as,

where lipschitz constant $\mathcal{L} < \infty$ is supposed to be large enough. Then, the time discrete numerical scheme defined by (3.3) is stable in $H^2$-norm when $\delta t < \frac{1}{\mu \mathcal{L}}$.

Proof. Subtraction equation (3.3) from

one can find the roundoff error $e^n = u^n-\hat{u}^n$ for $n = 0, 1, \ldots, \frac{T}{\delta t}$, as

where $u^n$ and $\hat{u}^n$ are the exact and approximate solutions of (3.3), respectively.

Applying the Lipschitz condition (3.10) and then simplifying, it follows that

Multiplying (3.12) by $e^{n+1}$ and integrating on $\Omega$, yield

Performing integration by parts, we obtain

Now, it follows from the Schwarz inequality $(|\langle v, w \rangle|\leq \|v\|\|w\|)$ that

where we used the inequality $vw\leq \frac{1}{2}(v^2+w^2)$. Rearranging (3.14), we have

Since $\delta t < \frac{1}{\mu \mathcal{L}}$ and $\mathcal{L}$ is assumed to be sufficiently large, one has

and then it follows from $1-\mu \delta t \mathcal{L} > 0$ that

Hence, one can find for $n = 0, 1, \ldots, \frac{T}{\delta t}-1$

The proof completes by taking the limit as $n\rightarrow \infty$,

Theorem 3.2. Under the assumptions of Theorem 3.1, the time discrete solution $\hat{u}^n$ is $H^2$-convergent to $u^n$.

Proof. Assume that $e^n = u^n-\hat{u}^n$ is the perturbation error. Since $\hat{u}^n$ is the approximate solution of (3.2) at the time step $n$ which satisfies initial and boundary conditions (3.1), it follows that $e^n$ satisfies (3.2)

Applying Lipschitz condition (3.10), one can write after simplification

Multiplying (3.20) by $e^{n+1}$, and integrating over $\Omega$

Using integration by parts, one has

It follows from the Schwarz inequality that

Further, by the Young's inequality ($ab\leq \frac{1}{2\varepsilon}a^2+\frac{\varepsilon}{2}b^2$) with ($\varepsilon = \delta t$) we obtain

By simplification of the above relation, we obtain

Since $\delta t < \frac{1}{\mu \mathcal{L}}$, it follows that

By repeating this relation for $n = 0, 1, \cdots, \frac{T}{\delta t}-1$, one can write

Because $e^0 = 0$, one can show that

Then, since

it follows from $\mathcal{R}\leq C\delta t$ that

It is obvious that $\delta t\rightarrow 0$ as $n\rightarrow \infty$ and then we obtain

This completes the proof.

4.   Numerical experiments

In this section, some numerical examples are presented to illustrate the validity and the merits of the new technique. The accuracy of the method has been measured by $L_2$-error i.e.,

where $i = n\delta t/2^m$, $n = 0, \ldots, 10T(2^{m-1})-1$. In all examples, we assume that the primal time step size is $\delta t: = 0.1$.

Example 4.1. Consider the Fisher's equation:

The exact solution is given in ^[29]

and the initial and boundary condition can be extracted by the exact solution.

The effects of the refinement level $J$, multiplicity parameter $r$ and time step size $\delta t$ on $L_2$-error are given in Table 1. Figure 1 is plotted to show the effect of time step size on the accuracy. As the time step size increases, it can be seen that the error decreases, and the approximate solution converges to the exact solution. Figure 2 illustrate the approximate solution and $L_{\infty}$-error taking $r = 3$, $J = 2$ and $m = 8$. Table 2 displays $L^2$-error using the presented method taking $r = 3$, $J = 2$, $\delta t = 0.1/2^m$, $m = 1, \ldots, 8$. The results have been compared with implicit $(\theta = 1)$ and explicit method $(\theta = 0)$.

Example 4.2. Consider the Fisher's equation:

The exact solution is given by ^[14]

and the boundary and initial conditions can be obtained by it.

Table 3 shows the effects of the refinement level $J$, multiplicity parameter $r$ and time step size $\delta t$ on $L_2$-error. Figure 3 is also provided for further observations. Figure 4 shows that the approximate solution converges to the exact solution as the time step size increases. The approximate solution and $L_{\infty}$-error is presented graphically for $r = 3$, $J = 2$ and $m = 8$ and the results are shown in Figure 5. The results prove the efficiency and accuracy of the proposed method. In Table 4, we compare the $L^2$-errors taking $r = 3$, $J = 2$ and $\delta t = 0.1/2^{m-1}$ at time $t = 1$ between the Crank-Nicolson method and implicit $(\theta = 1)$.

Example 4.3. Consider the Fisher's equation:

The exact solution is given by ^[14]

Figure 6 shows the approximate solution and $L_2$-error taking $r = 3$, $J = 2$ and $m = 8$. One can see the effect of the refinement level $J$, the multiplicity parameter $r$ and time step size, on $L_2$ error in Figure 7. We observe that with increasing the refinement level $J$ and the multiplicity parameter $r$ the $L_2$ error decreases. Table 5 displays $L^2$-error using the presented method taking $r = 3$, $J = 2$, $\delta t = 0.1/2^m$, $m = 1, \ldots, 8$. The results have been compared with implicit $(\theta = 1)$ and explicit method $(\theta = 0)$.

5.   Conclusion

Multiwavelets Galerkin method is applied to solve the Fisher's equation. After discretization of time using the Crank-Nicolson method, a system of ordinary differential equations arises at any time step. Then Multiwavelets Galerkin method is used to solve this system of equations. The result of applying the method is a nonlinear system of algebraic equations at any time step. By solving this system, one can find the approximate solution at any time. The convergence and stability analysis are investigated, and numerical simulations indicate that the proposed method gives a satisfactory approximation to the exact solution.

Acknowledgments

This project was supported by Researchers Supporting Project number (RSP-2020/210), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

The author declares no conflicts of interest in this paper.