# A Numerical Method for Solving Fractional Differential Equations

Authors: Yahong Wang; Haili Zhou; Liangcai Mei; Yingzhen Lin
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, we solve the fractional differential equations (FDEs) with boundary value conditions in Sobolev spaceHn0,1. The strategy is constructing multiscale orthonormal basis for Hn0,1 to get the approximation for the problems. The convergence of the method is proved, and it is tested on some numerical experiments; the tests show that our method is more efficient and accurate. The notion of numerical stability with respect to the condition number is introduced proving that the proposed method is numerically stable in this sense.


---

## Body

## 1. Introduction

FDEs have been a subject of interest not only among mathematicians but also among physicists and engineers. In fact, we can find numerous applications in economic system [1], fluid mechanics [2], biology [3], signal processing, dynamics of earthquakes, optics, electromagnetic waves, chaotic dynamics, polymer science, proteins, statistical physics, thermodynamics, neural networks, and so on [4]. A completely different and very novel application field is the area of mathematical psychology where fractional-order systems may be used to model the behaviour of human beings [5].FDEs are always with weakly singular kernel and more complicated than integer ones. Actually, in many cases, it is difficult to obtain the analytical solution. Therefore, the numerical methods are essential for approximation solution of many FDEs. Many approximations have sprung up recently, such as the numerical scheme by coupling of radial kernels and localized Laplace transform [6], the radial basis function (RBF)-based numerical scheme which uses the Coimbra variable time fractional derivative of order 0<αt,x<1 [7], the time-space numerical technique based on time-space radial kernels [8], the local meshless method based on Laplace transform [9], finite difference methods [10], finite element methods [11–13], spectral methods [14], shooting method [15], approximation formula [16], pseudo-spectral method [17], variational iteration method [18], Adomian decomposition method (ADM) [19, 20], trapezoidal methods [21], reproducing kernel method [22–26], and so on. The development, however, for efficient numerical methods to solve linear and nonlinear FDEs is still an important issue.In this paper, we are concerned with the approximation of FDEs as follows:(1)Dαux+a1xu′x+a0xux=fx,x∈0,1,α>0, (2)u0−α0u′0=γ1,u1+α1u′1=γ2.Here,(3)Dαux=1Γm−α∫0xx−tm−α−1umtdt,where m:=⌈α⌉,aixi=0,1,fx∈L20,1. The existence and uniqueness of the solution ux of problem (1) are established in Ref. [11].The following parts are structured as follows. In Section2, we consider the solution of (1) in Sobolev space, introduce a suitable definition of ε-approximate solution, and analyze the stability based on some assumptions. In Section 3, a set of orthonormal basis in H2 has been obtained which is used for constructing the ε-approximate solution unx which converges to the exact solution ux uniformly. In Section 4, the accuracy and convergence rate computed by the method we proposed show that our method is much more better, and the numerical simulations further verify our theoretical analysis. In the last section, we draw the conclusions.


## 2. Theoretical Analysis of the Approximate Solution for (1)

