A new reproducing kernel method for solving the fractional differential equations

Abstract

In this paper, we investigate the solution of fractional differential equations (FDEs) based on the reproducing kernel method. A reproducing kernel space is constructed using Legendre polynomials. Within this space, the ɛ-approximate method is applied to analyze the condition number of the matrix, and the stability and convergence of the proposed algorithm are examined. Finally, three numerical examples are presented to verify the effectiveness of the method. The obtained absolute errors are smaller than those reported in the reference literatures. This work provides an accurate and reliable numerical tool for solving FDEs.

Introduction

In recent years, there has been increasing focus on the study of fractional differential equations (FDEs). The primary reason for this is that fractional calculus operators can describe many engineering problems more accurately than traditional integer-order calculus operators.^1–3 FDEs have attracted significant attention from various scientific fields, including physics, earthquake dynamics, signal processing, fluid mechanics, ⁴ chaotic dynamics, biology, ⁵ electromagnetic waves, polymer science, thermodynamics, and so on. ⁶ FDEs are often characterized by weakly singular kernels and are more complex than integer-order differential equations. In many cases, it is difficult to obtain analytical solutions, prompting researchers to focus on approximate solutions and develop a variety of methods, ⁷ such as the finite difference method, ⁸ the local meshless method based on Laplace transforms, ⁹ the Laplace transform method, ¹⁰ variational iteration methods, ¹¹ spectral methods,^12,13 and the shooting method, ¹⁴ etc. An efficient method for solving a class of higher-order fractional differential equations with general boundary conditions is proposed by introducing appropriate basis functions. ¹⁵ Additionally, a small-parameter asymptotic method for solving fractional differential equations is proposed: by decomposing the field variable and the Riemann–Liouville integral, this method converts time-fractional differential equations into a system of two coupled integer-order partial differential equations, and solving the latter yields the solution to the original fractional differential equation. ¹⁶ To the best of our knowledge, the reproducing kernel space is widely recognized as an ideal framework for numerical analysis research. In previous work, the Taylors formula or Delta function was used to construct the reproducing kernel space,^17–20 which has been proved to be an effective tool to solving various kinds of differential equations.^21,22 Building on this foundation, the application scope of the reproducing kernel method has been further extended to solve fractional differential equations with periodic conditions in Hilbert space. ²³

In this paper, we mainly focus on the numerical solution of FDEs as follows

