Galerkin finite element method for nonlinear fractional differential equations

Abstract

In this paper, we study the existence, regularity, and approximation of the solution for a class of nonlinear fractional differential equations. In order to do this, suitable variational formulations are defined for nonlinear boundary value problems with Riemann-Liouville and Caputo fractional derivatives together with the homogeneous Dirichlet condition. We investigate the well-posedness and also the regularity of the corresponding weak solutions. Then, we develop a Galerkin finite element approach for the numerical approximation of the weak formulations and drive a priori error estimates and prove the stability of the schemes. Finally, some numerical experiments are provided to demonstrate the accuracy of the proposed method.

1 Introduction

In the current study, we consider a fractional order nonlinear boundary value problem as follows: Find u such that

$$ -_{0}^{\circ}{D^{s}_{x}} u(x) +g(x,u(x)) =f(x),\quad x\in {\Omega}:=[0,1], $$

(1)

$$ u(0)=0,\quad u(1)=0, $$

where s ∈ (1,2), and $_{0}^{\circ }{D^{s}_{x}}$ refers to either Riemann-Liouville or Caputo fractional derivatives which are detailed below. Furthermore, f and g are known functions chosen from the suitable function spaces defined in Section 2.

Developing the order of differentiation to any real number is an interesting question. To find an affirmative answer, some efforts have been done and different types of so-called fractional derivatives have been introduced [20]. It is notified that the fractional derivative is a concept with attractive applications in science and engineering. It appears in the anisotropic diffusion modeling anomalously for cardiac tissue in microscopic and macroscopic levels. Furthermore, fractional derivative models are certain instances of nonlocal models which are introduced in comparison with the classical ones [12].

Similar to the ordinary and partial differential equations, one follows two approaches in seeking solutions for fractional differential equations (FDEs): analytic and numeric solutions. The analytical methods such as the Fourier, Laplace and Mellin transform methods, and even Green function approach (fundamental solution) are available for some special types of FDEs [20, 30]. In practice, we have to apply numerical methods due to the lack of applicability of analytical methods for a wide range of FDEs. Hence, the study of the numerical approaches for them is of great importance.

It is worthy to mention that in spite of different definitions for fractional differential operators, most of them are defined by the Abel integral operator. Among the popular numerical approaches to Abel’s integral equation, one could mention the collocation method [7] and the Galerkin method based on piecewise polynomials [13, 32] where the different varieties of those approaches, in general perspective, spectral and projection methods, could be utilized to find approximation for FDEs [24, 30]. Converting FDEs to a suitable integral equation, and then solving it numerically with the abovementioned approaches [33] or directly solving it by finite difference method [25], are based on the adequate regularity assumptions of the strong solution which is not available in general [14, 18]. Here, we contribute to the study on the numerical solution of one-dimensional nonlinear FDEs which involve Riemann-Liouville or Caputo derivative by introducing an appropriate weak formulation and describing the Galerkin solution with convergence analysis in some appropriate functional spaces.

The fractional operator in (1) is nonlocal and g is a nonlinear function with respect to u, so the study of the existence, uniqueness, and regularity of solution and the numerical investigation are challenging. The existence of classical solution for the nonlinear FDEs is considered in [38]; and for the linear operator in one dimension, the regularity of the solution is investigated in [14]. In this work, we explore the issues of the existence and uniqueness with the aid of Browder-Minty method of the monotone operator theory.

In the recent literature, due to their applications in science and engineering [24], several types of numerical methods have been proposed for the approximation of FDEs. The theory and the numerical solution of a linear Riemann-Liouville and Caputo FDEs with two-point boundary condition have been extensively studied in [21, 27]. In those works, the fractional boundary value problem is reformulated appropriately in terms of the Volterra integral equation, and then the numerical approach is proceed by some suitable schemes such as the piecewise polynomial collocation and the spectral Galerkin methods. Spectral and pseudo-spectral methods are some of the interesting numerical approaches which are taken into consideration for the FDEs. Among the whole research on this area, we can mention [26, 35,36,37] wherein the Jacobi polynomials play a crucial role in the construction of the approximation. The attention to this class of orthogonal polynomials is a motivation for introducing a generalization of them with application in the numerical solution of FDEs [9].

In this work, we study the Galerkin finite element method for Riemann-Liouville and Caputo nonlinear fractional boundary value problems of Dirichlet type. The finite element method is a popular numerical approach in order to find an approximation for nonlinear differential equations [17] and [34]. The finite element solutions of quasi-linear elliptic problems with non-monotone operators have been considered in [1] and [16]. Furthermore, in [15] under the assumptions of the strong monotonicity and Lipschitz continuity of the corresponding second order nonlinear elliptic operator, a linear order approximation by finite element method has been obtained. In this paper, the investigated fractional nonlinear operators have the monotonicity and Lipschitz continuity properties, which are crucial for our analysis.

The reasonable energy space associated with the non-local operators is fractional Sobolev space. Also, due to the presence of the nonlinear term in (1), we utilize Musielak-Orlicz space in order to introduce a suitable functional space by intersection of the two mentioned spaces in a convenient way. Then, with the aid of the monotone operator theory, the coercivity of the nonlinear variational formulation along with the Riemann-Liouville and Caputo fractional derivatives is investigated. This approach leads to getting a unique weak solution which can be approximated by the finite element method. Two main features of this research are as follow:

• We study the existence and uniqueness issue of the weak solution for the main problem (1) utilizing the monotone operator theory in some Musielak-Orlicz and fractional Sobolev spaces.

• We develop the Galerkin finite element method for nonlinear FDEs and obtain a priori error estimates for the method based on the generalization of the Céa’s lemma.

We organize the reminder of the paper as follows: in Section 2, some introduction regarding the fractional calculus, the semi-linear monotone operator, and also the suitable functional spaces is briefly presented. Section 3 is devoted to the variational formulation of nonlinear boundary value problems along with Riemann-Liouville and Caputo derivatives. Furthermore, the regularity of the solution is studied in this section. The numerical approximation of the weak solution is examined by finite element approximation in Section 4 along with the full study of the existence and uniqueness issue of the discrete equations and the convergence and stability of the method. In numerical experiments section, some FDEs are solved by the finite element method. Finally, we provide some conclusions and further remarks for the future works.

2 Preliminaries

This section is devoted to some preliminaries to fractional calculus to provide an introduction to the problem considered in the paper. The energy space regarding fractional operators, fractional Sobolev spaces are introduced in this section. Then, in order to ensure the existence of the weak solution by monotonicity arguments, some preface to nonlinear functions on Musielak-Orlicz spaces is provided.

2.1 Fractional calculus

To make this paper self-contained, we recall the Riemann-Liouville and the Caputo fractional integral and derivatives from [20]. For any s > 0 with n − 1 < s < n, $n \in \mathbb {N}$, the right- and left-sided fractional integrals on the bounded interval [a,b] are as follows:

Left fractional integral operator is defined as

$$ ({~}_{a}{I^{s}_{x}}u)(x)= \frac{1}{\Gamma(s)}{{\int}_{a}^{x}}(x-y)^{s-1} u(y)\mathrm{d}y, $$

(2)

while the right fractional integral operator is given by

$$ ({~}_{x}{I^{s}_{b}}u)(x)= \frac{1}{\Gamma(s)}{{\int}_{x}^{b}}(y-x)^{s-1} u(y)\mathrm{d}y. $$

(3)

Left-sided Riemann-Liouville fractional derivative of the order s for the function u ∈ Hⁿ(Ω) can be defined as

$$ {~}^{R}_{a}{D_{x}^{s}} u=D^{n} _{0}I_{x}^{n -s}u, $$

(4)

where the operator Dⁿ denotes the classical derivative of the order n. The corresponding right-sided Riemann-Liouville fractional derivative is stated as

$$ {~}^{R}_{x}{D_{b}^{s}} u=(-1)^{n} D^{n} _{x}I_{b}^{n-s}u. $$

(5)

In addition, the left-sided Caputo derivative of the order s is given by

$$ {~}^{C}_{0}{D_{x}^{s}} u= _{0}{I_{x}^{s}} D^{n} u, $$

(6)

where the following relation defines the right-sided Caputo derivative

$$ {~}^{C}_{x}{D_{b}^{s}} u={(-1)^{n}}_{x}{I_{b}^{s}}D^{n} u. $$

(7)

From the above definitions, it is apparent that the Abel integral operator has a significant role in defining fractional derivatives.

2.2 Some properties of semi-linear operators

Let V be a real Banach space and V^∗ as its dual space. Let define the duality pairing

$$ V^{*}\times V \ni (y,x) \rightarrow \langle y,x \rangle =y(x), $$

which means that 〈y,x〉 is the value of a continuous linear functional y ∈ V^∗ on an element x ∈ V and ∥.∥ and ∥.∥_∗ are the norms associated with V and V^∗, respectively. Due to the central role of the monotonicity property, we present a formal definition for this concept.

Definition 1

([4, 22]) Let V be a separable Banach space. An operator F is called monotone on V if

$$ \langle Fx-Fy,x-y\rangle \geq 0, \quad \forall x,y\in V, $$

it is called strictly monotone if

$$ \langle Fx-Fy,x-y\rangle > 0, \quad \forall x,y\in V, \quad x\neq y, $$

