An Explicit Wavelet Method for Solution of Nonlinear Fractional Wave Equations

Abstract

An explicit method for solving time fractional wave equations with various nonlinearity is proposed using techniques of Laplace transform and wavelet approximation of functions and their integrals. To construct this method, a generalized Coiflet with N vanishing moments is adopted as the basis function, where N can be any positive even number. As has been shown, convergence order of these approximations can be N. The original fractional wave equation is transformed into a time Volterra-type integro-differential equation associated with a smooth time kernel and spatial derivatives of unknown function by using the technique of Laplace transform. Then, an explicit solution procedure based on the collocation method and the proposed algorithm on integral approximation is established to solve the transformed nonlinear integro-differential equation. Eventually the nonlinear fractional wave equation can be readily and accurately solved. As examples, this method is applied to solve several fractional wave equations with various nonlinearities. Results show that the proposed method can successfully avoid difficulties in the treatment of singularity associated with fractional derivatives. Compared with other existing methods, this method not only has the advantage of high-order accuracy, but it also does not even need to solve the nonlinear spatial system after time discretization to obtain the numerical solution, which significantly reduces the storage and computation cost.

1. Introduction

In some complex mechanical and physical processes, many empirical formulae are often expressed in the form of power-law functions, and corresponding relations are usually not the laws in forms of standard derivatives. These processes sometimes exhibit obvious properties of memory, heredity, and path dependency. Under such a situation, when the classical integer-form derivatives are used to quantitatively describe the above problems, there is often a need to construct complex nonlinear equations and introduce some artificial empirical parameters and assumptions that are inconsistent with the reality. These nonlinear models are usually very troublesome in theoretical analysis and numerical solution. Fractional calculus, including the fractional derivative and integral, becomes one of the important tools for mathematical modeling of complex mechanical and physical processes because they can succinctly and accurately describe the historical memory and spatially nonlocal correlation in these processes.

In the last few decades, applied scientists and engineers have realized that integrating fractional derivatives on the basis of conventional differential equations can provide an elegant and natural framework for analyzing various practical problems that need to be modeled with time memory and spatial nonlocal characteristics. In these applications, the corresponding fractional derivatives, unlike the conventional integer-order ones, have several different definitions, which include the Caputo fractional derivative (CFD), the Riemann–Liouville fractional derivative (RFD), the Gruenwald–Letnikov fractional derivative, the Riesz fractional derivative, and several others [1,2]. Among these definitions, the CFD and the RFD are most commonly used. Noting that the CFD is only suitable for differentiable functions and the RFD of a constant is not zero, Abu-Shady et al. [3,4] proposed a generalized fractional derivative (GFD) definition that could well eliminate these drawbacks.

An exact and unique geometric and physical interpretation of fractional calculus is still an unsolved open question, although there is a lot of research in this field [5,6,7]. However, a common recognition in fractional calculus is that the fractional derivative could give a power-law approximation of the local behavior of functions, and could be used in science and engineering to investigate the behavior of objects and systems that are characterized by the power-law non-locality and power-law long-term memory.

In recent years, fractional calculus has become an important tool to describe complex behaviors in a variety of scientific and engineering fields, although differences in their definitions exist. Until now, examples of applications have successfully included the fields of physics science [8,9,10,11,12,13,14,15], biological science [16,17,18], material science [19], control techniques [20], environmental science [21], economics [22], photovoltaics [23,24], electrochemistry [25], and multidisciplinary engineering [26], etc. Among them, a class of highly representative and effective applications of fractional derivative models is the fractional wave equation (FWE) established for the abnormal diffusion and wave problems that exist widely in the fields of viscoelasticity, quantum physics, and material physics.

In this study, we consider the following nonlinear time-FWE:

$$\{\begin{array}{l}\rho D_t^{\alpha }u+\lambda D_t^{\beta }u=\Delta u+f\left(u\right)+h\left(x,y,t\right),1<\alpha \le 2,\beta <\alpha \\ u\left(x,y,0\right)=g_0\left(x,y\right),u_t\left(x,y,0\right)=g_1\left(x,y\right),\hspace{1em}\hspace{1em}\hspace{1em}\left(x,y\right)\in \Omega ,\\ u\left(x,y,t\right)=q\left(x,y,t\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(x,y\right)\in \partial \Omega ,t>0,\end{array}$$

(1)

in which Ω = [0, a] × [0, a] is a bounded domain with boundary ∂Ω; α and β are the parameters describing the derivative of the fractional order with respect to time; ρ, λ, and μ are positive constants; and h(x, y, t) is a source term. As the CFD definition has been widely adopted in most engineering applications [8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26], we also consider the CFD for the terms $D_t^{\alpha }u\left(x,y,t\right)\mathrm{and}{D}_t^{\beta }u\left(x,y,t\right)$ in the time-FWE of Equation (1), which is

$$D_t^{\alpha }u\left(x,y,t\right)=\{\begin{array}{l}\frac{{\partial }^{n}u\left(x,y,t\right)}{\partial t^n}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\alpha =n,n\in N,\\ \frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_0^t\frac{{\partial }^{n}u\left(x,y,s\right)}{\partial s^n}\frac{1}{{\left(t-s\right)}^{\alpha +1-n}}}ds,n-1<\alpha <n.\end{array}$$

(2)

We note that the initial values of Equation (1) associated with the CFD can be similar to those of the conventional integer-order differential equations, which make the equation easier to deal with, and the physical significance and extensive practical applications easier to understand. Equation (1) can be reduced to several important cases. For example, if we put λ = 0, f(u) = 0 or u then we obtain the time FWE without or with damping [8,9,10]. If we set f(u) = u and β = α/2, then we can have the fractional telegraph equation [11,12,13,14]. When λ = 0 and f(u) = sin(u), the time fractional nonlinear sine-Gordon equation is obtained [27]. When λ = 0 and f(u) =−r₁u−r₂u²−r₃u³, then we get the time fractional nonlinear Klein–Gordon equation [15,27,28].

For most FWEs, it is usually difficult to obtain their analytical solutions. Even for numerical solutions, the FWEs seem to be much more difficult to accurately, reliably, and effectively deal with than those with standard integer orders due to the coupling between nonlinear and nonlocal effects in Equation (1), which poses significant obstacles to the development of effective solution methods. Despite these difficulties, it is still very urgent that we develop high-precision and efficient algorithms to accurately solve these equations. Many researchers have tried to propose various numerical methods to solve such types of nonlinear fractional wave equations. As representative examples, Huang et al. [8] constructed two types of finite difference schemes to solve the time FWEs and proved that both of the two schemes are convergent with first-order accuracy in the temporal direction and second-order accuracy in the spatial direction. Sun and Wu [29] proposed a fully discrete difference method for the solution of linear fractional diffusion-wave systems and analyzed the stability and convergence of the method. However, this approach adds complexity by introducing two new variables to transform the original equation into a lower-order system of equations. Du et al. [30] proposed another difference method with improved convergence order for the solution of linear FWEs. As the implementation of finite difference schemes is relatively simple, research on the use of such methods in solving FWEs is very rich [31]. However, there inevitably exists a contradiction between the nonlocal characteristics of the fractional derivatives and the inherent finite precision and local characteristics of the schemes of finite difference [32]. In order to improve the computational efficiency, Liang et al. [11] suggested a fast high order difference scheme based on an efficient sum-of-exponentials approximation to solve the time fractional differential equations. However, this is currently only applicable to some linear problems. There have also been many attempts at methods other than the difference discretization methodology. For examples, Bhrawy et al. [9] have proposed a spectral algorithm based on the Jacobi operational matrix to solve the second- and fourth-order FWEs, and Ghafoor et al. [33] adopted the wavelet collocation method to discretize spatial operators and only used the finite difference scheme to approximate the time fractional derivatives for the (1 + 1)- and (1 + 2)-dimensional time fractional diffusion-wave equations.

Most of the studies mentioned above have only applied to linear problems. For the solution of nonlinear fractional wave problems, Vong and Wang [28] proposed a high order finite difference scheme for a two-dimensional fractional Klein–Gordon equation with Neumann boundary conditions, which realized 4th order convergence for each spatial dimension. Doha et al. [34] introduced a new approach implementing a shifted Jacobi operational matrix in combination with the shifted Jacobi collocation technique for the numerical solution of nonlinear multi-term FWEs. The resulting nonlinear algebraic equations were solved using Newton’s iterative method. Dehghan et al. [27] proposed a numerical method for the solution of time fractional nonlinear sine-Gordon equations. Their work showed that the meshless method based on the radial basis functions and collocation approach are also suitable for the treatment of the nonlinear FWEs. However, such a method needs to solve a linear boundary-value problem at each time step during the iterative solution process, which can be computational expensive for high-dimensional problems. Lyu and Vong [15] suggested a difference scheme for nonlinear time fractional Klein–Gordon-type equations with second-order accuracy in time. In their approach, the nonlinearity was dealt with by a linearized scheme. However, such a treatment may not be valid for problems with derivative nonlinearity.

Despite the above progresses, there are still some critical issues in the numerical solution of the nonlinear FWEs. Nonlinearity and nonlocality natures pose major challenges to the accurate and efficient solution of these problems. More specifically, in the development of solution methods, it is still very difficult to avoid iterative processes and complicated inversion of spatial derivative operators at each time step to save computational cost, and to make the proposed methods valid for problems with stronger nonlinearities, not just those of the Klein– and sine-Gordon types.

Wavelet theory has mainly been developed in the last few decades, and has aroused great interest in many research fields of mathematics, physics, computer science, and engineering [35]. In recent years, researchers have made very effective progresses in using wavelet theory to solve FWEs. For example, Liu et al. [36] obtained exact solutions of several types of one-dimensional linear fractional partial differential equations by using an exact formula of Laplace inversion based on the generalized Coiflet (a class of wavelet) developed by Wang et al. [37]. Zhou et al. [38] proposed a generalized Coiflet–Galerkin method to solve nonlinear fractional vibration, diffusion, and wave differential equations. The efficiency of this method was justified by numerical solution of one-dimensional nonlinear examples. These methods [36,37,38] successfully adopted the quasi-interpolation property of wavelet approximation, which can expand relevant functions and nonlinear terms to explicit forms with coefficients being their single-value sampling. This property enables the effective solution of problems involving convolution integrals and nonlinear terms.

By noting the potential advantages of Coiflet-type basis function in solving nonlinear FWEs, we will adopt the Coiflet-based solution strategy to solve Equation (1). First, Equation (1) is transformed into a time Volterra-type integro-differential equation associated with a smooth time kernel to effectively avoid the treatment of any singular integral kernels. Then, by denoting various spatial derivatives of the unknown function as new functions, the integral relations between these new functions are applied in terms of wavelet approximation of multiple integrals [39,40,41], so that the original equation with spatial derivatives of various orders can be converted to a system of time-integral equations with discrete spatial nodal values of the highest-order spatial derivative. Since the error order of the adopted scheme of integral approximation based on wavelet theory is independent of the order of the integral, the error order caused by the whole discretization process becomes consistent with the error order of the direct approximation of a function, and, interestingly, irrelevant to the equation order [39,40,41]. By transforming the original equation into the integral-differential equation with a smooth integral kernel, we use the wavelet expansion in Equation (8) to approximate the convolution integral on time. On one hand, due to the excellent performance of the wavelet-based approach, the integral approximation still maintains a high convergence rate. On the other hand, we expect that the smooth kernel function has the property of being zero at time zero, which can eliminate the nodal values at the present time during the wavelet approximation of the convolution integral according to a previous study [32]. Finally, an explicit solution procedure based on the collocation method and the proposed algorithm on integral approximation can be established to solve the transformed nonlinear integro-differential equation.

The outline of the paper is as follows. In Section 2, we introduce the scheme of wavelet approximation for interval-bounded functions and their multiple integrals. In Section 3, we describe how the wavelet solution method of Equation (1) is established, and some numerical examples are examined to demonstrate the computational efficiency and accuracy of the proposed method in Section 4. Finally, Section 5 presents the concluding remarks.

2. Wavelet Approximation of Multiple Integrals in a Bounded Domain

Following our previous research [39,40,41,42,43], we consider the Coiflet approximation of f(x)$\in$L²([0, a]) as follows:

$$f(x)\approx {\mathrm{P}}^{m}f(x)={\displaystyle \sum _{k=0}^{a_m}f\left(k/2^m\right)}{\Phi }_{m,k}^{}(x),x\in [0,a]$$

(3)

in which m is the resolution level and a_m = [2^ma] is the integer part of 2^ma and the modified Coiflet basis

$${\Phi }_{m,k}^{}(x)=\{\begin{array}{l}\phi (2^{m}x+M_1-k)+{\displaystyle \sum _{j=-{\alpha }_2}^{-1}{T}_{L,k}(\frac{j}{2^m})\phi (2^{m}x+M_1-j)},0\le k\le {\alpha }_1\\ \phi (2^{m}x+M_1-k),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\alpha }_1+1\le k\le a_m-{\alpha }_2-1\\ \phi (2^{m}x+M_1-k)+{\displaystyle \sum _{j=a_m+1}^{a_m+{\alpha }_1}{T}_{R,a_m-k}(\frac{j}{2^m})\phi (2^{m}x+M_1-j)},a_m-{\alpha }_2\le k\le a_m\end{array}$$

(4)

where $M_1={\displaystyle {\int }_{-\infty }^{\infty }x\phi (x)dx}$ and [0, 3N − 1] are the first-order moment and compact support of the generalized Coiflet-type orthogonal scaling function φ(x), respectively. The scaling function φ(x) with N = 6 and M₁ = 7 is adopted in this study. This scaling function and its corresponding wavelet function are shown in Figure 1.

In addition, in Equation (4), the parameters ${\alpha }_1=M_1-1$ and ${\alpha }_2=3N-2-M_1$, the functions:

$$T_{L,j}(x)={\displaystyle \sum _{i=0}^{N-1}{2}^{im}\frac{{\zeta }_{0,i,j}}{i!}{\left(x\right)}^i}\mathrm{and}{T}_{R,j}(x)={\displaystyle \sum _{i=0}^{N-1}{2}^{im}\frac{{\zeta }_{1,i,j}}{i!}{(x-2^{-m}{a}_m)}^i}$$

(5)

in which the coefficients $\left\{{\zeta }_{0,i,k}\right\}$ and $\left\{{\zeta }_{1,i,k}\right\}$ can be obtained through the relations $P_0=\left\{{\zeta }_{0,i,k}\right\}={(I-B_0)}^{-1}{A}_0$, $P_1=\left\{{\zeta }_{1,i,k}\right\}={(I-B_1)}^{-1}{A}_1$ together with the matrices

$$\begin{array}{l}\hspace{1em}{A}_0=\{{\phi }^{(i)}(M_1-k)\},\hspace{1em}\hspace{1em}\hspace{1em}{A}_1=\{{\phi }^{(i)}(k+M_1)\},\\ B_0=\left\{{\displaystyle \sum _{l=-{\alpha }_1}^{-1}\frac{1}{i!}{l}^i{\phi }^{(j)}(-l+M_1)}\right\},B_1=\left\{{\displaystyle \sum _{l=1}^{M_1-1}\frac{1}{j!}{l}^j{\phi }^{(i)}(-l+M_1)}\right\}.\end{array}$$

(6)

The function values of the scaling function $\phi \left(x\right)$ and its derivatives ${\phi }^{\left(i\right)}\left(x\right)\left(i=1,2,\dots ,N-1\right)$ can be exactly obtained [42].

We define the nth order integral of the function f(x) as

$$f^{\left(-n\right)}\left(x\right)={\displaystyle {\int }_0^x{\displaystyle {\int }_0^{{\xi }_n}\dots }}{\displaystyle {\int }_0^{{\xi }_2}f\left({\xi }_1\right)d{\xi }_1}d{\xi }_2\dots d{\xi }_n.$$

(7)

Performing n-order integration to Equation (3) yields

$$f^{\left(-n\right)}\left(x\right)\approx {\displaystyle \sum _{k=0}^{a_m}f\left(k/2^m\right)}{\Phi }_{m,k}^{\left(-n\right)}\left(x\right),x\in [0,a]$$

(8)

where

$${\Phi }_{m,k}^{\left(-n\right)}(x)=\{\begin{array}{l}{\phi }_{m,k-M_1}^{\left(-n\right)}(x)+{\displaystyle \sum _{j=-{\alpha }_2}^{-1}{T}_{L,k}(\frac{j}{2^m}){\phi }_{m,j-M_1}^{\left(-n\right)}(x)},\hspace{1em}\hspace{1em}0\le k\le {\alpha }_1\\ {\phi }_{m,k-M_1}^{\left(-n\right)}(x),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\alpha }_1+1\le k\le a_m-{\alpha }_2-1\\ {\phi }_{m,k-M_1}^{\left(-n\right)}(x)+{\displaystyle \sum _{j=a_m+1}^{a_m+{\alpha }_1}{T}_{R,a_m-k}(\frac{j}{2^m}){\phi }_{m,j-M_1}^{\left(-n\right)}(x)}.a_m-{\alpha }_2\le k\le a_m\end{array}$$

(9)

The nth order integral of the scaling function basis can be defined as

$${\phi }_{m,k}^{\left(-n\right)}\left(x\right)={\displaystyle {\int }_0^x{\displaystyle {\int }_0^{{\xi }_n}\dots }}{\displaystyle {\int }_0^{{\xi }_2}\phi \left(2^{m}x-k\right)d{\xi }_1}d{\xi }_2\dots d{\xi }_n=\frac{1}{2^{mn}}\left({\phi }^{\left(-n\right)}\left(2^{m}x-k\right)-{\displaystyle \sum _{l=1}^n\frac{{\left(2^{m}x\right)}^{n-l}}{\left(n-l\right)!}{\phi }^{\left(-l\right)}\left(-k\right)}\right)$$

(10)

in which values of ${\phi }^{\left(-n\right)}\left(x\right)$ can be exactly calculated by the method in [41].

For a function f(x) to be smooth enough, the accuracy of the approximation (8) has been estimated as [41]

$${\Vert f^{\left(-n\right)}(x)-f_{{\mathrm{P}}_m}^{\left(-n\right)}(x)\Vert }_{{\mathrm{L}}^2\left[0,a\right]}\le P_{0,L}{2}^{-mN}$$

(11)

in which $P_{0,\mathrm{L}}$ is a constant independent of the resolution level m. From Equation (11), we can see that the convergence order is always N, independent of the integration order.

For the approximation of integrals of a two-dimensional function f(x,y)$\in$L²([0, a]²), Wang et al. [39,40] have suggested that

$$f^{\left(-n_1,-n_2\right)}(x,y)={\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}f(\frac{k}{2^m},\frac{l}{2^m}){\Phi }_{m,k}^{\left(-n_1\right)}(x){\Phi }_{m,l}^{\left(-n_2\right)}(y)}}+O(2^{-mN})$$