In order to present our algorithm, we consider solving operator equation (4) in Sobolev space as follows. Let ℒ:X⟶Y be linear and bounded, X be Sobolev space, Y=L2, u∈X,f∈Y,andℬii=1,2 be a bounded linear functional on X.(4)ℒu=f,ℬiu=γi.In the following analysis, we will use the operator normsℒ,ℬi [22]. Let ρ1,ρ2,…,ρk,… be a set of orthonormal basis for X. Combining the least square theory and the idea of residuals in the Sobolev space, we proposed the concept of ε-approximate solution for (4) as follows.Definition 1.
For∀ε>0, if ℒu−fY2+∑i=12ℬiu−γi2<ε2, u is called ε-approximate solution of (4).Theorem 1.
For∀ε>0,∃Nε, when n>Nε,(5)un=∑k=1nck∗ρk,is the ε-approximate solution for (4), whereck∗ satisfy(6)∑k=1nck∗ℒρk−fY2+∑i=12∑k=1nck∗ℬiρk−γi2=minck∑k=1nckℒρk−fY2+∑i=12∑k=1nckℬiρk−γi2.Proof.
Sinceu∈X is the solution of (4), ρ1,ρ2,…,ρk,… is a set of orthonormal basis for X, so u=∑k=1∞dkρk satisfy equation (4). For ∀ε>0,∃Nε, when n>Nε,v=∑k=1ndkρk satisfy v−uX<ε/ℒ2+∑i=12ℬi2.
In fact,(7)ℒv−fY2+∑i=12ℬiv−γi2=ℒv−ℒuY2+∑i=12ℬiv−ℬiu2≤ℒ2v−uX2+∑i=12ℬi2ℒ2v−uX2=ℒ2+∑i=12ℬi2v−uX2<ε2.
So, equation (6) ≤∑k=1ndkℒρk−fY2+∑i=12∑k=1ndkℬiρk−γi2<ε2.
From equation (6), we learn that ck∗ are the representations of ck which make J reach the minimum, where(8)Jc1,c2,…,cn=∑k=1nckℒρk−fY2+∑i=12∑k=1nckℬiρk−γi2.
Applying∂Jc1,c2,…,cn/∂cj=0 to get the minimum values, we obtain(9)∑k=1nckℒρk,ℒρjY+∑i=12∑k=1nckℬiρkℬiρj=ℒρj,fY+∑i=12ℬiρjγi.j=1,2,…,n.Theorem 2.
The solution of equation (9) is unique when ℒ is reversible.Proof.
We just prove that the solution of the homogeneous linear equations of (9) is unique. Multiply cj on both sides of the equation:(10)∑k=1nckℒρk,ℒρjY+∑i=12∑k=1nckℬiρkℬiρj=0.
Adding all equations together, we get(11)∑k=1nckℒρk,∑j=1ncjℒρjY+∑i=12∑k=1nckℬiρk∑j=1ncjℬiρj=0,which implies that(12)∑k=1nckℒρkY2=0.
Asℒ is reversible and ρkk=1n is a set of orthonormal basis, ck=0k=1,2,…,n.
Theorem1 and Theorem 2 show that the ε-approximate solution (5) exists and is unique as n>Nε is determined. Corollary 1 will prove that the algorithm for solving the ε-approximate solution (5) is stable with respect to changing n. Also, in Section 3, we will prove that ε-approximate solution (5) converges to the exact solution u uniformly.
According to (9), it can be rewritten by(13)Ac=b,where(14)b=ℒρk,fY+∑i=12ℬiρjγin×1,c=c1,c2,…,cnT.A=ℒρk,ℒρjY+∑i=12ℬiρkℬiρjn×n.Lemma 1.
Assumex=x1,x2,…,xnT,x=1,Ax=λx; then,(15)λ≤ℒ2+∑i=12ℬi2.Proof.
UsingA in formula (13) and assumption Ax=λx, then(16)λxk=∑j=1nℒρk,LρjY+∑i=12ℬiρkℬiρjxj=∑j=1nℒρk,xjℒρjY+∑i=12ℬiρkxjℬiρj=ℒρk,∑j=1nxjℒρjY+∑i=12ℬiρk∑j=1nxjℬiρj
Multiplyxkk=1,2,…,n on both sides of equation (16) and add all equations together, and one gets(17)λ∑k=1nxk2=∑k=1nxkℒρk,∑j=1nxjℒρjY+∑i=12∑k=1nxkℬiρk∑j=1nxjℬiρj=∑k=1nxkℒρkY2+∑i=12ℬi∑k=1nxkρk2≤ℒ2∑k=1nxk2+∑i=12ℬi2∑k=1nxk2.
Asx=1, that is, ∑k=1nxk2=1, λ=λ∑k=1nxk2≤ℒ2+∑i=12ℬi2.Lemma 2.
IfuX=1, then ℒuY≥1/ℒ−1.Proof.
Letℒu=g, according to the condition uX=1, then 1=ux=ℒ−1gx≤ℒ−1xgY, that is, ℒuY=gY≥1/ℒ−1.Corollary 1.
The method for getting theε-approximate solution (5) is stable.Proof.
In (5), ck∗ needs to be determined which is the element of the matrix c in (13). Combined with Lemma 1 and Lemma 2, it yields(18)λ≥∑k=1nxkℒρkY2=ℒ∑k=1nxkρkY2≥1ℒ−12.
So, the condition numbercondA2=λ1/λn≤ℒ2+∑i=12ℬi2/1/ℒ−12=ℒ2+∑i=12ℬi2ℒ−12.
It means that our algorithm is stable.
Note thatcondA2 denotes the condition number of matrix A and λ1 and λn represent the maximum eigenvalue and the minimum eigenvalue of A, respectively.


## 3. Constructing Orthonormal Basis for Sobolev SpaceHn0,1