it is called coercive

$$ \langle Fx,x \rangle > \gamma (\Vert x \Vert). \Vert x \Vert, \quad \text{where} \quad \gamma(s)\rightarrow \infty \quad \text{as}\quad s\rightarrow \infty, $$

it is called hemi-continuous if the real-valued function

$$ s\rightarrow \langle F(u+s.v), w\rangle, $$

is continuous on [0,1] for every fixed u,v and w ∈ V. Furthermore, in the terminology of the article [19], the operator $F:V\rightarrow V^{*}$ with the domain D = D(F) ⊂ V is hemi-continuous if for u ∈ D and w ∈ V, we get $F(u+t_{n}w)\rightarrow F(u),$ when the sequence t_n tends to 0.

In the following, we state the Browder-Minty theorem which is utilized to prove the existence and uniqueness of the weak solution.

Theorem 1

(Browder-Minty) Let V be a real reflexive Banach space and a hemi-continuous monotone operator $F:V\rightarrow V^{*}$ be coercive. Then, for any g ∈ V^∗, there exists a solution u^∗∈ V of the equation

$$ F(u)=g. $$

This solution is unique if F is a strictly monotone operator.

Proof

See the details of the proof in either [4] or [22, Chapter 19]. □

2.3 Functional spaces

To formulate an appropriate functional space so that the problem (1) is well-posed, we first recall the definition of fractional-order Sobolev spaces. As usual, the standard Lebesgue spaces are denoted by $L^{p}\left ({\Omega }\right )$, and their norms by $\left \Vert \cdot \right \Vert _{L^{p}\left ({\Omega }\right ) }$. For p = 2, the scalar product is denoted by $\left (u,v\right ) ={\int \limits }_{\Omega }u(x){v(x)}\mathrm {d}x$, and the norm by $\left \Vert \cdot \right \Vert =\left (\cdot ,\cdot \right )^{1/2}$. Let $\lbrace \lambda _{n}\rbrace _{n\in \mathbb {N}}$ be the set of all eigenvalues of the following boundary-value problem:

$$ \begin{array}{@{}rcl@{}} \begin{array}{ll} D^{2}u(x)&=-\lambda u(x), \quad x\in {\Omega} ,\\ \frac{du}{dt}(0)&=u(1)=0 , \end{array} \end{array} $$

(8)

where ϕ_n is an eigenfunction related to λ_n for $n \in \mathbb {N}$. Now, for $s \in \mathbb {R},$ a Hilbert scale H^s(Ω) is defined based on $\lbrace \phi _{n}\rbrace _{n\in \mathbb {N}}$ with the following scalar products and norms:

$$ (u,v)_{H^{s}({\Omega})}=\sum\limits_{n=1}^{\infty}\lambda^{s} (u,\phi_{n})(v,\phi_{n}),\quad u,v\in \text{span} \{\phi_{n}\}_{n\in \mathbb{N}}, $$

(9)

and

$$ \Vert u\Vert_{H^{s}({\Omega})}=(u,u)^{\frac{1}{2}}_{H^{s}({\Omega})}. $$

Let $\lambda _{n}:={\mu _{n}^{2}}$, then $\lbrace \mu _{n}^{-s}\phi _{n}\rbrace _{n\in \mathbb {N}}$ forms an orthonormal basis for H^s(Ω). It is well-known that $\{ \phi _{n} \}_{n \in \mathbb {N}}$ is an orthonormal basis for H⁰(Ω) = L²(Ω), so for u_n = (u,ϕ_n), we have $u={\sum }_{n=1}^{\ \infty } u_{n} \phi _{n}$. For any s ≥ 0, the fractional order Sobolev space is defined by the spectral properties of the operator (8) and the inner product (9) as follows:

$$ H^{s}({\Omega}):=\big\{ u \in L^{2}({\Omega}) \mid \sum\limits_{k=1}^{\ \infty} {\lambda^{s}_{k}} {u^{2}_{k}} < \infty \big \}, $$

(10)

for more details, see [3]. Another approach to define the fractional Sobolev space is using the definition of L^p(Ω) spaces along with the Slobodeckij semi-norm [10]. To our aim, it suffices to set p = 2 and let $\left \lfloor s\right \rfloor$ denote the largest integer for which $\left \lfloor s\right \rfloor \leqslant s$, and define $\lambda \in \left [ 0,1\right [$, by $s=\left \lfloor s \right \rfloor +\lambda$. For $s\in \mathbb {R}_{>0}\backslash \mathbb {N}$, we introduce the scalar product:

$$ \begin{array}{@{}rcl@{}} \left( \varphi,\psi\right)_{H^{s}({\Omega})} & :=&\sum\limits_{\alpha \leqslant\left\lfloor s\right\rfloor }{\left( D^{\alpha}\varphi,D^{\alpha }\psi\right) }\\ &&+{{\int}_{\Omega}{\int}_{\Omega}{\frac{\left( D^{\left\lfloor s\right\rfloor}\varphi(x)-D^{\left\lfloor s\right\rfloor}\varphi(y)\right) {\left( D^{\left\lfloor s\right\rfloor}\psi(x)-D^{\left\lfloor s\right\rfloor}\psi(y)\right) }} {|x-y|^{1+2\lambda}}\mathrm{d}x\mathrm{d}y}}, \end{array} $$

(11)

and the norm $\Vert \varphi \Vert _{H^{s}({\Omega })}:=\left (\varphi ,\varphi \right )_{H^{s}({\Omega })}^{1/2}$. For $s\in \mathbb {N}$, obviously the second term in (11) is ignored. Then, the Sobolev space $H^{s }\left ({\Omega }\right )$ is given by

$$ H^{s}\left( {\Omega}\right) :=\left\{ u\in L^{2}({\Omega})\mid\forall~0\leq k\leq\left\lfloor s\right\rfloor \quad u^{\left( k\right) }\in L^{2}\left( {\Omega}\right) \quad\text{and}\quad \Vert u\Vert_{H^{s} ({\Omega})}<\infty\right\} . $$

The dual space of H^s(Ω) is denoted by H^∗(Ω) := H^−s(Ω) and is equipped with the norm

$$ \Vert u\Vert_{H^{-s}({\Omega})}:=\sup\limits_{v\in H^{s}({\Omega})} \frac{(u,v)}{\Vert v\Vert_{H^{s}({\Omega})}}, $$

(12)

where (⋅,⋅) denotes the continuous extension of the L²-scalar product to the duality pairing 〈⋅,⋅〉 in H^−s(Ω) × H^s(Ω). Let $\tilde {H}^{s}({\Omega })$ be the set of functions in H^s(Ω) which are extended by zero to the whole domain $\mathbb {R}$. This space could also be defined by the intermediate space of the order s ∈ (0,1) given as:

$$ \begin{array}{@{}rcl@{}} {H_{0}^{s}}({\Omega})=[{H_{0}^{m}}({\Omega}), H^{0}({\Omega})]_{\theta}, \quad m\in \mathbb{Z}, \end{array} $$

(13)

where m(1 − 𝜃) = s and ${H^{m}_{0}}({\Omega })$ denotes the closure of $\mathcal {D}({\Omega })$ in H^m(Ω) [28]. In the latter, $\mathcal {D}({\Omega })$ stands for the set of all infinitely differentiable functions with compact support which is equipped with the locally convex topology. Indeed, for $\phi \in {H}_{0}^{s}({\Omega })$, ϕ and its derivatives of the order k ≤ m have the compact support property. For m = 1, we set $\tilde {H}^{s}({\Omega }):= H_{0}^{1-s}({\Omega })$. Let I_L (respectively I_R) denote half interval $(-\infty , b)$ (res. $(a, \infty )$), then $\tilde {H}^{s}_{L}({\Omega })$ (res. $\tilde {H}^{s}_{R}({\Omega })$) stands for the set of all u ∈ H^s(Ω) whose extension by zero denoted by $\tilde { u}$ are in H^s(I_L) (res. H^s(I_R)) [2, 18].

Theorem 2

([18]) Assume that n − 1 < s < n, for $n\in \mathbb {N}$. The operators $^{R}_{0}{D_{x}^{s}} u$ and $^{R}_{x}{D_{1}^{s}} u$ for $u\in \mathcal {D}({\Omega })$ can be extended continuously to operators with the same notations from $\tilde {H}^{s}({\Omega })$ to L²(Ω), i.e.,

$$ {\Vert^{R}_{0}}{D_{x}^{s}} u \Vert_{L^{2}(\mathbb{R})}\leq c\Vert u\Vert_{\tilde{H}^{s}({\Omega})}, $$

(14)

and

$$ {\Vert^{R}_{x}}{D_{1}^{s}} u \Vert_{L^{2}(\mathbb{R})}\leq c\Vert u\Vert_{\tilde{H}^{s}({\Omega})}. $$

(15)

The following theorem represents some profitable characteristics of the fractional differential and integral operators.

Theorem 3

The following statements hold:

• a) The integral operators $_{0}{I_{x}^{s}}$ and $_{x}{I_{1}^{s}}$ satisfy the semi-group property.

• b) For ϕ,ψ ∈ L²(Ω), $(_{0}{I_{x}^{s}}\phi , \psi )=(\phi , _{x}{I_{1}^{s}}\psi )$.