(12)

3. Solution Procedure

In this section, we establish a wavelet-integral-collocation method (WICM) to solve the nonlinear FWEs.

Noting that the integral is unlike its derivative, its wavelet-based approximation can ensure a constant high convergence order. Due to this interesting fact, we will first consider using Laplace transform to convert Equation (1) into a time-integral equation. To do so, we denote U(x, y, s) as the time Laplace transform of u(x, y, t), i.e., $\mathcal{L}$[u(x, y, t)] = U(x, y, s), where $\mathcal{L}$[$·$] is the operator representing the time Laplace transform. For the fractional derivative of u(x, y, t), we have

$$\mathcal{L}\left[D_t^{\alpha }u\left(x,y,t\right)\right]=s^{a}U\left(x,y,s\right)-{\displaystyle \sum _{i=0}^{\left[\alpha \right]}{s}^{\alpha -1-i}}{g}_i\left(x,y\right)$$

(13)

where [α] represents the maximal integer less than α.

Applying the time Laplace transform to Equation (1), we have [38]

$$U\left(x,y,s\right)-R\left(s\right)\mathcal{L}\left[p(x,y,t)\right]=G\left(x,y,s\right)$$

(14)

in which

$$\begin{array}{l}p(x,y,t)=\Delta u+f\left(u\right)+h\left(x,y,t\right),\\ R\left(s\right)=1/\left(\rho s^{\alpha }+\lambda s^{\beta }\right),\\ G\left(x,y,s\right)=R\left(s\right)\left(\rho {\displaystyle \sum _{i=0}^{\left[\alpha \right]}{s}^{\alpha -1-i}}{g}_i\left(x,y\right)+\lambda {\displaystyle \sum _{i=0}^{\left[\beta \right]}{s}^{\beta -1-i}}{g}_i\left(x,y\right)\right).\end{array}$$

