Hybrid Wavelet Methods for Nonlinear Multi-Term Caputo Variable-Order Partial Differential Equations

1  Introduction

Wavelet theory has emerged as a powerful tool in numerical analysis due to its inherent ability to represent functions with localized features. Wavelets are a special kind of function that forms the basis of L2(R) and can be generated by dilating and translating a mother wavelet. They have a compact support that is, wavelets are zero outside a finite interval. This is the reason that wavelets-based operational matrices contain many zero entries. Researchers have worked on different wavelets for solving different types of fractional differential equations. Awati et al. [1] utilized the Haar wavelet-based collocation method and shifted Chebyshev collocation method for the simulations of magnetohydrodynamic flow of micropolar nanofluid with stagnation point model. Rabiei and Razzaghi [2] solve the fractional optimal control problems by using the fractional-order Boubaker wavelets method. Bernoulli wavelet operational matrix method [3] is used to solve the proposed glucose-insulin regulatory model. Neutral delay differential equations are solved by using the Gegenbauer and Bernoulli wavelets method in [4]. For more detail about the wavelet methods, we refer the readers to [5–7] and [8–11].

Parallel to this, fractional calculus—particularly with constant-order derivatives—has gained considerable attention over the past decade, owing to its capacity to model memory and hereditary properties in various scientific and engineering systems. Constant-order fractional calculus is the study of various types of fractional-order differential and integral operators. Variable-order (VO) fractional calculus is the generalization of the constant-order fractional calculus in which the order of its operators is a function of time or space or both. In recent years, VO differential equations have gained more interest among researchers due to their many applications in real world phenomena. The VO differential equations can more accurately describe the complex behavior of many physical phenomena. They enable a more precise depiction of nonlinear phenomena or systems experiencing abrupt changes. Like fractional differential and integral operators, there are many definitions of variable-order differential and integral operators with singular and non-singular kernel.

Despite the evident advantages of VO models, they have received comparatively less attention than their constant-order counterparts. We have reviewed several research articles from the literature which are focussed on the VO fractional calculus. For example, a Bernoulli polynomials-based method is proposed in [12] for solving multi-term VO ordinary differential equations in Caputo sense. Legendre wavelet based method is utilized in [13] for the solution of nonlinear VO ordinary differential equations. In [14], an explicit finite difference scheme is utilized for the solution of linear and semi-linear VO differential equations. The focus of the authors in [15] is to study the existence and uniqueness of weak solutions of VO Laplacian equations with variable exponents. In [16], authors present a method for the solution of Atangana-Baleanu VO mobile-immobile advection-dispersion model. Shen et al. [17] solve the VO time fractional diffusion equation by utilizing the L1 approximation for CVO fractional time derivative and the finite difference formula for the second-order spatial derivative. They have also worked on the analysis of the method.

Since the kernel of the VO fractional operators has a variable exponent, obtaining analytical solutions is difficult, and these have not gained much attention from the researchers. Building on these foundations, this paper aims to develop a robust numerical method for solving VO fractional differential equations by leveraging wavelet-based techniques. By integrating the localized efficiency of wavelets with the modeling flexibility of variable-order operators, we aim to contribute a more accurate and computationally efficient approach for complex nonlinear multi-term CVO partial differential equation of following form

DCtγ(s,t)u(s,t)=∑j=1Jaj(s,t)DCsβj(s,t)u(s,t)+c(s,t)u(s,t)+g(u,DCsα(s,t)u)+f(s,t),u(s,0)=u0(s), 0≤s≤η, η∈R+,u(0,t)=A0(t), u(η,t)=C0(t), 0≤t≤T, T∈R+,(1)

where 0<nγ≤γ(s,t)≤mγ≤1, 0<α(s,t),βj(s,t)≤2, J∈N, and g(u,DCsα(s,t)u) is any nonlinear functional of u and DCsα(s,t)u. Here, DCtγ(s,t)u is a CVO fractional time derivative defined by

DCtγ(s,t)u(s,t)=1Γ(1−γ(s,t))∫0t(t−τ)−γ(s,t)∂∂τu(s,τ)dτ.(2)