• c) For any s > 0, the function x^s ∈ H^α(Ω), where $0\leq \alpha <s+\frac {1}{2}.$

• d) For any non-negative α,γ, the Riemann-Liouville integral operator I^α is a bounded map from $\tilde {H}^{\gamma }({\Omega })$ into $\tilde {H}^{\gamma +\alpha }({\Omega }).$

• e) The operators $_{0}^{R}{D_{x}^{s}}:\tilde {H}^{s}_{L}({\Omega })\rightarrow L^{2}({\Omega })$ and $_{x}^{R}{D_{1}^{s}}:\tilde {H}^{s}_{R}({\Omega })\rightarrow L^{2}({\Omega })$ are continuous.

• f) For any s ∈ (0,1) and $u\in \tilde {H}^{1}_{R}({\Omega })$, $_{0}^{R}{D_{x}^{s}}u=_{0}^{R}I_{x}^{1-s} u^{\prime }$. Meanwhile, $u\in \tilde {H}^{1}_{L}({\Omega })$, then $_{x}^{R}{D_{1}^{s}}u=-_{x}^{R}I_{1}^{1-s} u^{\prime }$.

Proof

The first item has been investigated in [20, Theorem 2.4]. The Fubini theorem proves item b [20, Lemma 2.7]. The proof of other items can be found in [18]. □

2.3.1 Nonlinear functions on Orlicz spaces

Throughout this paper, we require some important properties for the nonlinear part of (1). A suitable functional space to deal with the monotone operators with nonlinear terms is the Orlicz spaces or the generalized Orlicz space which is called Musielak-Orlicz space [2, 5]. We recall some necessary definitions and properties related to the mentioned spaces.

Definition 2

A function $A: \mathbb {R}\rightarrow \left [0,\infty \right ]$ is termed an N-function if

• a) It is even and convex;

• b) A(t) = 0 if and only if t = 0;

• c) $\lim _{t\rightarrow 0} \frac {A(t)}{t} = 0$, $\lim _{t\rightarrow \infty } \frac {A(t)}{t} =\infty$.

For more details on N-function, one can see [2, Chapter VIII], [31, Chapter I], and [6, Appendix].

Definition 3

([29, Chapter II]) Let (Ω,Σ,μ) be a measure space such that μ is σ-finite and complete. We say that a function $\varphi : {\Omega } \times \mathbb {R} \rightarrow \left [ 0, \infty \right ]$ is a Musielak-Orlicz function on Ω if

• a) φ(.,t) is measurable for all $t\in \left [ 0, \infty \right ]$;

• b) For a.e. x ∈Ω,φ(x,.) is non-trivial, even and convex;

• c) φ(x,.) is vanishing and continuous at 0 for a.e. x ∈Ω;

• d) φ(x,t) > 0, ∀t > 0 and $\lim _{t \rightarrow \infty } \varphi (x,t) = \infty$.

In this paper, Σ is the σ-algebra of subsets of Ω and μ denotes the Lebesgue measure on Σ. The complementary Musielak-Orlicz function $\tilde {\varphi }$ is defined by

$$ \tilde{\varphi}(x,t) : = \sup\big\{ s\vert t\vert - \varphi(x,s)~ |~ s>0 \big\}, $$

which is the same as the complementary of the Young function (a more general class of N-functions) defined in [6, 31].

Assumption A ([3]) Assume that a nonlinear function $g(x,t):{\Omega } \times \mathbb {R} \rightarrow \mathbb {R}$ satisfies the following properties:

$$ \left\{ \! \begin{aligned} &g(x,.) \text{continuous, odd, strictly monotone,} \quad & \text{a.e. on } {\Omega},\\ &g(x,0)=0, \quad \lim\limits_{t\rightarrow \infty} g(x,t)=\infty, & \text{a.e. on } {\Omega},\\ &g(.,t) \text{is measurable}, & \forall t \in \mathbb{R}. \end{aligned} \right. $$

Note that the inverse of the function g(x,.) exists due to the strictly monotone property. Let us denote it by $\tilde {g}(x,.).$ We define G(x,t) and $\tilde {G}(x,t)$ by

$$ G(x,t):={\int}_{0}^{\vert t\vert }g(x,s)\mathrm{d}s,\quad \quad \tilde{G}(x,t):={\int}_{0}^{\vert t\vert }\tilde{g}(x,s)\mathrm{d}s. $$

These functions are complementary Musielak-Orlicz functions that are N-functions with respect to the second variable.

Definition 4

Let G(x,.) be an N-function. We say this function satisfies the global (Δ₂)-condition if there exists a constant c ∈ (0,1] such that for a.e. x ∈Ω and for all t ≥ 0

$$ c t g(x,t)\leq G(x,t)\leq t g(x,t), $$

where the function g(x,t) satisfies in Assumption A.

Let ${\mathscr{M}}({\Omega })$ represent all real-valued measurable functions defined on Ω, then the Musielak-Orlicz space generated by G is defined as follows:

$$ L_{G}({\Omega}):=\Big\lbrace u \in \mathcal{M}({\Omega}) \mid G(.,u(.))\in L^{1}({\Omega})\Big\rbrace, $$

which means that the modular

$$ \begin{array}{ll} \rho_{G}&: \mathcal{M}({\Omega}) \rightarrow \left[0, \infty\right],\\ \rho_{G} (u)&:= {\int}_{\Omega}G(t,u(t))\mathrm{d}t, \end{array} $$

is measurable [5]. If G(x,.) and $\tilde {G}(x,.)$ satisfy Definition 4, then by Theorem 8.9 in [2], it is concluded that L_G(Ω) equipped with the Luxemburg norm;

$$ \Vert u\Vert_{G,{\Omega}}:=\inf\Big\lbrace m>0 \mid \rho_{G}(\frac{u}{m})\leq 1\Big\rbrace, $$

is a reflexive Banach space. The same result is valid for the space $L_{\tilde {G}}({\Omega })$ with the norm $\Vert .\Vert _{\tilde {G},{\Omega }}$. Moreover, in our analysis we need the generalized Hölder inequality given by

$$ \Big\vert {\int}_{\Omega} u(t)v(t) \mathrm{d}t \Big\vert \leq 2\Vert u\Vert_{G,{\Omega}} \Vert v \Vert_{\tilde{G},{\Omega}}, \quad \forall u \in L_{G}({\Omega}), ~~ \forall v\in L_{\tilde{G}}({\Omega}). $$

(16)

There are some proofs for this relation that one can see them in [23, Chapter 3] and [2]. Furthermore, the following important result:

$$ \lim\limits_{\Vert u\Vert_{G, {\Omega}}\rightarrow \infty} \frac{\rho_{G}(u)}{\Vert u\Vert_{G, {\Omega}}}=\infty, $$

(17)

which is obtained in [6] has a significant role in the applicability of the monotone operator theorems for our target.

Lemma 1

([6]) If g(.,u(.)) satisfies the (Δ₂)-condition, then for all u ∈ L_G(Ω) one can get $g(.,u(.))\in L_{\tilde {G}}({\Omega })$.

Now, we are ready to define the suitable function space which is appropriate for our problem.

Definition 5

Let 1 < s < 2. Under the Assumption A which has been fulfilled in Definition 4, consider the following reflective Banach space U as

$$ U:=U({\Omega}, G):=\lbrace \phi \in \tilde{H}^{\frac{s}{2}}({\Omega}) \mid G(.,\phi)\in L^{1}({\Omega}) \rbrace, $$

where equipped with the norm

$$ \Vert u\Vert_{U}:=\Vert u\Vert_{\tilde{H}^{\frac{s}{2}}({\Omega})}+\Vert u\Vert_{{G,{\Omega}}}. $$

(18)

We state the following lemma from [3] which is important in the investigation of the regularity of the solution.

Lemma 2

Let 0 ≤ s < 1 and assume that for M > 0, there exists a constant l_M such that g satisfies

$$ \vert g(x,u_{1}) -g(x,u_{2}) \vert\leq l_{M} \vert u_{1}- u_{2} \vert, \quad x,y \in {\Omega},\quad u_{i}\in \mathbb{R} \ \text{with} \ \vert u_{i} \vert \leq M. $$

(19)

Then, for $u \in H^{s}({\Omega })\cap L^{\infty }({\Omega })$, we have g(.,u(.)) ∈ H^s(Ω).

3 Variational formulation and regularity

In this section, we aim to work with an appropriate variational formulation to overcome the difficulty of dealing with the nonlinear and fractional terms of the main problem for both cases of the Riemann-Liouville and the Caputo derivatives separately. The nonlocal variational problems possess reduced order smoothing properties which are investigated in this section.

3.1 The Riemann-Liouville fractional operator

The appropriate variational formulation of the problem (1) in the linear case (with g(x,u) = 0) introduced in [18] is:

Find $u\in U=\tilde {H}^{\frac {s}{2}}({\Omega })$ such that

$$ \mathcal{A}(u,v)=(f,v), \quad v\in V=U, $$

(20)

where

$$ \mathcal{A}(u,v):=-(_{0}^{R}D^{\frac{s}{2}}_{x} u , _{x}^{R}D^{\frac{s}{2}}_{1} v), $$

(21)

and f ∈ L²(Ω).

Considering the above form for the fractional part has some pros in utilizing the nice properties indicated in Theorem 3 for the approximation procedure. Hence, for the nonlinear problem (1), the weak formulation is stated as follows:

Find u ∈ U satisfying

$$ \begin{array}{ll} \mathcal{L}(u,v)&:= \mathcal{A}(u,v)+ \mathcal{B}(u,v)=\langle f,v\rangle:=F(v),\quad \quad v \in U, \end{array} $$

(22)

where f ∈ U^∗, ${\mathscr{B}}(u,v):=(g(x,u),v)$ and U^∗ is the dual space of the reflexive Banach space U introduced in Definition 5 and their duality map is denoted by 〈.,.〉.

Theorem 4

Let 1 < s < 2 and $u\in \tilde {H}^{\frac {s}{2}}({\Omega }),$ the operator $\mathcal {A}(u,v)$ is coercive and monotone, i.e.,

$$ \exists ~c>0 \quad \text{s.t}\quad \mathcal{A}(u,u)\geq c\Vert u \Vert^{2}_{\tilde{H}^{\frac{s}{2}}({\Omega})}. $$

Proof

It is easily verified that

$$ \mathcal{A}(u,u)=-(_{0}^{R}{D^{s}_{x}}u,u). $$

Let us consider $Su(x):= -_{0}^{R}{D^{s}_{x}}u$ and borrow the notation S^εu for ε > 0 from [13] as follows:

$$ S^{\varepsilon} u(x)=\frac{-1}{\Gamma(2-s)}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}{{\int}_{0}^{x}}(x-t)^{1-s}e^{-\varepsilon(x-t)}u(t)\mathrm{d}t, \quad x>0. $$