(15)

Applying inverse time Laplace transform to Equation (14), we have

$$\begin{array}{l}u\left(x,y,t\right)-{\displaystyle {\int }_0^{t}r\left(t-\tau \right)}p\left(x,y,\tau \right)d\tau =g\left(x,y,t\right),\\ r\left(t\right)={\mathcal{L}}^{-1}\left[R\left(s\right)\right]=\frac{1}{\rho }{t}^{\alpha -1}{E}_{\alpha -\beta ,\alpha }\left(-\frac{\lambda }{\rho }{t}^{\alpha -\beta }\right),\\ g\left(x,y,t\right)={\mathcal{L}}^{-1}\left[G\left(x,y,s\right)\right]={\displaystyle \sum _{i=0}^{\left[\alpha \right]}{t}^i{E}_{\alpha -\beta ,1+i}\left(-\frac{\lambda }{\rho }{t}^{\alpha -\beta }\right)}{g}_i\left(x,y\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\frac{\lambda }{\rho }{\displaystyle \sum _{i=0}^{\left[\beta \right]}{t}^i{E}_{\alpha -\beta ,1+i}\left(-\frac{\lambda }{\rho }{t}^{\alpha -\beta }\right)}{g}_i\left(x,y\right).\end{array}$$

(16)

in which we have assumed α > β, E_a,b(z) is the generalized Mittag–Leffler-type function [44,45] defined as $E_{a,b}\left(z\right)={\displaystyle \sum _{k=0}^{\infty }{z}^k/\Gamma \left(ak+b\right)}$. The function g(x, y, t) can be obtained by the initial conditions in Equation (1), and Γ($·$) is the Gamma function. We note that interesting behaviors of Mittag–Leffler-type functions associated with the fractional calculus have been studied in-depth in [45]. It can be easy to verify that the integral kernel, r(t), has a very useful property of r(0) = 0.

Further we can have the following integro-differential equation

$$\{\begin{array}{l}u\left(x,y,t\right)-{\displaystyle {\int }_0^{t}r\left(t-\tau \right)}\left(\Delta u+f\left(u\right)+h\left(x,y,\tau \right)\right)d\tau =g\left(x,y,t\right),\\ u\left(x,y,t\right)=q\left(x,y,t\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(x,y\right)\in \partial \Omega ,t>0.\end{array}$$

(17a)

When using numerical methods to solve this equation, it usually needs to numerically approximate spatial derivatives in the equation, which will inevitably reduce the accuracy of the solution, especially for the cases with high-order derivatives and high spatial dimensions. In order to avoid such a situation, we consider the following treatment [39,40].

We first define $uij={\partial }^{i+j}u/\partial x^i\partial y^j\left(i,j=0,1,2\right)$ as new functions. Then, Equation (17a) can be rewritten as

$$\{\begin{array}{l}u\left(x,y,t\right)-{\displaystyle {\int }_0^{t}r\left(t-\tau \right)}\left(u20+u02+f\left(u\right)+h\left(x,y,\tau \right)\right)d\tau =g\left(x,y,t\right),\\ u\left(x,y,t\right)=q\left(x,y,t\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(x,y\right)\in \partial \Omega ,t>0.\end{array}$$

(17b)

For the two functions u20 and u02 in Equation (17b), applying Equation (12), we obtain

$$u20\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u20(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}(x){\Phi }_{m,l}\left(y\right)}}$$

(18a)

$$u02\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u02(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}(x){\Phi }_{m,l}\left(y\right)}}$$

(18b)

Performing integration to Equation (18) one and two times with respect to x, we have

$$u10\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u20(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}^{\left(-1\right)}(x){\Phi }_{m,l}\left(y\right)+u10\left(0,y,t\right)}}$$