We first present some notations that will be used in the following.(19)Hn0,1=ux|ux,u′x,…,un−1xis absolutelycontinuous,unx∈L20,1.The inner product is defined as follows:(20)u,vn=∑k=0n−1ukvk0+∫01unvndx,un2=u,un.By Ref. [22], Hn is the reproducing kernel space, so ux=u,Rxn, u′x=u,Rx′n, and Rx is the reproducing kernel of Hn.Specifically, for model (1) with condition (2), when the range of α is given, we will choose the appropriate space Hn0,1. Now we consider α∈1,2, so we choose space H20,1. According to the analysis in Section 2, we give the following notation:(21)Au=△1Γ2−α∫0xx−t1−αu″tdt+a1xu′x+a0xux.C1u=△u0−α0u′0,C2u=△u1+α1u′1,So, problem (1) with condition (2) can be rewritten as follows:(22)Au=f,Ciu=γi.i=1,2.In the following parts, we always assume that whenf∈L20,1, the solution of (22) exists and is unique.Now we try to give a set of multiscale orthonormal basis forH2. Considering the analysis in Section 2, we shall prove the boundedness of A and the uniform convergence of ε-approximate solution unx.By the mother wavelet,(23)φ10x=x,x∈0,12,1−x,x∈12,1,0,else.We construct the following scale functions inH010,1=u∈H10,1|u0=u1=0:(24)φikx=12i+12ix−k2i−1,x∈k2i−1,k+1/22i−1,2ik+12i−1−x,x∈k+1/22i−1,k+12i−1,0,else.i=1,2,⋯;0≤k≤2i−1−1.Theorem 3.
φ10,φ20,φ21,…,φik,…is a set of multiscale orthonormal basis forH01.Proof.
We only have to proveφ10,φ20,φ21,…,φik,… to be complete. For u∈H01, let u,φik1=0. That is, ∫01/2u′dx−∫1/21u′dx=2u1/2−u0−u1=2u1/2=0,sou1/2=0. Similarly, we can get uxik=uk/2i−1=0. Since ux is continuous, u≡0.
LetJ0ux=△∫0xutdt,J02ux=△∫0xJ0utdt.Theorem 4.
ρjxj=1∞=1,x,x2/2,J0φ10x,J0φ20x,J0φ21x,…,J0φikx,…is a set of orthonormal basis forH20,1.Proof.
It follows from Theorem3 that ρjxj=1∞ are orthonormal.
Foru∈H2, if u,12=u,x2=u,x2/22=0, then u0=u′0=u′1=0; if u,J0φik2=0, then u,J0φik2=∫01u′′φik′dx=0. In light of u∈H2, u∈H1; combining with the above results, u′∈H01. Then, u′,φik1=∫01u′′φik′dx=0. Applying Theorem 3, u′=0, so u≡0,andρjxj=1∞ is a set of complete system.Theorem 5.
The linear operatorA is bounded.Proof.
(25)AuL2=1Γ2−α∫0xx−t1−αu′′tdt+a1xu′x+a0xuxL2.
Now we mainly focus on proving the boundedness of1/Γ2−α∫0xx−t1−αu′′tdtL2. Based on the Holder¨ inequality, we obtain(26)1Γ2−α∫0xx−t1−αu″tdtL22=1Γ22−α∫01∫0xx−t1−αu″tdt2dx≤1Γ22−α∫01∫0xx−t1−αdt∫0xx−t1−αu′′2tdtdx≤M1∫01dx∫0xx−t1−αu′′2tdt=M1∫01u′′2tdt∫t1x−t1−αdx=M2∫011−t2−αu′′2tdt≤M2∫01u′′2tdt≤M2u22.
According to (25), u′x=<u,Rx′>2≤u2Rx′2, so u′2x≤Rx′22u22. It implies that a1xu′xL2≤M3u2. Similarly, a0xuxL2≤M4u2. Applying Cauchy–Schwarz inequality, one gets(27)AuL2≤M2+M3+M4u2,where Mii=1,2,3,4 is a constant, and the theorem holds true.Theorem 6.
Letu be the exact solution of (22); then, ε-approximate solution un converges to u uniformly.Proof.
Note thatunx−ux=un−u,Rx2≤un−u2Rx2, and for ∀ε>0, un−u2=A−1Aun−u2≤A−1‖Aun−u‖L2=A−1Aun−fL2<Def2.1ε⟶0, and it means un converges to u uniformly on 0,1.
By Corollary1, the algorithm we get for ε-approximate solution in H2 is stable.
Similarly, we can also construct orthonormal basis forHn to get the ε-approximate solution for (22) with higher-order fractional derivative. All the above theories hold true. For example, for H3, the orthonormal basis would be(28)1,x,x22,x36,J02φ10x,J02φ20x,J02φ21x,…,J02φikx,….
That is, our method could be widely used.