Using the Plancherel theorem, we get:

$$ (S^{\varepsilon} u,u)={\int}_{-\infty}^{\infty}(S^{\varepsilon} u\hat{)}(w) \overline{\hat{u}(w)}\mathrm{d}w, $$

(23)

where the notation $\hat {}$ refers to the Fourier transform. Let us introduce

$$ \hat{a}_{\varepsilon}(w)= w^{2} {\int}_{0}^{\infty}x^{1-s} e^{-(\varepsilon+\mathrm{i} w)x}\mathrm{d}x, $$

therefore, $(S^{\varepsilon } u\hat {)}(w)=\hat {u}(w)\hat {a}_{\varepsilon }(w)$. Now, regarding the principal value of the power function, we write $\hat {a}_{\varepsilon }(w)=w^{2}(\varepsilon +\mathrm {i} w)^{s-2}$. Hence,

$$ \text{Re}\hat{a}_{\varepsilon}(w)>\text{Re} \hat{a}_{0} (w)=\cos(\frac{\pi(2-s)}{2})\vert w\vert^{s}. $$

From (23) and the above equation, we have

$$ \begin{array}{ll} \text{Re}(u,S^{\varepsilon} u)&=\text{Re}({\int}_{-\infty}^{\infty}\vert \hat{u}(w)\vert^{2} \hat{a}_{\varepsilon}(w) \mathrm{d}w)\\ &\geq\cos(\frac{\pi(2-s)}{2}){\int}_{-\infty}^{\infty}\vert \hat{u}(w)\vert^{2} \vert w\vert^{s} \mathrm{d}w. \end{array} $$

(24)

From the contradiction argument investigated in [18, Lemma 4.2], we conclude the coerciveness of the operator $\mathcal {A}.$ According to this result, one can deduce the monotonicity of the operator $\mathcal {A}$ by the Definition 1, i.e.,

$$ \mathcal{A}(u,u-v)-\mathcal{A}(v,u-v)>0. $$

□

Remark 1

It should be noticed that the results of Theorem 4 are well-known [18]; to make the paper self-contained, we have provided another proof based on the interesting results for the Abel integral operator in [13].

In the next theorems, we assert that the results guarantee the existence of a unique weak solution for (1) along with the Riemann-Liouville derivative.

Theorem 5

Suppose that Assumption A holds for the function g(x,u) which satisfies the global (Δ₂)-condition. Consider the variational form (22); then, for every g ∈ U^∗ and 1 < s < 2, (1) has a unique weak solution.

Proof

Let u ∈ U be fixed. Regarding the Lemma 1, $g(.,u(.))\in L_{\tilde {G}}({\Omega }).$ Now, using Hölder inequality, we obtain:

$$ \vert \mathcal{L}(u,v)\vert\leq {\Vert_{0}^{R}}D^{\frac{s}{2}}_{x} u\Vert_{L^{2}({\Omega})}{\Vert_{x}^{R}}D^{\frac{s}{2}}_{1} v\Vert_{L^{2}({\Omega})}+2\Vert g(.,u)\Vert_{\tilde{G}} \Vert v\Vert_{U}. $$

(25)

Applying Theorem 2 and Lemma 1, the above inequality can be simplified as

$$ \vert \mathcal{L}(u,v)\vert\leq \big({\Vert_{0}^{R}}D^{\frac{s}{2}}_{x} u\Vert_{L^{2}({\Omega})}+2\Vert g(.,u)\Vert_{\tilde{G}}\big) \Vert v\Vert_{U}, $$

(26)

which means that the operator ${\mathscr{L}}$ is bounded. Since ${\mathscr{L}}(u,.)$ is linear with respect to the second variable, ${\mathscr{L}}(u,.)\in U^{*}$ for all u ∈ U. Due to the strict monotonicity property of g(x,.) and Theorem 4, one can conclude the monotonicity of the operator ${{\mathscr{L}}}$. With regard to the previous theorem, we have the coercivity of the operator $\mathcal {A}$. Under the assumption about the function g(x,.) and (17), we get that

$$ \lim\limits_{\Vert u\Vert_{G, {\Omega}}\rightarrow \infty} \frac{\rho_{G}(u)}{\Vert u\Vert_{G, {\Omega}}}=\infty, $$

therefore, ${\mathscr{L}}$ is coercive. Here, we want to show the hemi-continuity of the nonlinear monotone operator ${\mathscr{L}}.$ To this end, let $u, w\in \mathbb {R}$ and the sequence t_n tend to 0. The objective is to show that ${\mathscr{L}}(u+t_{n}w,v)$ tends to ${\mathscr{L}}(u,v).$ It is an evident fact that $(_{0}^{R}D^{\frac {s}{2}}_{x} (u+t_{n}w),_{x}^{R}D^{\frac {s}{2}}_{1} v)$ converges to $(_{0}^{R}D^{\frac {s}{2}}_{x} u,_{x}^{R}D^{\frac {s}{2}}_{1}v)$ when t_n tends toward 0. Regarding the continuity of the function g(x,.), the claim on the operator ${\mathscr{L}}$ being hemi-continuous is verified. Consequently, the existence of a unique weak solution is proved by utilizing the Browder-Minty Theorem 1. □

We proceed with the discussion on the regularity of the solution. In order to do this, the following theorem is stated.

Theorem 6

Consider (1) along with the Riemann-Liouville fractional derivative where the function g(x,u) satisfies the Lipschitz condition (19). Then, this equation has a solution which fulfills the nonlinear Volterra-Fredholm integral equation of the form:

$$ u(x)=x^{s-1}\big(_{0}{I^{s}_{x}}(g(.,u(.))-f(.))\big)(1)-~_{0}{I^{s}_{x}}\big(g(.,u(.))-f(.)\big)(x). $$

(27)

In addition, let u ∈ U be a weak solution of (1). Then, $u\in \tilde {U}$ where $\tilde {U}:=\big \lbrace u \in {H}^{\alpha }({\Omega }) \cap \tilde {H}^{\frac {s}{2}}({\Omega }) \mid G(.,u(.))\in L^{1}({\Omega }) \big \rbrace$ for $0\leq \alpha \leq s-\frac {1}{2}.$

Proof

According to the argument about converting the FDEs into integral equations in Chapter 5 of [11] and by adjusting the homogeneous Dirichlet boundary conditions, u(x) takes the following form:

$$ u(x)=w x^{s-1}-\frac{1}{\Gamma(s)}{{\int}_{0}^{x}} (x-y)^{s-1}\big(g(y,u(y))-f(y)\big)\mathrm{d}y, $$

(28)

where $w=\big (_{0}{I^{s}_{x}}(g(.,u(.))-f(.))\big )(1).$ Moreover, by means of Theorem 3 part (c), we have x^s− 1 ∈ H^α(Ω), for $0\leq \alpha <s-\frac {1}{2}.$ On the other hand, Lemma 2 and u ∈ U insure that $g(.,u(.))\in {H}^{\frac {s}{2}}({\Omega })$ which is a subset of L²(Ω). Hence, it achieves that g(.,u(.)) − f(.) ∈ L²(Ω). In addition, we conclude from part (d) of Theorem 3 that $_{0}{I^{s}_{x}}(g(.,u(.))-f(.))\in \tilde {H}^{s}({\Omega })$. Consequently, since $u\in \tilde {H}^{\frac {s}{2}}({\Omega })$, one can deduce that $u\in {H}^{\alpha }({\Omega }) \cap \tilde {H}^{\frac {s}{2}}({\Omega }) ,$ for $0\leq \alpha <s-\frac {1}{2}.$□