$\left\{\begin{array}{cc}{D}^{\alpha }u(x)+a_1(x)u^{\prime }(x)+a_0(x)u(x)=f(x),\hfill & x\in (\mathrm{0,1}),\alpha >0,\hfill \\ u(0)-{\alpha }_0{u}^{\prime }(0)={\gamma }_0,\hfill \\ u(1)+{\alpha }_1{u}^{\prime }(1)={\gamma }_1,\hfill \end{array}\right.$

(1.1)

here, D ^α u(x) denotes the Caputo fractional derivative. a _i (x) (i = 0, 1), f(x) ∈ L² [0, 1]. The existence and uniqueness of the solution u(x) of equation (1.1) are established in ²⁴ .

The remaining part of this paper is organized as follows: In the Preliminaries section, the construction of a basis for the reproducing kernel space using Legendre polynomials is presented. An efficient technique based on the ε-approximate method for solving (1.1) is then introduced, accompanied by a theoretical analysis of the uniqueness, stability, and convergence of the approximate solution for homogeneous equations. Finally, Numerical results are presented to demonstrate the validity of the proposed method.

Preliminaries

Caputo operator theory

Definition 2.1

The Caputo fractional order derivative is defined by ¹³

$D^{a}u(x)=\frac{1}{\mathrm{\Gamma }(n-a)}{\int }_0^x{(x-t)}^{n-a-1}{u}^{(n)}(t)\mathrm{d}t,a>0,$

where Γ(⋅) denotes the well-known Gamma function, and n − 1 ≤ a < n, $n\in \mathbb{N}$ .

The following properties are satisfied by the operator D ^a for n − 1 ≤ a < n

$D^a{t}^k=\frac{\mathrm{\Gamma }(k+1)}{\mathrm{\Gamma }(k+1-a)}{t}^{k-a},k\in \mathbb{N},k\ge \lceil a\rceil$

where ⌈a⌉ denotes the smallest integer greater than or equal to a. For more properties of fractional derivatives and integrals, see for example,^25–27

Basis of reproducing kernel space based on Legendre polynomials

The well-known Legendre polynomials are defined on the interval [−1, 1] and its recurrence formula

$\left\{\begin{array}{c}{L}_0(z)=1,\hfill \\ L_n(z)=\frac{1}{2^{n}n!}\frac{d^n}{d{z}^n}{(z^2-1)}^n,n=\mathrm{1,2},\cdots \hfill \end{array}\right.z\in [-\mathrm{1,1}].$

Let z = 2x − 1, we can get the following formulation in the interval [0, 1]

$\left\{\begin{array}{c}{L}_0(x)=1,\hfill \\ L_n(x)=\frac{1}{n!}\frac{d^n}{d{x}^n}{(x^2-x)}^n,n=\mathrm{1,2},\cdots \hfill \end{array}\right.x\in [\mathrm{0,1}].$

The Legendre polynomials has following properties

${\int }_0^1\sqrt{2n+1}{L}_n(x)\sqrt{2m+1}{L}_m(x)dx={\delta }_{mn}.$

The first seven terms of the polynomials of (2n + 1)^1/2L _n (x) are listed in the Table 1

Table 1.

The first seven polynomials of (2n + 1)^1/2L _n (x).

n	(2n + 1)1/2L n (x)
0	1
1	31/2 (−1 + 2x)
2	51/2 (1 − 6x + 6x2)
3	71/2 (−1 + 12x − 30x2 + 20x3)
4	3 (1 − 20x + 90x2 − 140x3 + 70x4)
5	111/2 (−1 + 30x − 210x2 + 560x3 − 630x4 + 252x5)
6	131/2 (1 − 42x + 420x2 − 1680x3 + 3150x4 − 2772x5 + 924x6)

.

Theorem 2.1

Let

${\overline{W}}_m=Span\{L_0(x),\sqrt{3}{L}_1(x),\sqrt{5}{L}_2(x),\dots ,\sqrt{2m+1}{L}_m(x)\},$

the inner product ${\overline{W}}_m$ is given

$⟨u(x),v(x)⟩={\int }_0^{1}u(x)v(x)dx,u(x),v(x)\in {\overline{W}}_m.$

Then

$R(x,y)=R_y(x)=\sum _{i=0}^m(2i+1)L_i(x)L_i(y)$

is reproducing kernel of ${\overline{W}}_m$ .

Proof. Using, ²⁸ we can prove that ${\overline{W}}_m$ is a reproducing kernel Hilbert space. Below is the proof that for $\forall u(y)\in {\overline{W}}_m$ , R _x (y) is a reproducing kernel function of ${\overline{W}}_m$ .

Let

$u(x)=\sum _{i=0}^m{a}_i\sqrt{2i+1}{L}_i(x),$

we have

$\begin{array}{ll}\hfill ⟨u(x),R_y(x)⟩& =⟨\sum _{i=0}^m{a}_i\sqrt{2i+1}{L}_i(x),\sum _{j=0}^m(2j+1)L_j(x)L_j(y)⟩\hfill \\ \hfill & =\sum _{i=0}^m{a}_i\sqrt{2i+1}⟨L_i(x),\sum _{j=0}^m(2j+1)L_j(x)L_j(y)⟩\hfill \\ \hfill & =\sum _{i=0}^m{a}_i\sqrt{2i+1}{L}_i(y)\hfill \\ \hfill & =u(y),\hfill \end{array}$

so R _x (y) is the reproducing kernel of ${\overline{W}}_m$ . □

Using ²⁹ and the reproducing kernel of ${\overline{W}}_m$ , we can get the new reproducing kernel spaces ${\stackrel{}{W}}_m$ and reproducing kernel function R _m (x, y). Next, we list some reproducing kernel spaces:

• Space $W_2=\{u(x)\mid u(x)\in {\overline{W}}_2,u(0)=0\}$ , W₂ has the same inner product as ${\overline{W}}_2$ . Its reproducing kernel

$R_2(x,y)=R(x,y)-\frac{R(x,0)R(0,y)}{R(\mathrm{0,0})}=4xy\left(12-15y+5x(-3+4y)\right).$

• • Space $W_3=\{u(x)\mid u(x)\in {\overline{W}}_3,u(0)=0,u(1)=0\}$ , W₃ has the same inner product as ${\overline{W}}_3$ . Its reproducing kernel

$R_3(x,y)={\overline{R}}_3(x,y)-\frac{{\overline{R}}_3(x,1){\overline{R}}_3(1,y)}{{\overline{R}}_3(\mathrm{1,1})}=60(-1+x)x(-1+y)y\left(4-7y+7x(-1+2y)\right),$

where ${\overline{R}}_3(x,y)=R(x,y)-R(x,0)R(0,y)/R(\mathrm{0,0})$ .

• • Space $W_4=\{u(x)\mid u(x)\in {\overline{W}}_4,u(0)=u^{\prime }(0)=u^{″}(0)=0\}$ , W₄ has the same inner product as ${\overline{W}}_4$ . Its reproducing kernel

$R_4(x,y)={\overline{\overline{R}}}_4(x,y)-\frac{\frac{{\partial }^2{\overline{\overline{R}}}_4(x,0)}{\partial y^2}\frac{{\partial }^2{\overline{\overline{R}}}_4(0,y)}{\partial x^2}}{\frac{{\partial }^4{\overline{\overline{R}}}_4(\mathrm{0,0})}{\partial x^2\partial y^2}}=8x^3{y}^3\left(56-63y+9x(-7+8y)\right),$

where ${\overline{\overline{R}}}_4(x,y)={\overline{R}}_4(x,y)-\frac{\partial {\overline{R}}_4(x,0)}{\partial y}\frac{\partial {\overline{R}}_4(0,y)}{\partial x}/\frac{{\partial }^2{\overline{R}}_4(\mathrm{0,0})}{\partial x\partial y},{\overline{R}}_4(x,y)=R(x,y)-\frac{R(x,0)R(0,y)}{R(\mathrm{0,0})}$ .

An efficient technique based on ɛ-approximate method for equation

Theoretical analysis of the algorithm for equation

Let

$v(x)=u(x)+(-\hspace{0.17em}{\alpha }_{0}u(x)-{\gamma }_0)(1-x)+({\alpha }_1{u}^{\prime }(x)-{\gamma }_1),$

then after homogenization, equation (1.1) are converted to the following form

$\left\{\begin{array}{c}{D}^{\alpha }v(x)+a_1(x)v^{\prime }(x)+a_0(x)v(x)=g(x),x\in (\mathrm{0,1}),\hfill \\ v(0)=v(1)=0,\hfill \end{array}\right.$

(3.1)

where $g(x)=f(x)+D^{\alpha }\left((-\hspace{0.17em}{\alpha }_{0}u(x)-{\gamma }_0)(1-x)+({\alpha }_1{u}^{\prime }(x)-{\gamma }_1)\right)+a_1(x){\left((-\hspace{0.17em}{\alpha }_{0}u(x)-{\gamma }_0)(1-x)+({\alpha }_1{u}^{\prime }(x)-{\gamma }_1)\right)}^{\prime }+a_0(x)\left((-\hspace{0.17em}{\alpha }_{0}u(x)-{\gamma }_0)(1-x)+({\alpha }_1{u}^{\prime }(x)-{\gamma }_1)\right)$ .

According to (3.1), a linear operator $\mathcal{P}:W_m[\mathrm{0,1}]\to L^2[\mathrm{0,1}]$ is defined by

$\mathcal{P}v(x)=D^{\alpha }v(x)+a_1(x)v^{\prime }(x)+a_0(x)v(x),\forall v(x)\in W_m[\mathrm{0,1}].$

(3.2)

It has the following properties.

Theorem 3.1

The operator $\mathcal{P}$ defined in (3.2) is a bounded and linear operator.

Proof. Obviously, $\mathcal{P}$ is linear. By applying Cauchy-Schwartzs inequality, we have

$\Vert \mathcal{P}v(x){\Vert }_{L^2}\le \Vert D^{\alpha }v(x){\Vert }_{L^2}+\Vert a_1(x)v^{\prime }(x){\Vert }_{L^2}+\Vert a_0(x)v(x){\Vert }_{L^2}.$

(3.3)

Let v(x) ∈ W _m [0, 1], by the reproducibility property of the reproducing kernel function R _m (x, y) ∈ W _m , there exist positive constants S _i (i = 0, 1, 2, ⋯ ), such that

$|v^{(i)}(x)|=|⟨v(x),\frac{{\partial }^{(i)}{R}_m(\cdot ,x)}{\partial x^{(i)}}⟩|\le \Vert \frac{{\partial }^{(i)}{R}_m(\cdot ,x)}{\partial x^{(i)}}\Vert \Vert v(x){\Vert }_{W_m}\le S_i\Vert v(x){\Vert }_{W_m}.$

By employing the aforementioned formula and utilizing Cauchy–Schwarzs inequality, a direct calculation can be performed to ascertain the existence of positive constants C and S, such that

$|D^{\alpha }v(x)|\le \frac{1}{\mathrm{\Gamma }(m-\alpha )}{\int }_0^x|{(x-t)}^{m-\alpha -1}||v{(x)}^{(m)}(t)|dt\le \frac{M\Vert v(x){\Vert }_{W_m}}{\mathrm{\Gamma }(m-\alpha )}{\int }_0^x|{(x-t)}^{m-\alpha -1}|dt\le \frac{CS\Vert v(x){\Vert }_{W_m}}{\mathrm{\Gamma }(m-\alpha )}.$

Then

$\Vert D^{\alpha }v(x){\Vert }_{L_2}^2\le {\int }_0^1\frac{C^2{S}^2{\Vert v(x)\Vert }_{W_m}^2}{{\mathrm{\Gamma }}^2(m-\alpha )}dx\le M_1^2\Vert v(x){\Vert }_{W_m}^2.$

(3.4)

Similarly, we have

$\Vert a_1(x)v^{\prime }(x){\Vert }_{L_2}\le M_2\Vert v(x){\Vert }_{W_m},\Vert a_0(x)v(x){\Vert }_{L_2}\le M_3\Vert v{\Vert }_{W_m}.$

(3.5)

Applying (3.3), (3.4), and (3.5), one gets

$\Vert \mathcal{P}v(x){\Vert }_{L^2}\le (M_1+M_2+M_3)\Vert v(x){\Vert }_{W_m},$

where M _i (i = 1, 2, 3) is a constant, so the theorem holds true. □

Let ${\{x_i\}}_{i=0}^{\infty }$ be a set of nodes in the interval [0,1], so we construct the bases ${\varphi }_i^m(x)\in W_m$ , m = 2, 3, ⋯,

${\varphi }_i^m(x)=\mathcal{P}{R}_m(y,x){|}_{y=x_i},i=\mathrm{0,1,2},\dots ,\infty .$

In the following, we solve equation (3.1) via ɛ-approximate method by combining the least square theory and the residual idea.

Definition 3.1

For ∀ɛ > 0, if $\Vert \mathcal{P}{v}_n(x)-g{\Vert }_{L^2}<\epsilon$ , then v _n (x) is called ɛ-approximate solution of (3.1).

According to (1.4), when k = m − 1 (m = 2, 3, ⋯ ), there is $v(x),{\varphi }_i^m\in W_m$ , let

$v_n(x)=\sum _{i=0}^n{c}_i{\varphi }_i^m(x),$

then v _n (x) is the ɛ-approximate solution for (3.1), and have

$G(c_0,c_1,c_2,\dots ,c_n)=\Vert \mathcal{P}{v}_n(x)-g{\Vert }_{L^2}^2=\Vert \sum _{i=0}^n{c}_i\mathcal{P}{\varphi }_i^m(x)-g{\Vert }_{L^2}^2<{\epsilon }^2,$

(3.6)

where c _i (i = 0, 1, 2, …, n) are constants.

The ɛ-approximate method for equation (3.1) is as follows: seeking $c_0^{*},c_1^{*},c_2^{*},\dots ,c_n^{*}$ such that

$G(c_0^{*},c_1^{*},c_2^{*},\dots ,c_n^{*})=\Vert \sum _{i=0}^n{c}_i^{*}\mathcal{P}{\varphi }_i^m(x)-g{\Vert }_{L^2}^2\triangleq min\Vert \sum _{i=0}^n{c}_i\mathcal{P}{\varphi }_i^m(x)-g{\Vert }_{L^2}^2.$

(3.7)

The coefficients $c_i^{*}(i=\mathrm{0,1,2},\dots ,n)$ are identified as the minimum value of G (c₀, c₁, c₂, …, c _n ). Therefore, $c_i^{*}$ satisfies the following

$\frac{\partial }{\partial c_i^{*}}G(c_0^{*},c_1^{*},c_2^{*},\dots ,c_n^{*})=0,$

which implies that

$\sum _{i=0}^n{c}_i^{*}{〈\mathcal{P}{\varphi }_i^m(x),\mathcal{P}{\varphi }_j^m(x)〉}_{L^2}={〈\mathcal{P}{\varphi }_j^m(x),g^{*}〉}_{L^2},j=\mathrm{0,1,2},\dots ,n.$

(3.8)

thus, $c_0^{*},c_1^{*},c_2^{*},\dots ,c_n^{*}$ can be determined by

$AX=d,$

where

$A={(〈\mathcal{P}{\varphi }_i^m(x),\mathcal{P}{\varphi }_j^m{(x)〉}_{L^2})}_{(n+1)\times (n+1)},X={(x_0,x_1,x_2,\dots ,x_n)}^T,d={(〈\mathcal{P}{\varphi }_j^m(x),g^{*}〉)}_{(n+1)\times 1}^T.$

Theorem 3.2

The numerical scheme (3.7) has a unique solution provided that $\mathcal{P}$ is reversible.

Proof. It is sufficient to prove the linear system AX = 0 has only zero solution. Assume that

${⟨\sum _{i=0}^n\mathcal{P}{\varphi }_i^m(x)x_i,\mathcal{P}{\varphi }_j^m(x)⟩}_{L^2}=0,j=\mathrm{0,1,2},\dots ,n.$

(3.9)

Multiplying both sides of (3.9) by x _j and a cumulative summation yields that

${⟨\sum _{i=0}^n\mathcal{P}{\varphi }_i^m(x)x_i,\sum _{j=0}^n\mathcal{P}{\varphi }_j^m(x)x_j⟩}_{L^2}=0,$

which leads to

$\mathcal{P}\left(\sum _{i=0}^n{\varphi }_i^m(x)x_i\right)=0.$

By the linear independence of ${\{{\varphi }_i^m\}}_{i=0}^n$ and reversibility of the operator $\mathcal{P}$ , we immediately obtain that x _i = 0, i = 0, 1, 2, …, n. □

Stability and convergence analysis

We will first discuss the stability of our proposed method. The condition number of matrix A is defined as follows

$cond(A)=|\frac{{\mu }_{\mathit{max}}}{{\mu }_{\mathit{min}}}|,$

where μ _max and μ _min are the maximum and minimum eigenvalues of A, respectively.

Lemma 3.1

The eigenvalue of A obtained by (3.8) is bounded by $\Vert \mathcal{P}\Vert$ .

Proof. Let μ is an eigenvalue of $A\left(A={({⟨\mathcal{P}{\varphi }_i^m(x),\mathcal{P}{\varphi }_j^m(x)⟩}_{L^2})}_{(n+1)\times (n+1)}\right)$ , then there exists a vector $X\left(X={(x_0,x_1,x_2,\dots ,x_n)}^T\right)$ , such that AX = μX. Thus

$\mu x_i=\sum _{j=0}^n{a}_{ij}{x}_j=\sum _{j=0}^n{〈\mathcal{P}{\varphi }_i^m(x),\mathcal{P}{\varphi }_j^m(x)〉}_{L^2}{x}_j={〈\mathcal{P}{\varphi }_i^m(x),\sum _{j=0}^n\mathcal{P}{\varphi }_j^m(x)x_j〉}_{L^2},\forall i=\mathrm{0,1,2},\dots ,n.$

Summing x _i from 0 to n in above formula, we derive that

$\mu =\mu \sum _{i=0}^n{x}_i^2={⟨\sum _{i=0}^n\mathcal{P}{\varphi }_i^m(x)x_i,\sum _{j=0}^n\mathcal{P}{\varphi }_j^m(x)x_j⟩}_{L^2}=\Vert \mathcal{P}\left(\sum _{i=0}^n{\varphi }_i^m(x)x_i\right){\Vert }_{L^2}^2\le \Vert \mathcal{P}{\Vert }^2\Vert \sum _{i=0}^n{\varphi }_i^m(x)x_i{\Vert }_{W_m}^2=\Vert \mathcal{P}{\Vert }^2.$

□

Lemma 3.2

If $\mathcal{P}$ is a reversible and bounded operator, then ${\mathcal{P}}^{-1}$ is bounded ³⁰

Lemma 3.3

For v ∈ W _m , if $\Vert v{\Vert }_{W_m}=1$ , then $\Vert \mathcal{P}v{\Vert }_{L^2}\ge 1/\Vert {\mathcal{P}}^{-1}\Vert$ .

Proof. Let $\mathcal{P}v=g$ , according to the condition $\Vert v{\Vert }_{W_m}=1$ , one gets

$1=\Vert v{\Vert }_{W_m}=\Vert {\mathcal{P}}^{-1}g{\Vert }_{W_m}\le \Vert {\mathcal{P}}^{-1}{\Vert }_{W_m}\Vert g{\Vert }_{L^2},$

that is

$\Vert \mathcal{P}v{\Vert }_{L^2}=\Vert g{\Vert }_{L^2}\ge \frac{1}{\Vert {\mathcal{P}}^{-1}\Vert }.$

□

Theorem 3.3

The numerical scheme for getting the ɛ-approximate solution of (3.7) is stable.

Proof. Denote μ be the eigenvalue of A obtained by (3.7), then we have

$\mu =\Vert \sum _{i=0}^n{x}_i\mathcal{P}{\varphi }_i^m(x){\Vert }_{L^2}^2.$

Combining the fact that $\Vert \sum _{i=0}^n{x}_i{\varphi }_i^m(x){\Vert }_{W_m}=1$ and Lemma 3.3, yields

$\mu =\Vert \mathcal{P}\sum _{i=0}^n{x}_i{\varphi }_i^m(x){\Vert }_{L^2}^2\ge \frac{1}{\Vert {\mathcal{P}}^{-1}{\Vert }^2}.$

So, the condition number

$cond(A)=|\frac{{\mu }_{\mathit{max}}}{{\mu }_{\mathit{min}}}|\le \frac{\Vert \mathcal{P}{\Vert }^2}{\frac{1}{\Vert {\mathcal{P}}^{-1}{\Vert }^2}}=\Vert \mathcal{P}{\Vert }^2\Vert {\mathcal{P}}^{-1}{\Vert }^2.$

It means that the condition number is bounded, so our algorithm is stable. ³¹ □

Next, we will provide the convergence analysis of the presented scheme.

Theorem 3.4

Let v ∈ W _m [0, 1] be the exact solution of (3.1), then approximate solution v _n ∈ W _m [0, 1] converges to v uniformly.

Proof. Note that

$|v(x)-v_n(x)|=|{⟨v(x)-v_n(x),R_m(\cdot ,x)⟩}_{W_m}|\le \Vert v(x)-v_n(x){\Vert }_{W_m}\Vert R_m(\cdot ,x){\Vert }_{W_m},$

and for ∀ɛ > 0, have

$\Vert v-v_n{\Vert }_{W_m}=\Vert {\mathcal{P}}^{-1}\mathcal{P}(v-v_n){\Vert }_{W_m}\le \Vert {\mathcal{P}}^{-1}\Vert \Vert \mathcal{P}(v-v_n){\Vert }_{L^2}=\Vert {\mathcal{P}}^{-1}\Vert \Vert g-\mathcal{P}{v}_n{\Vert }_{L^2}<\epsilon \to 0.$

It means v _n converges to v uniformly on interval [0,1]. □

Similarly, we can proof each order derivative $v_n^{\prime }(x),v{″}_n(x),\cdots$ of ɛ-approximate solution for (3.1) uniformly converge to v′(x), v′′(x), ⋯ respectively.

Numerical results

In this section, three numerical examples are presented to demonstrate the effectiveness of the proposed algorithm.

Example 4.1

When α = 0.5

$D^{\alpha }u(x)-a_0(x)u(x)=f(x),x\in [\mathrm{0,1}],u(0)=0,$

where

$f(x)=\left\{\begin{array}{c}\frac{8x^{\frac{3}{2}}}{3\sqrt{\pi }},x\in [\mathrm{0,0.5}],\hfill \\ \frac{8x^{\frac{3}{2}}}{3\sqrt{\pi }}-1,x\in (0.5,1].\hfill \end{array}\right.a_0(x)=\left\{\begin{array}{c}0,x\in [\mathrm{0,0.5}],\hfill \\ \frac{1}{x^2},x\in (0.5,1].\hfill \end{array}\right.$

The exact solution is u(x) = x². We use the basis of W₂ to construct the ɛ-approximate solution u _n (x). Error $e^n\triangleq \sqrt{\sum _{i=0}^n\left(u(x_i)-u_n(x_i)\right)}$ (x _i = ih, i = 0, 1, 2, …, n, h = 1/n). The results are presented in Table 2, compared with Ref. 32.

Table 2.

Results of e ⁿ compared with Ref. 32 for Example 1.

n	e n	e n (ref. [32])
32	1.75 × 10−12	7.56 × 10−4
64	5.76 × 10−12	1.30 × 10−4
128	5.57 × 10−12	2.28 × 10−5
256	2.56 × 10−12	4.08 × 10−6
512	2.70 × 10−11	7.54 × 10−7

Example 4.2

Consider the example 2 in Ref. 32 with α = 1.9

$\left\{\begin{array}{c}{D}^{\alpha }u(x)-(2x+6)u(x)=f(x),x\in (\mathrm{0,1}),\hfill \\ u(0)-\frac{1}{\alpha -1}{u}^{\prime }(0)=\frac{\alpha -4}{\alpha -1},u(1)+u^{\prime }(1)=3\alpha +8,\hfill \end{array}\right.$

where f(x) = −1.82736 + 14.7159x^0.1 − 4.88079x^0.9 − 22.9339x^1.1 − 10.9209x^2.1 − 2 (3 + x) (1 + 3x + x^1.9 − 7x² + x^2.8 + 4x³ + x⁴), the exact solution u(x) = x ^α + x^2α−1 + 1 + 3x − 7x² + 4x³ + x⁴. For α = 1.9, we use the basis of W₃ to get the ɛ-approximate solution u _n (x), the error estimates are all discussed in Ref. 32 and our method. Let $e^N\triangleq \underset{1\le i\le N}{\mathit{max}}|\left(u(x_i)-u_N(x_i)\right)|(x_i=i/N,i=\mathrm{0,1,2},\dots ,N)$ . Table 3 presents a comparison of our results with those reported in Ref. 32. Furthermore, the comparison between the methods proposed in this paper and the results from Ref. 33 and 34 is presented in Table 4. As can be seen from Tables 3 and 4, when n takes different values, the approximate solutions converge rapidly to the exact solutions with small errors.

Table 3.

Results of e ^N compared with Ref. 32 for Example 2.

N	e N	e N (ref. [32])
18	5.09 × 10−4	3.98 × 10−2
36	5.08 × 10−4	8.48 × 10−3
72	5.09 × 10−4	1.54 × 10−3
144	5.09 × 10−4	1.41 × 10−4

Table 4.

Results of e ^N compared with Ref. 33 and 34 for Example 2.

N	e N	e N (ref. [33])	e N (ref. [34])
64	7.14 × 10−4	1.13 × 10−1	6.64 × 10−2
128	6.85 × 10−4	5.88 × 10−2	3.35 × 10−2
256	8.57 × 10−4	3.03 × 10−2	1.68 × 10−2
512	7.38 × 10−4	1.56 × 10−2	8.43 × 10−3

Example 4.3

Consider the following example with 2 ≤ α ≤ 3,

$\left\{\begin{array}{c}{D}^{\alpha }u(x)+u(x)=\frac{6!}{\mathrm{\Gamma }(4.5)x^{3.5}}+x^6,x\in (\mathrm{0,1}),\hfill \\ u(0)=0,u^{\prime }(0)=0,u^{″}(0)=0,\hfill \end{array}\right.$

with the solution u(x) = x³.

In the Table 5, we give errors of each derivative of different n when α = 2.3, where

Table 5.

Results of errors of each derivative for Example 3.

n	$e_0^n$	$e_1^n$	$e_2^n$	CPU(s)
10	2.85 × 10−16	2.27 × 10−15	1.62 × 10−14	13.734
20	9.81 × 10−16	1.43 × 10−15	4.92 × 10−15	16.063
40	7.22 × 10−15	3.91 × 10−14	4.69 × 10−14	25.094

$e_k^n\triangleq \sqrt{\sum _{i=0}^n\left(u^{(k)}(x_i)-u_n^{(k)}(x_i)\right)},k=\mathrm{0,1,2}$ . When the error precision meets the requirement of 10⁻¹⁴, the time consumed is also recorded in the Table 5.

Figures 1 and 2 show that the absolute-errors corresponds to different the number of basis n and different values of fractional parameter α.

Figure 1.

The Absolute-errors with α = 2.3 for different the number of basis n.

Figure 2.

The Absolute-errors with n = 10 for different values of fractional parameter α.

Conclusion

This paper proposes an efficient numerical scheme for solving FDEs by integrating the reproducing kernel space theory with the ɛ-approximate method. The method proposed in this paper has the following advantages. First, in the absence of an exact solution, the approximate solution obtained by this method exhibits good agreement with the theoretical exact solution. Both the approximate solution and its derivatives converge uniformly to the exact solution. The numerical examples further verify the effectiveness and superiority of the proposed method. Second, the results demonstrate that the errors of the scheme are smaller than those reported in the reference literature. Furthermore, the computational time of this method is relatively short. Third, this method boasts strong applicability in physics, including transonic multiphase flows and the electric circuit model of fractional stiff differential equations.^35,36 Moreover, it is extensible to solving fractional integro-differential equations, stochastic fractional differential equations and high-dimensional systems.