(19a)

$$u01\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u02(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}(x){\Phi }_{m,l}^{(-1)}\left(y\right)}}+u01\left(x,0,t\right)$$

(19b)

and

$$u\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u20(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}^{\left(-2\right)}(x){\Phi }_{m,l}\left(y\right)+xu10\left(0,y,t\right)+q\left(0,y,t\right),}}$$

(20a)

$$u\approx {\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u02(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}(x){\Phi }_{m,l}^{(-2)}\left(y\right)}}+yu01\left(x,0,t\right)+q\left(x,0,t\right).$$

(20b)

Substituting x = a and y = a into Equation (20) we can derive

$$u10\left(0,y,t\right)\approx 1/a\left(q\left(a,y,t\right)-q\left(0,y,t\right)-{\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u20(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,k}^{\left(-2\right)}(a){\Phi }_{m,l}\left(y\right)}}\right)$$

(21a)

$$u01\left(x,0,t\right)\approx 1/a\left(q\left(x,a,t\right)-q\left(x,0,t\right)-{\displaystyle \sum _{k=0}^{a_m}{\displaystyle \sum _{l=0}^{a_m}u02(\frac{k}{2^m},\frac{l}{2^m},t){\Phi }_{m,l}\left(x\right){\Phi }_{m,l}^{\left(-2\right)}(a)}}\right).$$

(21b)

Substituting Equations (18)–(21) into Equation (17b) and then considering the collocation method by taking $x=k^{\prime }/2^m=k^{\prime }\Delta h$ and $y=l^{\prime }/2^m=l^{\prime }\Delta h$ for $k\prime ,l\prime =0,1,\dots ,a_m$, where ∆h represent grid size in both of the x and y directions, the integral Equation (17b) can finally be discretized as

$$\mathrm{A}u20-{\displaystyle {\int }_0^{t}r\left(t-\tau \right)}\left(u02+u20+f\left(u20\right)+h(\tau )\right)d\tau \approx g\left(t\right)-q_1(t),$$

(22)

in which $u02=\mathbf{C}u20{+\mathrm{B}}^{-1}\left(q_1(t)-q_2(t)\right)$, f(u20) indicates that the nonlinear term acts on each coordinate of the vector u20, and the matrices C = B⁻¹A with

$$\begin{array}{l}A=\left\{a_{op}={\Phi }_{m,k}^{\left(-2\right)}\left(\frac{k^{\prime }}{2^m}\right){\Phi }_{m,l}\left(\frac{l^{\prime }}{2^m}\right)-\frac{k^{\prime }}{2^{m}a}{\Phi }_{m,k}^{\left(-2\right)}\left(a\right){\Phi }_{m,l}\left(\frac{l^{\prime }}{2^m}\right)\right\},\\ B=\left\{b_{op}={\Phi }_{m,k}\left(\frac{k^{\prime }}{2^m}\right){\Phi }_{m,l}^{\left(-2\right)}\left(\frac{l^{\prime }}{2^m}\right)-\frac{l^{\prime }}{2^{m}a}{\Phi }_{m,k}\left(\frac{k^{\prime }}{2^m}\right){\Phi }_{m,l}^{\left(-2\right)}\left(a\right)\right\}.\end{array}$$

(23)

The vectors $u20={\left\{u{20}_p^{}=u20\left(k/2^m,l/2^m\right)\right\}}^{\mathrm{T}}$, $g={\left\{g_p=g\left(k/2^m,l/2^m\right)\right\}}^{\mathrm{T}}$, $h={\left\{h_p=h\left(k/2^m,l/2^m\right)\right\}}^{\mathrm{T}}$, and

$$\begin{array}{l}{q}_1\left(t\right)={\{q_{1,o}=\frac{k^{\prime }}{2^{m}a}\left(q\left(a,\frac{l^{\prime }}{2^m},t\right)-q\left(0,\frac{l^{\prime }}{2^m},t\right)\right)+q\left(0,\frac{l^{\prime }}{2^m},t\right)\}}^{\mathrm{T}},\\ q_2\left(t\right)={\{q_{2,o}=\frac{l^{\prime }}{2^{m}a}\left(q\left(\frac{k^{\prime }}{2^m},a,t\right)-q\left(\frac{k^{\prime }}{2^m},0,t\right)\right)+q\left(\frac{k^{\prime }}{2^m},0,t\right)\}}^{\mathrm{T}},\end{array}$$