Remark 2

Due to the presence of the singular term x^s− 1, it is apparent that the best possible regularity of the solution (1) can occur in H^α(Ω). The similar argument about the regularity of the linear form of (1) reported in Theorem 4.4 and Remark 4.5 of the interesting work [18] verifies the above claim. In that work, ${H}^{\alpha }({\Omega }),~ 0\leq \alpha <s-\frac {1}{2}$ is displayed by H^s− 1+α(Ω), where $1-\frac {s}{2}\leq \alpha <\frac {1}{2}$, to show the presence of the singular term better.

3.2 The Caputo fractional operator

As discussed in [18], the difference between the variational formulation of the Caputo and the Riemann-Liouville equations is their admissible test spaces. It means that the variational formulation of (1) along with the Caputo derivative is:

Find u ∈ U such that

$$ \begin{array}{ll} \mathcal{L}(u,v)&:= \mathcal{A}(u,v)+ \mathcal{B}(u,v)=\langle f,v\rangle,\quad v \in V, \end{array} $$

(29)

where

$$ U:=\lbrace\phi \in \tilde{H}^{\frac{s}{2}}({\Omega})\mid G(.,\phi(.))\in L^{1}({\Omega})\rbrace, $$

and

$$ V:=\lbrace\phi \in \tilde{H}^{\frac{s}{2}}({\Omega})\mid (x^{1-s},\phi)=0 \rbrace. $$

(30)

In order to define the appropriate test space V for the Caputo case, we assume that ϕ^∗(x) = (1 − x)^1−s, which belongs to $\tilde {H}^{\frac {s}{2}}({\Omega })$, and apparently for any ϕ ∈ U, we have $\mathcal {A}(\phi ,\phi ^{*})=0.$ In the Caputo fractional derivative case, we set $V=\text {span} \big \lbrace \tilde {\phi }_{i}(x)=\phi _{i}(x)-\gamma _{i}(1-x)^{s-1}~|~ i=0, \dots , N\big \rbrace$ where

$$ \gamma_{i}=\frac{(x^{1-s},\phi_{i}(x))}{(x^{1-s},(1-x)^{s-1})}, $$

(31)

and ϕ_i ∈ U. We will elucidate the above argument in the next section. Note that both Theorems 4 and 5 are valid for the operators including Caputo derivative which have the same variational formulations for the Riemann-Liouville counterpart. Now, we discuss the regularity of the solution by the following theorem.

Theorem 7

Let us consider (1) with the Caputo fractional derivative for 1 < s < 2 in which the function g(x,u) satisfies the Lipschitz condition (19) and f ∈ H^α(Ω) so that $\alpha + s \in (\frac {3}{2},2)$ and $\alpha \in [0,\frac {1}{2})$. Then, this equation has a solution which fulfills the following nonlinear Volterra-Fredholm integral equation

$$ u(x)= x _{0}{I^{s}_{x}}\big(g(.,u(.))-f(.)\big)(1)-~_{0}{I^{s}_{x}}\big(g(.,u(.))-f(.)\big)(x). $$

(32)

In addition, let u(x) ∈ U be a weak solution of (1). Then, $u\in \tilde {U}$ where $\tilde {U}:=\Big \{ \phi (x) \in {H}^{\alpha +s}({\Omega }) \cap \tilde {H}^{\frac {s}{2}}({\Omega }) \mid G(.,u(.))\in L^{1}({\Omega }) \Big \}$.

Proof

According to [11, Theorem 6. 43], u(x) has the following form:

$$ u(x)=w x-\frac{1}{\Gamma(s)}{{\int}_{0}^{x}} (x-y)^{s-1}\big(g(y,u(y))-f(y)\big)\mathrm{d}y, $$

(33)

where $w=\big (_{0}{I^{s}_{x}}(g(.,u(.))-f(.))\big )(1)$ is determined by adjusting the boundary condition. On the other hand, Lemma 2 and u ∈ U imply that $g(.,u(.))\in {H}^{\frac {s}{2}}({\Omega })$ which is a subset of H^α(Ω). Therefore, g(.,u(.)) − f(.) ∈ H^α(Ω). Hence, by part (d) of Theorem 3, we have $_{0}{I^{s}_{x}}(g(.,u(.))-f(.))\in {H}^{\alpha +s}({\Omega })$. Moreover, by means of Theorem 3 part (c), we have $x\in \tilde {H}^{\beta }({\Omega }),$ for $0\leq \beta <\frac {3}{2}.$ Consequently, from the inclusion argument and $u\in \tilde {H}^{\frac {s}{2}} ({\Omega })$, we can conclude that $u\in {H}^{\alpha +s}({\Omega }) \cap \tilde {H}^{\frac {s}{2}}({\Omega })$ where $0\leq \alpha <\frac {1}{2}.$□

Remark 3

As observed in Theorems 6 and 7, owing to the existence of the intrinsic singular term x^s− 1 in the solution representation, the solution of the differential equation with the Riemann-Liouville derivative has less regularity in comparison with the Caputo fractional counterpart. In fact, the best possible regularity in the Riemann-Liouville fractional derivative case belongs to $\tilde {H}^{s-1+\alpha }({\Omega })$. It is worthy to note that the superiority for the Caputo fractional derivative case comes from the fact that the function under the Caputo derivative is supposed to be twice differentiable.

In the following remark, the stability of the variational formulations is shown.

Remark 4

If we assume that $f \in \tilde { H}^{-\frac {s}{2}}({\Omega }) \hookrightarrow U^{*}$ and take u = v in the relation (22), then by using the results of Theorem (4) and the relation ${\mathscr{B}}(u,u) \geq 0$, we see that there is a constant c > 0 such that $c \Vert u\Vert ^{2}_{\tilde { H}^{\frac {s}{2}}({\Omega })} \leq \mathcal {A}(u,u) + {\mathscr{B}}(u,u)$. Note that $\langle f,u \rangle _{U^{*}, U} =\langle f,u \rangle _{\tilde { H}^{-\frac {s}{2}}, \tilde { H}^{\frac {s}{2}}}$, so we get that

$$ c \Vert u\Vert^{2}_{\tilde{ H}^{\frac{s}{2}}({\Omega})} \leq \vert \langle f,u \rangle \vert\leq \Vert f\Vert_{\tilde{ H}^{-\frac{s}{2}}({\Omega})} \Vert u \Vert_{\tilde{ H}^{\frac{s}{2}}({\Omega})}, $$

which means that there is a constant C > 0 such that $\Vert u\Vert ^{2}_{\tilde { H}^{\frac {s}{2}}({\Omega })} \leq C \Vert f\Vert _{\tilde { H}^{-\frac {s}{2}}({\Omega })}$.

4 Finite element approximation

To find an approximate solution, we discretize the continuous problem (22) by a Galerkin finite element method. In order to do this, a piecewise polynomial finite element method is introduced over the interval Ω = [0,1]. Let us define $\mathbb {P}_{r}({\Omega })$ as the space of univariate polynomials of the degree less than or equal to r, for positive integer r and χ_h be a uniform mesh partition on Ω, given by

$$ 0=x_{0}<x_{1}<\dots<x_{N-1}<x_{N}=1, \quad N \in \mathbb{N}, $$

(34)

with fixed mesh size h_i = x_i − x_i− 1. The set χ_h induces a mesh $\mathcal {T}_{h}=\{\tau _{i} | 1\leq i \leq N\}$ on Ω, where τ_i = [x_i− 1,x_i]. The length of a subinterval $\tau \in \mathcal {T}_{h}$ is denoted by h_τ and the maximal mesh width by $h:=\max \limits \left \{ h_{\tau }:\tau \in \mathcal {T}_{h}\right \}$. We choose standard continuous and piecewise polynomial function space of the degree $r\in \mathbb {N}$ on [0,1] defined by

$$ S_{\mathcal{T}}^{r}({\Omega}):= \{v\in C({\Omega}):\left. v\right\vert_{\tau}\in\mathbb{P}_{r}\left( \tau\right) ,\forall\tau\in\mathcal{T}\}. $$

(35)

The nodal points are given by

$$ \mathcal{N}_{r}:= \left\{ \xi_{i,j}:=x_{i-1}+j\frac{x_{i}-x_{i-1}}{r}\quad1\leq i\leq N\text{, }0\leq j\leq r-1\right\} \cup\left\{ 1\right\}. $$

We choose the usual standard Lagrange basis functions $b_{i,j}^{(r)}$ of $S_{\mathcal {T}}^{r}({\Omega })$. Now, with these piecewise functions, one can define the discrete admissible space which is a subspace of $S^{r}_{\mathcal {T}}({\Omega })\bigcap {H^{1}_{0}}({\Omega })$ denoted by A_h. Particularly, we focus on the linear elements in the numerical experiments. Let $\mathcal {I}_{h}$ be the Lagrange interpolation operator mapping into A_h. We denote the finite element test and trial spaces U_h for the Riemann-Liouville fractional derivative with A_h which is described above. In order to investigate the Caputo fractional derivative case, we consider finite-dimensional set U_h = A_h as the trial space. In addition, to construct a suitable test space V_h, let $V_{h}=\text {span} \Big \lbrace \tilde {\phi }_{i}(x) ~|~ i=0, 1, \dots , N \Big \rbrace$ where