## 4. Numerical Simulations

This section shows the error estimates, convergence rates, and condition numbers of the algorithm we proposed above. The results are obtained by Mathematica 10.4; compared with Ref. [15] and Ref. [16], our method is more accurate and converges faster.Example 1.
Consider the example in Ref. [15] with α=1.9.(29)−Dαux−2x+6ux=fx,x∈0,1,u0−1α−1u′0=γ0,u1+u′1=γ1,where the exact solution ux=xα+x2α−1+1+3x−7x2+4x3+x4, f(x) = −1.82736 + 14.7159x0.1 − 4.88079x0.9 − 22.9339x1.1 − 10.9209x2.1 − 2(3 + x)(1 + 3x + x1.9 − 7x2+x2.8+4x3+x4). For this problem, we used the orthonormal basis proposed in Theorem 4 for H2 to construct the ε-approximate solution unx. The error estimates and convergence rates are all discussed in Ref. [15] and our paper. But we add the condition number condA2 to illustrate the stability discussed in Corollary 1.
Table1 shows the results of our algorithm, where n is the number of orthonormal basis, and also the number of points. The computing methods of en and rn in this paper are the same as eN and rN in Ref. [15]. Table 2 shows the results of Ref. [15], where eN=△max0≤i≤Nunxi−uxi, xi=ih,0≤i≤N,h=1/N, rN=△log2eN/e2N.
The error is much smaller in Table1 with smaller n than the results in Table 2. Obviously, Ref. [15] is first-order convergent, but the convergence rate of our method is more than second-order accuracy. Besides, we give the reference value of stability, that is, the condition number condA2.Table 1
Our method.
nenrncondA2183.98E−21349.75368.48E−32.2331379.08721.54E−32.4651385.591441.41E−42.7001387.90Table 2
Ref. [15].
NeNrN646.64E−21283.35E−20.9942561.68E−20.9985128.43E−21.000Example 2.
Consider example 2 in Ref. [16] with α=2.5.(30)Dαux+ux=6!Γ4.5x3.5+x6,x∈0,1,u0=0,u′0=0,u″0=0,with the solution ux=x6. Ref. [16] used an approximation formula, and the error reached 10−3. Here, α=2.5, so we use the orthonormal basis of H3 to construct the ε-approximate solution unx, and the result is shown in Table 3 compared with Ref. [16], where en=∑i=1nunxi−uxi2,xi=ih,i=1,…,n,h=1/n.
Table3 shows that our approximate solution converges to the exact solution quickly with not very large n, and the error of the present scheme is smaller.
In the above examples,fx is continuous. In order to illustrate the feasible range of the algorithm, we give the following example where fx is not continuous, as shown in Figure 1.Table 3
Comparison ofen.
nOur methodRef. [16]108.15E−38.73E−2151.00E−35.13E−2251.29E−42.54E−2507.82E−59.49E−3Figure 1
fx for Example 3.Example 3.
Consider the following example withα=0.5.(31)Dαux−a0xux=fx,x∈0,1,u0=0,where(32)fx=8x3/23π,x∈0,0.5,8x3/23π−1,x∈0.5,1.a0x=0,x∈0,0.5,1x2,x∈0.5,1.
The exact solution isux=x2. Here, α=0.5, so we use the orthonormal basis of H01 to construct the ε-approximate solution unx. Error en and convergence rate rn are shown in Table 4.Table 4
Results ofen and rn.
nenrn327.56E−4641.30E−42.541282.28E−52.512564.08E−62.485127.54E−72.44


## 5. Conclusion

In this paper, we focused on analyzing and obtaining theε-approximate solution for FDEs in Sobolev space Hn0,1. The sets of orthonormal basis for H01,H2,H3 were constructed and applied in the three numerical examples. It is obvious that our method is stable, converges faster, and has higher accuracy. In fact, in a similar way, we can also construct the orthonormal basis for Hnn>3, that is, we can solve higher-order FDEs.



---

*Source: 3778016-2022-06-16.xml*