(24)

and the subscripts $o=\left(a_m+1\right)k^{\prime }+l^{\prime }$, $p=\left(a_m+1\right)k+l$, and $k,l,k\prime ,l\prime =0,1,\dots ,a_m$.

Letting K(t) = Cu20 + u20 + f(u20) + B⁻¹(q₁(t) − q₂(t)) + h(t), we have

$$\mathbf{A}u20-{\displaystyle {\int }_0^{t}r\left(t-\tau \right)K\left(\tau \right)d\tau }\approx g\left(t\right)-q_1(t),$$

(25)

Treating $r\left(t-\tau \right)K\left(\tau \right)$ as a function of $\tau$ in the interval [0, t], applying Equation (8) to approximate the convolution integral of this function, ${\displaystyle {\int }_0^{t}r\left(t-\tau \right)K\left(\tau \right)d\tau }$, we can have

$${\displaystyle {\int }_0^{t}r\left(t-\tau \right)K\left(\tau \right)d\tau }\approx {\displaystyle \sum _{k=0}^{[t{2}^j]}r\left(t_{[t{2}^j]-k}\right)}K\left(t_k\right){\Phi }_{j,k}^{\left(-1\right)}\left(t\right).$$

Setting $t=t_i=i/2^j$ = iΔt, we can obtain

$$u20\left(t_i\right)\approx A^{-1}{\displaystyle \sum _{k=0}^{i-1}r\left(t_{i-k}\right)\mathbf{K}\left(t_k\right){\Phi }_{m,k}^{\left(-1\right)}\left(t_i\right)}+A^{-1}\left(g\left(t_i\right)-q_1(t_i)\right).$$

(26)

in which the property r(0)=0 has been considered. It can be seen from Equation (26) that the solution u20 can be directly obtained step-by-step as the index i increases. Then the unknown function u can be reconstructed by the relation u = Au20 + q₁(t).

4. Numerical Results

In this section, we study the accuracy and stability of the proposed method by solving some examples of nonlinear wave equations. Error norms are defined as follows:

$$L_{\infty }\left(t\right)=\underset{1\le i\le N}{\max }\left|e_i\right|,L_2\left(t\right)=\sqrt{{\displaystyle \sum _{i=1}^N\left|{\left(e_i\right)}^2\right|}},\mathrm{RMS}\left(t\right)=L_2\left(t\right)/\sqrt{\Theta }$$

(27)

in which $e_i={\left(u_{\mathrm{exact}}\right)}_i-{\left(u_{\mathrm{approx}}\right)}_i$, u_approx and u_exact are exact and approximated solutions, respectively, and $\Theta$ the total number of grids in space.

Example 1.

We first consider the following time fractional diffusion-wave equation [8]

$$\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u}{\partial x^2}+\sin \left(\pi x\right),x\in [0,1]\\ u\left(x,0\right)=0,u_t\left(x,0\right)=0,\\ \left(0,t\right)=u\left(1,t\right)=0.\end{array}$$

(28)

The exact solution of this fractional differential equation is $u\left(x,t\right)=\frac{1}{{\pi }^2}\left[1-E_{\alpha ,1}\left(-{\pi }^2{t}^{\alpha }\right)\right]\sin \left(\pi x\right).$

We take Δh = 1/16 and Δt = 1/1024 to solve Equation (28). Comparison between the exact and approximate solutions of Equation (28) with different α are shown in Figure 2. From this figure, one can see that the numerical solution is consistent with the exact solution.

Table 1. Error norms of numerical solutions at t=1 for Example 1 with different α.

α	Scheme 1 [8] (Δh = 1/20 Δt = 1/10,000)	Scheme 2 [8] (Δh = 1/20 Δt = 1/10,000)	Scheme 1 [8] (Δh = 1/10,000 Δt = 1/1280)	Scheme 2 [8] (Δh = 1/10,000 Δt = 1/1280)	WICM (Δh = 1/16 Δt ≈ 1/1000)
1.3	1.7 × 10−4	1.7 × 10−4	1.4 × 10−5	3.9 × 10−6	3.2 × 10−6
1.5	2.3 × 10−4	2.3 × 10−4	1.8 × 10−5	1.2 × 10−5	2.2 × 10−6
1.8	3.2 × 10−4	3.2 × 10−4	1.6 × 10−4	6.0 × 10−5	1.8 × 10−7

shows the error norms of the numerical solutions obtained by the WICM and two finite difference schemes [8] named as Schemes 1 and 2. The standard backward Euler was used in the time discretization of the Schemes 1 and 2, while the central difference and Crank–Nicolson technique were used in the spatial discretization [8]. From Table 1, one can see that solutions obtained by using the present wavelet method with much coarser space–time meshes can have a much better numerical accuracy than those achieved from the finite difference.

Example 2.

We consider the following FWE with damping [9,10]

$$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}+\frac{\partial u\left(x,t\right)}{\partial t}=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+q\left(x,t\right),x\in \left[0,1\right],$$

(29)

whose exact solution is $u\left(x,t\right)=t^{2}x\left(1-x\right)$with the homogenous initial and boundary conditions and $q\left(x,t\right)=\frac{2x\left(1-x\right)}{\Gamma \left(3-\alpha \right)}{t}^{2-\alpha }+2tx\left(1-x\right)+2t^2.$

This problem is solved by the WICM with Δh = 1/16 and Δt = 1/1024. The approximate solution and its absolute error for α=1.3 are shown in Figure 3. We can see that the numerical solution obtained by the WICM is in good agreement with the analytical solution. The absolute errors of the approximate solutions for different α are shown in

Table 2. The absolute errors of wavelet solution of Example 2 for different α.