The CVO fractional space derivatives, DCsβj(s,t)u(s,t), are defined similarly. We have reviewed several research articles related to the variable-order partial differential equations, we will mention a few of them here. We found that most of the authors [18–24] focussed their attention on only γ(t), t∈ (0,1), order time derivative but integer order space derivatives partial differential equations. Some authors [25–28] consider the γ(s,t), 0≤s≤η, t∈ (0,1), order time derivative but integer order space derivatives. In [29,30], authors focussed their intention on the solution of variable-order partial differential equations in both VO space, β(s), and time, γ(t), derivative. Since the above problem (1) is a nonlinear multi-term CVO fractional partial differential equation in which the orders of both space and time derivatives are the functions of both time and spatial variables. Numerical techniques for solving these Eq. (1) are still in the early stages of development. This is the reason that we focussed our intention on the development of more efficient and accurate methods for the solution of (1). Our major contributions to this work are as follows:

•   Zhang et al. [18] developed an exponential-sum-approximation (ESA) technique for γ(t), t∈ (0,1), order time derivative, and utilized it for the solution of diffusion equations with VO time and integer order spatial derivative. We extended the ESA technique for γ(s,t), 0≤s≤η, 0≤t≤T. The purpose of this extension is to get a fast approximation of DCtγ(s,t)u(s,t).

•   To the best of our knowledge, this work is the first to propose a hybrid strategy for Caputo-type fractional differential equations with variable-order derivatives in both time and space, combining an efficient ESA-based scheme for the time component with CAS wavelet operational matrices for handling nonlinear space variable-order derivatives.

•   To handle the nonlinearities, we designed the interpolation technique for the above problem (1). The purpose of using the interpolation technique for the treatment of nonlinear terms is to reduce the computational costs of the methods as compared to the quasilinearization technique and Adomian polynomials.

•   We proposed two efficient method, the L1-CAS method and the fast-CAS method, for the solution of (1), by uniting the L1 approximations or the ESA technique for DCtγ(s,t)u(s,t), CAS wavelet variable-order operational matrices for CVO space derivatives, DCsβ1(s,t)u(s,t), DCsβ2(s,t)u(s,t),⋯, DCsβJ(s,t)u(s,t), and the interpolation technique for the treatment of nonlinearities, given in Section 3.

•   We performed the theoretical analysis of the proposed methods and provided the numerical simulations to illustrate the theoretical results. The obtained numerical results are thoroughly examined through both tabular and graphical representations.

This paper is organized as follows: Section 2 provides the detail about the CAS wavelet, the function approximations by the CAS wavelet series, and the construction of its operational matrices of variable-order integration. In Section 3, we discuss the algorithms of the L1-CAS method, the fast-CAS method, and their extension to the nonlinear multi-term CVO partial differential equation. Section 4 is dedicated to conducting the analysis of the proposed methods. In Section 5, we present two numerical examples aimed at evaluating the efficiency, reliability, and accuracy of our methods. In Section 6, we conclude our work.

2  CAS Wavelets

The cosine and sine (CAS) wavelets represent a distinct category of wavelet basis functions, originating from the cosine and sine functions. The CAS wavelets possess orthogonality to each other, streamlining computations and improving function representation. The classical CAS wavelets for interval [0,1) is defined as [31]