$$ \tilde{\phi}_{i}(x)=\phi_{i}(x)-\gamma_{i}(1-x)^{s-1}, $$

and γ_i is given by (31). Finally, the discrete variational formulation released from (22) and (29) is:

Find u_h ∈ A_h such that

$$ \mathcal{L}(u_{h},v_{h})=F(v_{h}), \quad \forall v_{h}\in V_{h}. $$

(36)

We notice that in the approximation procedure, one can use the property (f ) of Theorem 3 which means that for the computation of $_{0}^{R} D_{x}^{\frac {s}{2}}u_{h}$, one can utilize the relation:

$$ \begin{array}{@{}rcl@{}} {~}_{0}^{R} D_{x}^{\frac{s}{2}} \phi_{i} =_{0}^{R} I_{x}^{1-\frac{s}{2}} \phi_{i}&=&\frac{1}{\Gamma(1-\frac{s}{2})}{{\int}_{0}^{x}}(x-t)^{-\frac{s}{2}}\phi_{i}^{\prime}(t)\mathrm{d}t\\ &=&\frac{1}{\Gamma(1-\frac{s}{2})}{{\int}_{0}^{x}}(x-t)^{-\frac{s}{2}}(\frac{\chi_{[x_{i-1},x_{i}]}}{h_{i}}-\frac{\chi_{[x_{i},x_{i+1}]}}{h_{i+1}} )\mathrm{d}t\\ &=&\frac{1}{\Gamma(1-\frac{s}{2})}\Big[h^{-1}_{i}\big((x-x_{i-1})_{+}^{1-\frac{s}{2}} -(x-x_{i})_{+}^{1-\frac{s}{2}} \big)\\ &&- h^{-1}_{i+1}\big((x-x_{i})_{+}^{1-\frac{s}{2}} -(x-x_{i+1})_{+}^{1-\frac{s}{2}} \big)\Big], \end{array} $$

(37)

where $a_{+}=\max \limits \{a,0\}$, and analogously for $_{x}^{R} D_{1}^{\frac {s}{2}}u$, we apply

$$ \begin{array}{@{}rcl@{}} {~}_{x}^{R} D_{1}^{\frac{s}{2}} \phi_{i} =-_{x}^{R} I_{1}^{1-\frac{s}{2}} \phi_{i}&=&- \frac{1}{\Gamma(1-\frac{s}{2})}{{\int}_{x}^{1}}(x-t)^{-\frac{s}{2}}\phi_{i}^{\prime}(t)\mathrm{d}t\\ &=&\frac{1}{\Gamma(1-\frac{s}{2})}{{\int}_{x}^{1}}(x-t)^{-\frac{s}{2}}(\frac{\chi_{[x_{i-1},x_{i}]}}{h_{i}}-\frac{\chi_{[x_{i},x_{i+1}]}}{h_{i+1}} )\mathrm{d}t\\ &=&\frac{1}{\Gamma(1-\frac{s}{2})}\Big[h^{-1}_{i}\big((x_{i}-x)_{+}^{1-\frac{s}{2}} -(x_{i-1}-x)_{+}^{1-\frac{s}{2}} \big)\\ &&- h^{-1}_{i+1}\big((x_{i+1}-x)_{+}^{1-\frac{s}{2}} -(x_{i}-x)_{+}^{1-\frac{s}{2}} \big)\Big]. \end{array} $$

(38)

Therefore, the term $\mathcal {A}(\phi _{i}, \phi _{j})= \big (_{0}^{R} D_{x}^{\frac {s}{2}}\phi _{i}, _{x}^{R} D_{1}^{\frac {s}{2}}\phi _{j} \big )$ can be derived by the above arguments for the Riemann-Liouville derivative case. For the Caputo fractional derivative, we have $\mathcal {A}(\phi _{i},\tilde { \phi _{j}})= \big (_{0}^{R} D_{x}^{\frac {s}{2}}\phi _{i}, _{x}^{R} D_{1}^{\frac {s}{2}}\tilde { \phi _{j}} \big )$ which can be simplified as

$$ \big(_{0}^{R} D_{x}^{\frac{s}{2}}\phi_{i} , _{x}^{R} D_{1}^{\frac{s}{2}}\phi_{j} \big) - \gamma_{j} \big(_{0}^{R} D_{x}^{\frac{s}{2}}\phi_{i} , _{x}^{R} D_{1}^{\frac{s}{2}}(1-x)^{s-1} \big). $$

The first term can be computed using the relations (37) and (38) and the second term is disappeared because

$$ \begin{array}{ll} \big({~}_{0}^{R} D_{x}^{\frac{s}{2}}\phi_{i} , _{x}^{R} D_{1}^{\frac{s}{2}}(1-x)^{s-1}\big)& = - \big(_{0} I_{x}^{1-\frac{s}{2}}\phi^{\prime}_{i} , _{x}^{R} D_{1}^{\frac{s}{2}}(1-x)^{s-1}\big)\\ &= c_{\alpha} \big(\phi^{\prime}_{i} , _{0} I_{x}^{1-\frac{s}{2}}(1-x)^{\frac{s}{2}-1}\big)\\ &= c_{\alpha}\big(\phi^{\prime}_{i} , 1\big)\\ & =0, \end{array} $$

(39)

where c_α is a constant depending on α.

4.1 Convergence analysis

This section is devoted to the study of the approximate solution achieved in the previous section. For this end, we consider the existence and uniqueness issue for the discrete equation. In addition, we find an appropriate priori error bound.

Theorem 8

The discrete problem (36) has a unique solution.

Proof

The existence of the discrete solution can be verified through the Browder-Minty theorem with the same argument pursued in Section 3. For the uniqueness issue, let u₁ and u₂ be finite element solutions of (22). Hence,

$$ \begin{array}{@{}rcl@{}} 0= \mathcal{L}(u_{1},v_{h})- \mathcal{L}(u_{2},v_{h})&=&\mathcal{A}(u_{1},v_{h})-\mathcal{A}(u_{2},v_{h})+\mathcal{B}(u_{1},v_{h})-\mathcal{B}(u_{2},v_{h})\\ &=&-{\int}_{\Omega} {~}_{0}{~}^{R}D_{x}{~}^{\frac{s}{2}} (u_{1}-u_{2})(x) _{x}^{R}D_{1}{~}^{\frac{s}{2}}v_{h}(x)\mathrm{d}x\\ &&+{\int}_{\Omega} \left( g(x,u_{1}(x))-g(x,u_{2}(x))\right)v_{h}(x)\mathrm{d}x. \end{array} $$

(40)

Since the operator $\mathcal {A}(u,v)$ is coercive, for v = u₁ − u₂, we get that

$$ - {\int}_{\Omega} {~}_{0}{~}^{R}D_{x}{~}^{\frac{s}{2}} (u_{1}-u_{2})(x) _{x}{~}^{R}D_{1}{~}^{\frac{s}{2}}(u_{1}-u_{2})(x)\mathrm{d}x\geq c\Vert u_{1}-u_{2} \Vert^{2}_{\tilde{H}{~}^{\frac{s}{2}}({\Omega})}. $$

Using the above equation and (40), one can conclude that

$$ \begin{array}{@{}rcl@{}} c\Vert u_{1}-u_{2} \Vert^{2}_{\tilde{H}^{\frac{s}{2}}({\Omega})}&+{\int}_{\Omega} \left( g(x,u_{1}(x))-g(x,u_{2}(x))\right)(u_{1}(x)-u_{2}(x))\mathrm{d}x\\ &\leq - {\int}_{\Omega}{~}_{0}^{R}D_{x}^{\frac{s}{2}} (u_{1}-u_{2})(x) _{x}^{R}D_{1}^{\frac{s}{2}}(u_{1}-u_{2})(x)\mathrm{d}x\\ &+{\int}_{\Omega} (g(x,u_{1}(x))-g(x,u_{2}(x)))(u_{1}(x)-u_{2}(x))\mathrm{d}x\\ &= 0. \end{array} $$

Since g(x,u) is a monotone function with respect to the second variable, thereby the above inequality, we conclude the uniqueness of the approximate solution. □

Lemma 3

Let $\mathcal {T}_{h}$ be a uniform mesh on Ω. For real numbers s,m with $m\geq \frac {s}{2}$, and also $S_{\mathcal {T}}^{r} ({\Omega })$ with an integer r ≥ 0, we define $\hat {r}=\min \limits \lbrace r+1,m\rbrace -\frac {s}{2}$. Then, there is a constant c > 0 depending on s, m, r, and $\mathcal {T}_{h}$ such that

$$ \min\limits_{v_{h} \in S_{\mathcal{T}}^{r}({\Omega})} \Vert u-v_{h} \Vert_{H^{\frac{s}{2}}({\Omega})}\leq ch^{\hat{r}}\Vert u\Vert_{H^{m}({\Omega})}, $$

(41)

for all $u\in H^{\frac {s}{2}}({\Omega })\cap H^{m}({\Omega }).$ Particularly, for r = 1 and $\phi \in H^{\frac {s}{2}}({\Omega })\cap H^{\gamma }({\Omega })$ where $\gamma ={\min \limits } \lbrace 2,m \rbrace$ and $\hat {r}=\gamma -\frac {s}{2}$, one can deduce that