(x, t)	α = 1.1	α = 1.3	α = 1.5	α = 1.7	α = 1.9
(1/8, 1/8)	6.0 × 10−6	2.1 × 10−6	6.5 × 10−7	8.9 × 10−7	3.7 × 10−6
(2/8, 2/8)	1.6 × 10−5	5.1 × 10−6	1.4 × 10−6	1.9 × 10−6	1.0 × 10−5
(3/8, 3/8)	2.4 × 10−5	7.0 × 10−6	1.7 × 10−6	2.3 × 10−6	1.5 × 10−5
(4/8, 4/8)	2.8 × 10−5	7.5 × 10−6	1.6 × 10−6	1.9 × 10−6	1.5 × 10−5
(5/8, 5/8)	2.6 × 10−5	6.5 × 10−6	1.2 × 10−6	1.1 × 10−6	1.1 × 10−5
(6/8, 6/8)	2.0 × 10−5	4.7 × 10−6	7.5 × 10−7	3.9 × 10−7	5.3 × 10−5
(7/8, 7/8)	1.1 × 10−5	2.4 × 10−6	3.3 × 10−7	1.8 × 10−8	1.1 × 10−6

. From Figure 3 and Table 2, one finds that the WICM is very efficient and accurate in solving this problem.

Example 3.

We consider the 1D time fractional telegraph equation [11,12,13,14]

$$\begin{array}{l}\frac{{\partial }^{2\alpha }u\left(x,t\right)}{\partial t^{2\alpha }}+\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}+u\left(x,t\right)=\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right),\hspace{1em}\hspace{1em}x\in [0,1]\hspace{1em}\\ u\left(x,t\right)|{}_{t=0}=\frac{\partial u\left(x,t\right)}{\partial t}|{}_{t=0}=0,\end{array}$$

(30)

with the Dirichlet boundary conditions extracted from the exact solution $u\left(x,t\right)=t^2\sin \left(x\right).$

Figure 4 shows the comparison between the exact solution and the present wavelet solution with Δh = 1/16 and Δt = 1/1024 for α = 0.9. In addition,

Table 3. Error norms of wavelet solution at different times t for Example 3.

t	α = 0.7		α = 0.9
	L∞	RMS	L∞	RMS
0.25	6.4 × 10−6	1.8 × 10−5	1.6 × 10−5	4.1 × 10−5
0.5	5.7 × 10−6	1.6 × 10−5	1.2 × 10−5	3.1 × 10−5
0.75	5.0 × 10−6	1.4 × 10−5	2.5 × 10−6	7.0 × 10−6
1	2.1 × 10−5	4.5 × 10−5	4.3 × 10−6	1.2 × 10−5

gives the error norms of the wavelet solutions at several times for α = 0.7 and 0.9, respectively.

Example 4.

We consider the 1D time fractional Klein–Gordon and sine-Gordon equations [15]:

$$\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u}{\partial x^2}-f\left(u\right)+P(x,t),x\in \left[0,1\right]\\ u\left(x,0\right)=g_1\left(x\right),u_t\left(x,0\right)=g_2\left(x\right),\end{array}$$

(31)

with two types of nonlinear term f(u)

$$\begin{array}{l}\mathrm{case}1f\left(u\left(x,t\right)\right)=2u^3\hspace{1em}\hspace{1em}\hspace{1em}\left(\mathrm{Klein}-\mathrm{Gordon}\right),\\ \mathrm{case}2f\left(u\left(x,t\right)\right)=\sin \left(u\right)\hspace{1em}\left(\sin -\mathrm{Gordon}\right).\end{array}$$

(32)

The initial and Dirichlet boundary conditions, as well as the source item P(x, t) are compatible with the exact solution u(x, t) = sin(πx)(t4 + 1).

The error norms of the numerical solutions obtained, respectively, by the proposed WICM and linearized finite difference method (L-FDM) and the classical one (C-FDM) [15] are presented in

Table 4. Error norms of numerical solutions with similar time step for Example 4 with different α.

α	Case 1			Case 2
	L-FDM [15] ∆h = 1/80 Δt = 0.001	C-FDM [15] ∆h= 1/160 Δt = 0.001	WICM ∆h = 1/16 Δt ≈ 0.001	L-FDM [15] ∆h = 1/80 Δt = 0.001	C-FDM [15] ∆h = 1/160 Δt = 0.001	WICM ∆h = 1/16 Δt ≈ 0.001
1.2	6.7 × 10−5	1.7 × 10−5	4.0 × 10−5	1.5 × 10−4	3.5 × 10−5	2.6 × 10−5
1.5	8.1 × 10−5	2.0 × 10−5	7.9 × 10−5	1.4 × 10−4	3.3 × 10−5	1.1 × 10−5
1.8	1.0 × 10−4	2.5 × 10−5	1.2 × 10−5	1.6 × 10−4	4.0 × 10−5	1.5 × 10−5

and

Table 5. Error norms of numerical solutions with similar spatial grid for Example 4 with different α.

α	Case 1				Case 2
	L-FDM [15] ∆t = 1/1000 ∆h = 1/20	C-FDM [15] ∆t = 1/5000 ∆h = 1/20	WICM ∆t = 1/512 ∆h = 1/16	WICM ∆t = 1/4096 ∆h = 1/16	L-FDM [15] ∆t = 1/1000 ∆h = 1/20	C-FDM [15] ∆t = 1/5000 ∆h = 1/20	WICM ∆t = 1/512 ∆h = 1/16	WICM ∆t = 1/4096 ∆h = 1/16