ψp,q(s)={2ν2CASq(2νs−p),p2ν≤s<p+12ν,0,elsewhere,

The CAS wavelets for interval [0,η) is given as

ψp,q(s)={2νηCASq(2νηs−p),ηp2ν≤s<ηp+12ν,0,elsewhere,(3)

where CASq(s)=cos⁡(2qπs)+sin⁡(2qπs), and ν=0,1,2,⋯, and p=0,1,2,⋯,2ν−1, is the level of resolution and translation parameter, respectively, q∈Z.

These wavelets for all integers ν and p produce an orthonormal basis of L2(R). These wavelets have compact support, that is, they are zero outside [ηp−12ν,ηp2ν). These wavelets are valuable for approximating solutions to differential equations, especially in scenarios where conventional methods may be computationally demanding or unfeasible.

The CAS wavelets are orthonormal, that is,

∫0ηψp,q(s)ψp,q′(s)ds=δqq′,

which can be proved by using the transformation t=sη in [31,32].

2.1 Function Approximations

Since the CAS wavelets for all ν and p forms an orthonormal basis of L2[0,η), therefore, any function u(s)∈L2[0,η) can be represented by the infinite series of the CAS wavelets as

u(s)=∑p=0∞∑q∈Zlp,qψp,q(s).

We can use the finite sum of basis functions, ψp,q(s), to approximate u(s) as

u(s)≈∑p=02ν−1∑q=−QQlp,qψp,q(s)=LTΨ(s),

where M=2ν(2Q+1), Q∈N, and L and Ψ(s) are vectors of length M×1, which are given below:

L=[l0,−Q,l0,−Q+1,⋯,l0,Q,l1,−Q,l1,−Q+1,⋯,l1,Q,⋯,l2ν−1,−Q,l2ν−1,−Q+1,⋯,l2ν−1,Q]TΨ(s)=[ψ0,−Q(s),ψ0,−Q+1(s),⋯,ψ0,Q(s),ψ1,−Q(s),ψ1,−Q+1(s),⋯,ψ1,Q(s),⋯,ψ2ν−1,−Q(s),ψ2ν−1,−Q+1(s),⋯,ψ2ν−1,Q(s)]T.

Expand the truncated series of the CAS wavelets at the collocation points, sk=η2k−12M, k=1,2,⋯,M, we get the CAS wavelets matrix ΨM×M as

ΨM×M=[Ψ(s1),Ψ(s2),⋯,Ψ(sM)].

For ν=1, Q=1, η=2, we get

Ψ6×6=(−0.5176−1.41421.9319000000−0.5176−1.41421.93191.41421.41421.41420000001.41421.41421.41421.9319−1.4142−0.51760000001.9319−1.4142−0.5176).

2.2 Operational Matrices

In this section, our focus is on the construction of CAS wavelet operational matrices of variable-order fractional integration, which enable us to transform the CVO differential equation into a matrix form. Also, the purpose of constructing the operational matrices is to make the calculations fast, because operational matrices contain many zero entries.

The CAS wavelets operational matrix of variable-order integration

Let η0=ηp−12ν and η1=ηp2ν. Apply the fractional integral operator of order γ(s)>0 on the CAS wavelets as

I0sγ(s)ψp,q(s)=1Γ(γ(s))∫0s(s−τ)γ(s)−1ψp,q(τ)dτ, η0≤s<η1,=2νηΓ(γ(s))∫η0s(s−τ)γ(s)−1CASq(2νητ−(p−1))dτ,=2νηΓ(γ(s))(∫η0s(s−τ)γ(s)−1cos⁡(2ν+1qπητ−2qπ(p−1))dτ+∫η0s(s−τ)γ(s)−1sin⁡(2ν+1qπητ−2qπ(p−1))dτ).(4)

Let us denote  I0sγ(s)ψp,q(s)=Fp,qγ(s)(s), p=0,1,2,⋯,2ν−1, q={−Q,−Q+1,⋯,Q}, Q∈N. We can compute (4) by using the global adaptive technique command, “integral”, at different values of the parameters p and q in MATLAB. In vector form, we can represent (4) as

Fγ(s)(s)=[F0,−Qγ(s)(s),F0,−Q+1γ(s)(s),⋯,F0,Qγ(s)(s),F1,−Qγ(s)(s),F1,−Q+1γ(s)(s),⋯,F1,Qγ(s)(s),⋯,F2ν−1,−Qγ(s)(s),F2ν−1,−Q+1γ(s)(s),⋯,F2ν−1,Qγ(s)(s)]T.(5)

To get the CAS wavelets operational matrix of variable-order integration, we will expand (5) at the collocation points, sk=η2k−12M, k=1,2,⋯,M, as

FM×Mγ(s)=[Fγ(s1)(s1),Fγ(s2)(s2),⋯,Fγ(sM)(sM)].(6)

For γ(s)=sin⁡(s)+exp⁡(−s)−12, ν=1, Q=1, η=2, we obtain

F6×6γ(s)=(0.0769−0.8856−0.04350.12690.06380.04790000.1000−0.72680.02650.63681.05621.38311.31291.13500.89730000.41730.94031.40720.85390.2597−0.3951−0.0572−0.0407−0.03380000.55270.3749−0.4454).

The CAS wavelets operational matrix of variable-order integration for boundary value problems

To effectively solve the boundary value problems, it is mandatory to utilize another important operational matrix of variable order integration. To get this matrix, let us apply the γ(s) order integral on (3) from 0 to η as

I0ηγ(s)ψp,q(s)|s=η=1Γ(γ(η))∫0η(η−τ)γ(η)−1ψp,q(τ)dτ,=2νηΓ(γ(η))(∫η0η1(η−τ)γ(η)−1cos⁡(2ν+1qπητ−2qπ(p−1))dτ+∫η0η1(η−τ)γ(η)−1sin⁡(2ν+1qπητ−2qπ(p−1))dτ).(7)

Let us denote I0ηγ(s)ψp,q(s)|s=η=bp,qγ(η), γ(η)>0, p=0,1,2,⋯,2ν−1, q=[−Q,−Q+1,⋯,0,1,⋯,Q−1,Q], Q∈N. In vector form, we can represent (7) as

BM×1γ(η)=[b0,−Qγ(η),b0,−Q+1γ(η),⋯,b0,Qγ(η),b1,−Qγ(η),b1,−Q+1γ(η),⋯,b1,Qγ(η),⋯,b2ν−1,−Qγ(η),b2ν−1,−Q+1γ(η),⋯,b2ν−1,Qγ(η)]T.

Since BM×1γ(η) does not depend on any variable, this means that for sk=η2k−12M, k=1,2,⋯,M, we get same BM×1γ(η) at each collocation points. Therefore, we have

BM×Mγ(η)=[BM×1γ(η),BM×1γ(η),⋯,BM×1γ(η)].(8)

For γ(η)=sin⁡(η)+exp(−η)−12, ν=1, Q=1, η=2, we have

B6×1γ(η)=[0.0415, 0.5879, 0.7300, 1.5917, −0.0300, 0.0752].

3  Algorithms

This section is devoted to the development of two methods, one is based on the fast algorithm and the CAS wavelet technique, and the second is based on the L1-approximations and the CAS wavelet technique for the solution of linear and nonlinear multi-term CVO fractional partial differential equations.

Consider the following general form of linear multi-term CVO fractional partial differential equation

DCtγ(s,t)u(s,t)=∑j=1Jaj(s,t) DCsβj(s,t)u(s,t)+c(s,t)u(s,t)+f(s,t),u(s,0)=u0(s), 0≤s≤η,u(0,t)=A0(t), u(η,t)=C0(t), 0≤t≤T,(9)

where  0<βj(s,t)≤2, 0<γ(s,t)≤1, J∈N,DCtγ(s,t)u(s,t) represents the CVO time derivative of order γ(s,t) of u(s,t). In this work, we will evaluate DCtγ(s,t)u(s,t) by the L1-approximation or the fast algorithm. For CVO space derivatives, DCsβ1(s,t)u(s,t), DCsβ2(s,t)u(s,t),⋯, DCsβJ(s,t)u(s,t), we will utilize the CAS wavelet technique.

3.1 The Fast-CAS Method

In this subsection, we will work to save the memory and computational time by constructing an efficient algorithm for the CVO time derivative, DCtγ(s,t)u(s,t), 0<γ(s,t)≤1, and utilizing the CAS wavelet operational matrices technique. Discretize the time domain [0,T], T∈R+ as tr=rΔt, r=0,1,2,⋯,R∈N, and Δt=TR. For space domain [0,η], η∈R+, we have sk=η2k−12M, k=1,2,⋯,M. Consider the CVO time derivative (2) at (sk,tr) and split it in terms of local and history part, we have

DCtγ(sk,tr)u(sk,tr)=Hk,r+Lk,r,(10)

where Hk,r=1Γ(1−γ(sk,tr))∫0tr−1(tr−τ)−γ(sk,tr)∂∂τu(sk,τ)dτ and Lk,r=1Γ(1−γ(sk,tr))∫tr−1tr(tr−τ)−γ(sk,tr)∂∂τu (sk,τ)dτ are called the history and local part, respectively. For Lk,r, we utilize the forward difference formula for the approximation of ∂∂τu(sk,τ) on [tr−1,tr], and integrate the remaining part to have

Lk,r≈u(sk,tr)−u(sk,tr−1)Δtγ(sk,tr)Γ(2−γ(sk,tr)).

For Hk,r, we have

Hk,r=T−γ(sk,tr)Γ(1−γ(sk,tr))∫0tr−1(tr−τT)−γ(sk,tr)∂∂τu(sk,τ)dτ.(11)

The kernel (tr−τT)−γ(sk,tr) can be approximated by using the Lemma 2 of [18] as

|(tr−τT)−γ(sk,tr)−∑i=Nl+1NUΩk,r,ie−ξi(tr−τT)|≤(tr−τT)−γ(sk,tr)ϵ

where the quadrature weights and exponents are properly chosen by following the procedure [18,33], and are given as

Ωk,r,i=heγ(sk,tr)ihΓ(γ(sk,tr)), ξi=eih, ϵ=(ΔtT)2,h=2πlog⁡(3)+mγlog⁡(1/cos⁡1)+log⁡(1/ϵ),Nl=⌈1hnγ(log⁡(ϵ)+log⁡(Γ(1+mγ)))⌉,NU=⌊1h(log⁡(TΔt)+log⁡(log⁡(1/ϵ))+log⁡(nγ)+0.5)⌋,(12)

where 0<nγ≤γ(s,t)≤mγ≤1. Then the Eq. (11) can be written as

Hk,r≈T−γ(sk,tr)Γ(1−γ(sk,tr))∑i=Nl+1NUΩk,r,iVk,r,i,

where Vk,r,i=∫0tr−1eξi(tr−τT)∂∂τu(sk,τ)dτ. For r=1, Vk,1,i=0, and for r=2,3,⋯,R,Vk,r,i can be calculated by using the forward difference formula for ∂∂τu(sk,τ), as

Vk,r,i=∫0tr−2eξi(tr−τT)∂∂τu(sk,τ)dτ+∫tr−2tr−1eξi(tr−τT)∂∂τu(sk,τ)dτ,≈e−ξiΔtTVk,r−1,i+u(sk,tr−1)−u(sk,tr−2)Δt∫tr−2tr−1eξi(tr−τT)dτ,=e−ξiΔtTVk,r−1,i+Te−ξiΔtT−e−2ξiΔtTξiΔt(u(sk,tr−1)−u(sk,tr−2)).(13)

Substitute Hk,r and Lk,r in (10) to get the fast evaluation formula for the CVO time derivative of order γ(s,t) as

DFtγ(sk,tr)u(sk,tr)=T−γ(sk,tr)Γ(1−γ(sk,tr))∑i=Nl+1NUΩk,r,iVk,r,i+u(sk,tr)−u(sk,tr−1)Δtγ(sk,tr)Γ(2−γ(sk,tr)),(14)

where Vk,r,i is given in (13). For r=1, we have

DFtγ(sk,t1)u(sk,t1)=u(sk,t1)−u(sk,t0)Δtγ(sk,t1)Γ(2−γ(sk,t1)).(15)

Let er(s)=1Δtγ(s,tr)Γ(2−γ(s,tr)),ajr(s)=aj(s,tr) and dr(s)=c(s,tr)−er(s). Utilize the fast algorithm (14) for the evaluation of Caputo variable-order time derivative, DCtγ(s,t)u(s,t), at t=tr, r=1,2,3,⋯,R, in Eq. (9), we get

∑j=1Jajr(s)DCsβj(s,tr)u(s,tr)+dr(s)u(s,tr)=−f(s,tr)−er(s)u(s,tr−1)+T−γ(s,tr)Γ(1−γ(s,tr))∑i=Nl+1NUΩs,r,iVs,r,i,u(s,t0)=u0(s), u(0,tr)=A0(tr), u(η,tr)=C0(tr),(16)

where Vs,r,i=e−ξiΔtTVs,r−1,i+Te−ξiΔtT−e−2ξiΔtTξiΔt(u(s,tr−1)−u(s,tr−2)). Let

zjr(s)={s1−βj(s,tr)ηΓ(2−βj(s,tr)),0<βj(s,tr)≤1,0,1<βj(s,tr)≤2.

Apply the CAS wavelet technique on the semi discretized Eq. (16), along the boundary conditions. Firstly, we will approximate DCs2u(s,tr) by using the CAS wavelet series as

DCs2u(s,tr)≈∑p=02ν−1∑q=−QQlp,qrψp,q(s).(17)

Apply the Riemann-Liouville integral operator of order 2 on (17) and utilize the boundary conditions, we get

u(s,tr)≈∑p=02ν−1∑q=−QQlp,qr I0s2ψp,q(s)−sη∑p=02ν−1∑q=−QQlp,qr I0η2ψp,q(η)+sη(C0(tr)−A0(tr))+A0(tr).(18)

Now, apply the Caputo derivative of order βj(s,tr) on (18) to get

DCsβj(s,tr)u(s,tr)≈∑p=02ν−1∑q=−QQlp,qr I0s2−βj(s,tr)ψp,q(s)−zjr(s)∑p=02ν−1∑q=−QQlp,qr I0η2ψp,q(η)+zjr(s)(C0(tr)−A0(tr)).(19)

Let er′(s)=sηdr(s),ojr(s)=ajr(s)zjr(s) and gJr(s)=−(∑j=1Jajr(s)zjr(s)+er′(s))(C0(tr)−A0(tr))−dr(s)A0(tr)−f(s,t)−er(s)u(s,tr−1)+T−γ(s,tr)Γ(1−γ(s,tr))∑i=Nl+1NUΩs,r,iVs,r,i. Put the values of DCsβj(s,tr)u(s,tr) and u(s,tr) from Eqs. (19) and (18), respectively, in (16) to get

∑p=02ν−1∑q=−QQlp,qr(∑j=1Jajr(s) I0s2−βj(s,tr)ψp,q(s)−(∑j=1Jojr(s)+er′(s)) I0η2ψp,q(η)+dr(s) I0s2ψp,q(s))=gJr(s).(20)

Let yjr(s)=ojr(s)+er′(s)J. Evaluate (20) at the collocation points to obtain

LrT(∑j=1JAjrFM×M2−βj(s)+DrFM×M2−∑j=1JYjrBM×M2)=GJr,(21)

where GJrM×1=[gJr(η12M),gJr(η32M),⋯,gJr(η2M−12M)]T, and the matrices FM×M2−βj(s,tr),FM×M2, and BM×M2 are given in Section 2.2. The diagonal matrices Ajr,Dr, and Yjr are given below:

Ajr=[ajr(s1)0⋯00ajr(s2)⋯0⋮⋮⋱⋮00⋯ajr(sM)], Dr=[dr(s1)0⋯00dr(s2)⋯0⋮⋮⋱⋮00⋯dr(sM)],Yjr=[yjr(s1)0⋯00yjr(s2)⋯0⋮⋮⋱⋮00⋯yjr(sM)].

Solve the system (21) to get LrT at each value of tr and put it in (18) to get the approximate solution u(s,tr) at t=tr, r=1,2,⋯,R. See Algorithm 1 for details.

3.2 The L1-CAS Method

To get the L1 approximation of (2), we evaluate (2) at (s,tr) as

DCtγ(s,tr)u(s,tr)=1Γ(1−γ(s,tr))∫0tr(tr−τ)−γ(s,tr)∂∂τu(s,τ)dτ.

We will utilize the forward difference approximation to ∂∂τu(s,τ) over the interval [tj,tj+1], 0≤j≤r−1 as ∂∂τu(s,τ)≈u(s,tj+1)−u(s,tj)Δt, we obtain

DL1tγ(s,tr)u(s,tr)=1Γ(1−γ(s,tr))∑j=0r−1u(s,tj+1)−u(s,tj)Δt∫tjtj+1(tr−τ)−γ(s,tr)dτ,=Δt−γ(s,tr)Γ(2−γ(s,tr))∑j=0r−1Ej,r(u(s,tj+1)−u(s,tj)),or=Δt−γ(s,tr)Γ(2−γ(s,tr))(Er−1,ru(s,tr)−∑j=1r−1(Ej,r−Ej−1,r)u(s,tj)−E0,ru(s,t0)),(22)