$$ \min\limits_{v_{h} \in U_{h}} \Vert \phi-v_{h} \Vert_{H^{\frac{s}{2}}({\Omega})}\leq ch^{\gamma-\frac{s}{2}}\Vert \phi \Vert_{H^{\gamma}({\Omega})}. $$

(42)

Moreover, if ϕ ∈ H^γ(Ω) ∩ V where V is defined in (30), the following relation holds

$$ \min\limits_{v_{h} \in V_{h}} \Vert \phi-v_{h} \Vert_{H^{\frac{s}{2}}({\Omega})}\leq ch^{\gamma-\frac{s}{2}}\Vert \phi\Vert_{H^{\gamma}({\Omega})}, $$

(43)

Proof

The relation

$$ \inf_{v\in U_{h}}\Vert u-v\Vert_{H^{\alpha}({\Omega})}\leq \Vert u- \mathcal{I}_{h}u\Vert_{H^{\alpha}({\Omega})}, \quad 0\leq \alpha \leq 1, $$

and the similar argument for finding an error estimation of the standard Lagrange finite element for the integer order Sobolev space lead to an error bound for the interpolation error in the intermediate spaces (41) and the special cases (42) and (43); for more details, see [8, 18]. □

Theorem 9

Assume that u is the exact solution of (1) and u_h is the approximate solution of the variational formulation (22) or (29). Then

$$ \Vert u-u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}\leq Ch^{\gamma-\frac{s}{2}} \Vert u\Vert_{H^{\gamma}({\Omega})}, $$

(44)

where γ differs for the nonlinear boundary value problems with Caputo or Riemann-Liouville fractional derivative. For the case of Caputo differential operator, γ is equal to s. In addition, γ belongs to the interval $[\frac {s}{2}, s-\frac {1}{2}]$ for the Riemann-Liouville fractional counterpart.

Proof

Consider u_h ∈ U_h is the solution of finite element space of (1) which satisfies the following formulation:

$$ \mathcal{A}(u_{h},v_{h})+\mathcal{B}(u_{h},v_{h})=\langle f,v_{h}\rangle, \quad \quad v_{h}\in V_{h}. $$

Next, by subtracting the above equation and (22), we get that

$$ \mathcal{A}(u,v)-\mathcal{A}(u_{h},v_{h})+\mathcal{B}(u,v)-\mathcal{B}(u_{h},v_{h})=\langle f,v\rangle-\langle f,v_{h}\rangle. $$

(45)

Consider the projection operator $\mathcal {P}_{h}:H^{\frac {s}{2}}({\Omega })\rightarrow U_{h}$ defined by

$$ \mathcal{A}(u,v_{h})=\mathcal{A}(\mathcal{P}_{h}u,v_{h}). $$

(46)

Now by adding and subtracting $\mathcal {P}_{h}u$, we have

$$ u-u_{h}=(u-\mathcal{P}_{h}u)+(\mathcal{P}_{h}u-u_{h}):=\xi+\eta. $$

(47)

Then, (45) can be rewritten as follows:

$$ \mathcal{A}(u,v)-\mathcal{A}(\mathcal{P}_{h}u,v_{h})+\mathcal{A}(\mathcal{P}_{h}u,v_{h})-\mathcal{A}(u_{h},v_{h})+\mathcal{B}(u,v)-\mathcal{B}(u_{h},v_{h})=\langle f,v\rangle-\langle f,v_{h}\rangle, $$

(48)

therefore, from (46) and setting v = v_h, we get that

$$ \mathcal{A}(\mathcal{P}_{h}u,v_{h})-\mathcal{A}(u_{h},v_{h})+\mathcal{B}(u,v_{h})-\mathcal{B}(u_{h},v_{h})=0, $$

(49)

or, regarding the bilinearity of the operator $\mathcal {A}$,

$$ \mathcal{A}(\mathcal{P}_{h}u-u_{h},v_{h})+\mathcal{B}(u,v_{h})-\mathcal{B}(u_{h},v_{h})=0. $$

(50)

Letting v_h = η, we have

$$ \mathcal{A}(\eta,\eta)=\mathcal{B}(u_{h},\eta)-\mathcal{B}(u,\eta). $$

(51)

Next, by the coercivity of $\mathcal {A}$, there is a constant c₀ such that

$$ \mathcal{A}(\eta, \eta)\geq c_{0}\Vert \eta \Vert^{2}_{\tilde{H}^{\frac{s}{2}}({\Omega})}. $$

(52)

On the other hand,

$$ \begin{array}{ll} \vert\mathcal{B}(u_{h},\eta)-\mathcal{B}(u,\eta)\vert &=\vert\big(g(.,u)-g(.,u_{h}), \eta\big)\vert\\ &=\vert\Big(g(.,u)-g(.,\mathcal{P}_{h}u)+g(.,\mathcal{P}_{h}u)-g(.,u_{h}), \eta\Big)\vert\\ &\leq \Vert g(.,u)-g(.,\mathcal{P}_{h}u)\Vert \Vert \eta \Vert+\Vert g(.,\mathcal{P}_{h}u)-g(.,u_{h})\Vert\Vert \eta \Vert\\ &\leq l_{M}\Vert \eta \Vert_{H^{\frac{s}{2}}({\Omega})}(\Vert \xi \Vert_{H^{\frac{s}{2}}({\Omega})} +\Vert \eta \Vert_{H^{\frac{s}{2}}({\Omega})}), \end{array} $$

(53)

where l_M is verified in (19). If η≠ 0 and c₀ > l_M, then one can deduce from (52) and (53) that

$$ \Vert \eta\Vert_{H^{\frac{s}{2}}({\Omega})} \leq \frac{l_{M}}{c_{0}-l_{M}}\Vert \xi\Vert_{H^{\frac{s}{2}}({\Omega})}. $$

(54)

Consequently by (47), (54), and Lemma 3, we get that

$$ \Vert u-u_{h} \Vert_{H^{\frac{s}{2}}({\Omega})} \leq C h^{\gamma-\frac{s}{2}}\Vert u\Vert_{H^{\gamma}({\Omega})}. $$

(55)

It is conspicuous that for each derivative case with different regularity, γ should be different. Using the regularity Theorems 6 and 7 for f ∈ L²(Ω), one can deduce that $\gamma \in [\frac {s}{2}, s-\frac {1}{2}]$ for the Riemann-Liouville fractional derivative operator, and γ = s in the Caputo fractional one. □

Remark 5

Let f ∈ H^α(Ω) where α varies in the interval $\left [0,\frac {1}{2}\right )$ and $\alpha +s>\frac {3}{2}$. If α + s < 2, then Theorem 7 yields that the γ arising in the error estimation given for the Caputo fractional operator is equal to α + s which means that $\Vert u-u_{h}\Vert _{H^{\frac {s}{2}}({\Omega })}=\mathcal {O}(h^{\alpha +\frac {s}{2}})$. Otherwise, γ = 2 and $\Vert u-u_{h}\Vert _{H^{\frac {s}{2}}({\Omega })}=\mathcal {O}(h^{2-\frac {s}{2}})$.

We finalize this section with a remark on the stability of the Galerkin finite element method.

Remark 6

By means of the triangular inequality, it is easily seen that

$$ \Vert u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}\leq \Vert u-u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}+\Vert u\Vert_{H^{\frac{s}{2}}({\Omega})}. $$

Using Theorem 2, we obtain that

$$ \Vert u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}\leq Ch^{\gamma-\frac{s}{2}}\Vert u\Vert_{H^{\gamma}({\Omega})}+\Vert u\Vert_{H^{\frac{s}{2}}({\Omega})}, $$

(56)

where C is a constant and γ is equal to s, if one considers the Caputo differential operator and γ belongs to $[\frac {s}{2}, s-\frac {1}{2}]$ for the Riemann-Liouville fractional operator case. Therefore, we get an upper bound for the approximate solution in the case of Caputo fractional derivative:

$$ \Vert u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}\leq Ch^{\frac{s}{2}}\Vert u\Vert_{H^{s}({\Omega})}+\Vert u\Vert_{H^{\frac{s}{2}}({\Omega})}, $$

(57)

and the following one in the case of Riemann-Liouville fractional derivative operator

$$ \Vert u_{h}\Vert_{H^{\frac{s}{2}}({\Omega})}\leq C\Vert u\Vert_{H^{\frac{s}{2}}({\Omega})}. $$

(58)

So, the final results on the stability discussion are obtained by the argument given in Remark 4.

5 Numerical illustrations

The numerical experiments are employed to exhibit the applicability of the Galerkin finite element method for the fractional nonlinear boundary value problems with Caputo and Riemann-Liouville derivatives. The experiments are implemented in $\textit {Mathematica}^{{\circledR }}$ software platform. We report the absolute error along with the numerical and theoretical rates of convergence for some examples which satisfy the assumptions considered in the previous sections. Furthermore, the numerical algorithm is examined for some examples with the absence of the mentioned assumptions.

In general, for the numerical experiment with the Galerkin method, one of the main issues is the approximation of the integrals. In our examination, the Galerkin finite element solution is obtained from the fully discrete weak from:

Find u_h ∈ A_h such that

$$ \mathcal{L}_{h}(u_{h},v_{h})=F(v_{h}), \quad \forall v_{h}\in V_{h}. $$

(59)

In the above nonlinear system of equations, we have utilized Gauss-Kronrod quadrature formula to compute the integrals that need to be evaluated numerically. This happens mainly for the integrals involving the nonlinear term. Furthermore, to solve the nonlinear system, we employ the Newton iteration method. In order to do this, consider the bilinear form N_h(u_h;.,.) defined on $S^{r}_{\mathcal {T}}({\Omega }) \times S^{r}_{\mathcal {T}}({\Omega })$ by

$$ N_{h}(u_{h};w_{h},v_{h})=(_{0}^{R}D^{\frac{s}{2}}_{x} w_{h},~_{x}^{R}D^{\frac{s}{2}}_{1} v_{h})+(g_{u}(.,u_{h})w_{h},v_{h}). $$

The Newton method for approximating u_h by a sequence $\{{u^{k}_{h}}\}_{k \in \mathbb {N}}$ in $S^{r}_{\mathcal {T}}({\Omega })$ can be written as

$$ N_{h}({u^{k}_{h}};u^{k+1}_{h}-{u^{k}_{h}},v_{h})=F(v_{h})-\mathcal{L}_{h}({u^{k}_{h}},v_{h}), \quad \forall v_{h}\in S^{r}_{\mathcal{T}}({\Omega}), $$

where ${u^{0}_{h}}\in S^{r}_{\mathcal {T}}({\Omega })$ is an initial guess chosen by the steepest gradient algorithm.

Through tables and figures, we notify the experimental and the possible theoretical convergence rates (the reported numbers in the parentheses) for the finite element approximation of the nonlinear boundary value problems with the Riemann-Liouville and Caputo fractional derivatives in L² and $H^{\frac {s}{2}}$-norms.

Example 1

Consider the fractional derivative (1) with g(x,u(x)) = 3xu³(x). The right-hand side function f(x) is chosen such that

• (a) The exact solution with Riemann-Liouville fractional derivative is $u(x)= \frac {1}{\Gamma (s+1)}(x^{s-1}-x^{s})$;

• (b) The exact solution is $u(x)=\frac {1}{\Gamma (s+1)}(x-x^{s})$, where the derivative operator is of Caputo type.

This problem satisfies the assumptions introduced in Section 4.1. Therefore, the convergence rates of the Caputo and the Riemann derivative cases in term of $H^{\frac {s}{2}}$ error norm are $O(h^{\frac {s}{2}})$ and $O(h^{\frac {s-1}{2}})$, respectively. This argument directly follows from Theorem 9 for the function f(x) that belongs to L²(Ω). Tables 1 and 2 report the L² and $H^{\frac {s}{2}}$ error norms for different s ∈ (1,2). Moreover, we have provided some figures to exhibit both the theoretical and the practical rates of convergence. Figures 1 and 2 display the absolute errors for the Caputo and the Riemann-Liouville fractional derivatives with the above nonlinear term. In figures, the dashed lines show the theoretical convergence rate.

Table 1 The absolute error in L² and $H^{\frac {s}{2}}$-norms for different value $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for Example 1 with the Caputo fractional derivative operator

Table 2 The absolute error in L² and $H^{\frac {s}{2}}$-norms for different values $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for Example 1 with the Riemann-Liouville fractional operator

Fig. 1

Plots of the absolute error in L² and $H^{\frac {s}{2}}$-norms in logarithmic scale for Example 1 with the Caputo fractional derivative

Fig. 2

Plots of the absolute error in L² and $H^{\frac {s}{2}}$-norms in logarithmic scale for Example 1 with the Riemann-Liouville fractional derivative

Example 2

In this example, we discuss the approximation of (1) with $g(x,u(x))=\sin \limits (x)u^{5}(x)$. The right-hand side function f(x) is chosen such that

• (a) The exact solution is $u(x)= \frac {\Gamma (1.5)}{\Gamma (s+1.5)}(x^{s-1}-x^{s+0.5})$ for the Riemann-Liouville case.

• (b) The exact solution is $u(x)= \frac {\Gamma (1.5)}{\Gamma (s+1.5)}(x-x^{s+0.5})$ when (1) entailes the Caputo fractional derivative.

By a similar reasoning as Example 1, it is seen that the assumptions discussed in the theoretical parts hold. Therefore, we expect $O(h^{\frac {s}{2}})$ and $O(h^{\frac {s-1}{2}})$ convergence rates for the Caputo and the Riemann-Liouville fractional differential operators, respectively. This claim is verified by the numerical results reported in Tables 3 and 4 which exhibit the absolute errors in L² and $H^{\frac {s}{2}}$-norms for different s ∈ (1,2).

Table 3 The absolute error in L² and $H^{\frac {s}{2}}$-norms for different values $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for Example 2 with the Caputo fractional operator

Table 4 The absolute errors in L² and $H^{\frac {s}{2}}$-norms for different values $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for Example 2 with the Riemann-Liouville fractional operator

Example 3

Consider the nonlinear Riemann-Liouville fractional differential (1) with $g(x,u(x))=x\exp (u(x))$. The right-hand side function f(x) is chosen such that the exact solution u(x) is

$$ u(x)= \frac{1}{\Gamma(s+2)}(x^{s-1}-x^{s+1})-\frac{2}{\Gamma(s+3)}(x^{s-1}-x^{s+2}). $$

Table 5 reports the absolute error in L² and $H^{\frac {s}{2}}$-norms for different s ∈ (1,2) with the above nonlinear term.

Table 5 The absolute error in L² and $H^{\frac {s}{2}}$-norms for different values $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for Example 3 with the Riemann-Liouville fractional derivative

Example 4

We present the nonlinear Caputo fractional differential equation (1) with g(x,u(x)) = (u(x) − x)², where the exact solution is $u(x)=\frac {\Gamma (\frac {3}{4})}{\Gamma (s+\frac {3}{4})}(x-x^{s-\frac {1}{4}})$.

Table 6 reports the absolute errors in L² and $H^{\frac {s}{2}}$-norms for different s ∈ (1,2) with the above nonlinear term.

Table 6 The absolute error in L² and $H^{\frac {s}{2}}$-norms for different values $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$ and mesh size $h=\frac {1}{2^{k}\times 10}$ for the Caputo fractional operator for Example 4

Example 5

As the final example, we deal with the linear Caputo fractional differential equation (g(x,u(x)) = 0) by considering f(x) = x^𝜃 to belong H^α(Ω) for $\alpha \in \left [ 0,\theta +\frac {1}{2}\right )$ and $\theta \in \{\frac {-1}{3},\frac {-1}{4},\frac {-1}{5} \}.$ The exact solution for different values of 𝜃 is u(x) = c_𝜃(x^s− 1 − x^s+𝜃) with $c_{\theta }=\frac {\Gamma (\theta +1)}{\Gamma (s+\theta +1)}.$

Table 7 displays the theoretical and numerical rates of convergence in $H^{\frac {s}{2}}$-norm for different $s=\frac {7}{4},\frac {3}{2},\frac {4}{3}$. This problem is satisfied the Remark 5. For different values of α and s, the value of γ is given by $\min \limits \{\alpha +s,2\}$. For instance, for $s=\frac {7}{4}$ and $\theta =\frac {-1}{5},$ then γ = 2 and the rate of convergence is $O(h^{2-\frac {s}{2}})=O(h^{1.125}).$

Table 7 A comparison between theoretical and numerical convergence rates in $H^{\frac {s}{2}}$-norm for Example 5 with the Caputo fractional derivative

6 Conclusion and future studies

In this paper, we have studied the Lagrange finite element method for a class of semi-linear FDEs of the Riemann-Liouville and the Caputo types. To this aim, a weak formulations of the problems have been introduced in the suitable function spaces constructed by considering the fractional Sobolev and Musielak-Orlicz spaces due to the presence of the nonlinear term. In addition, the existence and uniqueness issues of the weak solution together with its regularity are discussed. The weak formulation is discretized by the Galerkin method with piecewise linear polynomials basis functions. Finding an error bound in $H^{\frac {s}{2}}$-norm is considered for the Riemann-Liouville and Caputo fractional differential equations. Different examples with the varieties of the nonlinear terms have been examined and the absolute errors are reported in L² and $H^{\frac {s}{2}}$-norms.

The nature of the nonlinearity and also the fractional essence of the problem cause low-order convergence of the method. In order to improve the approach for this class of FDEs, one can apply the idea of splitting method, where the solution is separated into regular and singular parts; this is possible by utilizing the Taylor expansion of the nonlinear operator and the finite element method accompanied by a quasi-uniform mesh. Also, as discussed in the numerical experiments section, the integrals in the obtained nonlinear system are discretized by a suitable quadrature method. Surveying the effect of quadrature method in finite element approximation and a priori error estimation is an idea for the future studies. As reported in the numerical section, we have observed the absolute errors in L²-norm which are sharper than the errors in $H^{\frac {s}{2}}$-norm. An interesting question for a further study is how to obtain an appropriate theoretical error bound in L²-norm.