1.2	1.1 × 10−3	3.3 × 10−6	8.9 × 10−5	2.2 × 10−6	2.4 × 10−3	7.1 × 10−6	5.9 × 10−5	4.8 × 10−6
1.5	1.3 × 10−3	3.9 × 10−6	2.2 × 10−5	8.9 × 10−7	2.3 × 10−3	6.8 × 10−6	1.8 × 10−6	8.4 × 10−8
1.8	1.7 × 10−3	5.0 × 10−6	4.3 × 10−5	1.0 × 10−7	2.7 × 10−3	8.0 × 10−6	1.9 × 10−6	4.3 × 10−8

. The second-order and fourth-order difference schemes in spatial discretization are used in the L-FDM and C-FDM, respectively [15]. As shown in Table 4, when the time step is Δt ≈ 0.001, the WICM with much coarser spatial grid can achieve even better accuracy than both of the L-FDM and C-FDM. Moreover, when the spatial grid sizes are ∆h = 1/16 for the WICM and 1/20 for the L-FDM and C-FDM, the former also can use a much larger time step to obtain a better accuracy compared with the latter, as shown in Table 5.

Example 5.

We consider the 2D time fractional Klein–Gordon and sine-Gordon equations [27,46]:

$$\{\begin{array}{l}\frac{{\partial }^{\alpha }u\left(x,y,t\right)}{\partial t^{\alpha }}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}-f\left(u\right)+P(x,y,t),\left(x,y\right)\in {[0,1]}^2,\\ u\left(x,y,0\right)=0,u_t\left(x,y,0\right)=0,\end{array}$$

(33)

where we consider two cases

$$\begin{array}{c}\mathrm{Case}1:f\left(u\right)=u^3,P\left(x,y,t\right)=\left[\frac{2}{\Gamma \left(3-\alpha \right)}{t}^{2-\alpha }+2t^2\right]\sin \left(x+y\right)+{\left[t^2\sin \left(x+y\right)\right]}^3;\\ \mathrm{Case}2:f\left(u\right)=\sin \left(u\right),P\left(x,y,t\right)=\left[\frac{2}{\Gamma \left(3-\alpha \right)}{t}^{2-\alpha }+2t^2\right]\sin \left(x+y\right)+\sin \left[t^2\sin \left(x+y\right)\right].\end{array}$$

(34)

From Equations (33) and (34), we can easily to verify that these two cases have the same exact solution

u(x, t) = t²sin(x + y).

(35)

The first case refers to the Klein–Gordon equation, and the second case refers to that of the sine-Gordon equation.

Table 6. Error norms of the wavelet solution at t = 1 with various time steps for Example 6.

	∆t(∆h = 1/16)	α = 1.25		α = 1.75
		L∞	RMS	L∞	RMS
Case 1	1/28	2.6 × 10−4	1.9 × 10−4	6.9 × 10−5	4.0 × 10−5
	1/210	4.7 × 10−5	3.4 × 10−5	1.1 × 10−5	6.2 × 10−6
	1/212	8.5 × 10−6	6.0 × 10−6	1.8 × 10−6	1.0 × 10−6
Case 2	1/28	2.6 × 10−4	1.9 × 10−4	7.4 × 10−5	4.3 × 10−5
	1/210	4.7 × 10−5	3.4 × 10−5	1.2 × 10−5	6.7 × 10−5
	1/212	8.5 × 10−6	6.0 × 10−6	1.9 × 10−6	1.1 × 10−6

and

Table 7. Error norms of the wavelet solution at t = 1 with various spatial grid sizes for Example 6.

	∆h	α = 1.3	α = 1.5	α = 1.7	α = 1.9
		L∞	L∞	L∞	L∞
Case 1	1/16	2.0 × 10−4	6.9 × 10−5	5.9 × 10−5	1.4 × 10−4
	1/32	3.4 × 10−5	8.6 × 10−6	8.1 × 10−6	2.9 × 10−5
	1/64	5.7 × 10−6	1.1 × 10−6	1.2 × 10−6	6.4 × 10−6
Case 2	1/16	2.0 × 10−4	6.9 × 10−5	5.5 × 10−5	1.4 × 10−4
	1/32	3.4 × 10−5	8.6 × 10−6	7.6 × 10−6	2.2 × 10−5
	1/64	5.7 × 10−6	1.1 × 10−6	1.1 × 10−6	4.8 × 10−6

show the relationships between the error norms of the proposed wavelet solutions and the grid sizes. It can be seen that the error decreases rapidly as the spatial grid size or the time step decreases, implying good properties in convergence and stability of the WICM. Table 6 and Table 7 also demonstrate that the proposed WICM can have better accuracy than several existing methods under the more coarse grid. For example, the maximum absolute error of the wavelet solution with ∆h = 1/16 and ∆t = 1/256 is L_∞ = 2.6 × 10⁻⁴ for the Klein–Gordon problem with α = 1.25, which is much smaller than L_∞=1.4288 × 10⁻³ and 6.1192 × 10⁻⁴ for the RBF meshless approach [27] and local meshless method [46] with ∆h = 1/21 and ∆t = 1/320, respectively.

5. Conclusions

An efficient WICM method based on the generalized Coiflet wavelet is proposed to solve the time FWEs with various nonlinearities. By using the Laplace transform, the original FWE is first transformed into the time Volterra-type integro-differential equation. Then, the WICM is developed to solve the resulting integro-differential equation. Such an approach has successfully avoided the difficulty in the treatment of singularity associated with the fractional derivatives, and the time-consuming task of solving the nonlinear spatial system at each time step. Numerical results of five benchmark examples show that, compared with several existing methods, the proposed wavelet method can obtain more accurate solutions in the case of using a coarser grid, which reflects the excellent characteristics of the WICM in terms of convergence and stability.

For possible future studies, we can consider the applicability of the proposed WICM in the solution of nonlinear wave equations with fractional derivatives under different definitions. We can also combine the wavelet approximation algorithm of functions bounded on irregular domains [47,48] into the WICM to solve real engineering problems with irregular spatial shapes/domains.