January 01, 1970
Fractional boundary value problems are often used to model complex systems and processes characterized by memory effects and anomalous diffusion. In this paper, we consider fractional boundary value problems involving the Riesz–Caputo operator, which is particularly suited for modeling physical phenomena exhibiting symmetric diffusive effects. We provide an integral representation of the solution to prove existence and uniqueness of the fractional differential problem. We introduce a B-spline collocation method to approximate the solution of the problem and provide a convergence analysis, with both theoretical insights and numerical experiments.
Keywords: Fractional boundary value problem, Riesz–Caputo operator, B-splines, Collocation method, Convergence analysis.
MSC 2020: 26A33, 35R11, 65M70.
Fractional differential equations extend the classical theory of differential equations by incorporating derivatives of non-integer order, providing a powerful tool for modeling memory and hereditary effects in various physical and engineering problems.
The growing interest in fractional calculus comes from its ability to more accurately describe complex phenomena where the traditional integer-order models would fail. In contrast to ordinary differential equations, which imply local interactions,
fractional differential equations show the ability to capture long-range dependencies and anomalous diffusion, making them highly relevant in fields such as viscoelasticity, fluid dynamics, control theory, signal processing, finance, and many others [1]–[5].
In this paper we will study nonlinear fractional boundary problems that have the Riesz-Caputo derivative as the fractional operator. Among various fractional operators, the Riesz-Caputo operator plays a crucial role due to its unique ability to model
diffusion processes with directional symmetry [6], such as anomalous transport in porous media and turbulence in fluid mechanics. The symmetric property of the
Riesz-Caputo operator was also exploited in [7] to efficiently evaluate the Riesz-Caputo derivative of B-splines, and later applied in [8] to the numerical solution of a linear convection-diffusion problem. Fractional differential equations involving the Riesz-Caputo operator are, however, a rather new research
subject, so analytical results on existence and uniqueness of their solutions as well as numerical methods for their approximations are yet limited. In [9], [10], comprehensive studies of differential problems with the Riemann-Liouville and the Caputo fractional derivatives are provided. Pedas and Tamme have studied both theoretically and numerically different
types of linear fractional differential problems based on the Caputo operator, see [11]–[13] and references there in. Few years later, the same authors dealt with the nonlinearity for a fractional initial value problem [14]. The general theory used to prove the convergence of the numerical method is the one originally developed by Vainikko [15], [16]. Non linear fractional differential problems based on the classical Caputo operator have been considered also in [17], [18]. Our aim is to extend these results to the Riesz-Caputo operator.
The paper is organized as follows: In Section 2 we introduce the mathematical model along with some useful definitions and few fundamental Lemmas describing analytical properties of the fractional integral and differential
operators involved. In the same section we also prove existence and uniqueness of the solution under suitable hypotheses. In Section 3 we describe the numerical method and provide a formal convergence analysis. The theoretical
findings are also demonstrated by numerical experiments in Section 4. Finally in Section 5 we outline the conclusions.
We consider the fractional boundary value problem: \[\label{eq:main95eq}
\left \{ \begin{array}{ll} D_{RC}^\alpha u(x)=f(x,u(x)), \;\;\;x \in [0,L], \\ \\
a_i u^{(i)}(0)+b_iu^{(i)}(L)=c_i, \;\;i=0,..,{\lfloor \alpha \rfloor} ,
\end{array}
\right.\tag{1}\] where \(1<\alpha<2\). Here, \(D_{RC}^\alpha u(x)\) denotes the Riesz-Caputo fractional derivative defined as [9], [19], [20] \[\label{eq:RCdef}
D_{RC}^\alpha u(x) \mathrel{\vcenter{:}}= \frac{1}{2}(D^\alpha_{0,x} + (-1)^m \;D^\alpha_{x,L})u(x),\tag{2}\] where the left and right Caputo derivatives are defined by: \[\begin{align}
\tag{3}
D^\alpha_{0,x} u(x) &\mathrel{\vcenter{:}}= \frac{1}{\Gamma(\tilde{m})} \int_0^x \frac{u^{(m)}(\xi)}{(x-\xi)^{1-\tilde{m}}} d\xi, \\
\tag{4}
D^\alpha_{x,L} u(x) &\mathrel{\vcenter{:}}= \frac{(-1)^m }{\Gamma(\tilde{m})} \int_x^L \frac{u^{(m)}(\xi)}{(\xi-x)^{1-\tilde{m}}} d\xi,
\end{align}\] with \(m= \lceil \alpha \rceil\), \(\tilde{m}=m-\alpha\) and
\(\Gamma\) the Euler gamma function.
Note that in definition 2 we use the constant \(\frac{1}{2}\) in front of the combination of the left and right Caputo operator because in this case the Riesz-Caputo derivative agrees with the
ordinary derivative when \(\alpha\) is an integer. However other choices can be found in literature (see e.g. [21]).
We start by analyzing some properties of the Riesz-Caputo operator.
Lemma 1. Let \(\alpha \in \mathbb{R}\), with \(1<\alpha<2\) and let \(D_{RC}^\alpha (\cdot)\) denote the Riesz-Caputo operator defined in 2 . Let \(p_i(x)= k x^i\), \(i=0,1\), \(k \in \mathbb{R}\). Then, \[\label{eq:prop95constant} D_{RC}^\alpha (p_i)=0, \qquad i=0,1.\qquad{(1)}\]
Proof. To compute the Riesz-Caputo derivative, we have to evaluate the integrals in eqs. 3 4 . Since \(1<\alpha<2\), in the integrand the second derivative appears, which is zero for both \(p_0(x)\) and \(p_1(x)\). ◻
Now, we introduce the fractional left and right Riemann-Liouville integrals of order \(\mu\): \[\label{eq:int95left} {I}_{0,x}^{(\mu)} u(x) := \frac{1}{\Gamma(\mu)} \int_0^x \frac{u(\xi)}{(x-\xi)^{1-\mu}} d\xi,\tag{5}\] \[\label{eq:int95right} {I}_{x,L}^{(\mu)} u(x) :=\frac{1}{\Gamma(\mu)} \int_x^L \frac{u(\xi)}{(\xi-x)^{1-\mu}} d\xi,\tag{6}\] for \(\mu>0\), and we set \({I}_{0,x}^{0}u = {I}_{x,L}^{0}u = u\), i.e. \({I}_{0,x}^{0} = {I}_{x,L}^{0} =I\), are the identity operators. The above defined operators are linear and bounded from \(L^\infty(0,L)\) into \(C(0,L)\) [16], [22], and respect the following properties:
Lemma 2. [10] If \(Re \{\alpha\} >0\) and \(u(x) \in L_p(0,L), \;1 \leq p \leq \infty\), then \[\begin{align} D^\alpha_{0,x} {I}_{x,0}^{(\alpha)} u(x)=u(x),\\ D^\alpha_{x,L} {I}_{x,L}^{(\alpha)} u(x)=u(x). \end{align}\]
Lemma 3. [10] Let \(Re \{\alpha\} >0\). If \(u(x) \in C^m[0,L]\), then the following relations hold: \[\begin{align} \label{eq:ILDL} I_{0,x}^\alpha\,D_{0,x}^\alpha u(x)=u(x)-\sum_{l=0}^{m-1} \frac{u^{(l)}(0)}{l!}x^l,\\ \label{eq:IRDR} I_{x,L}^\alpha\,D_{x,L}^\alpha u(x)=u(x)-\sum_{l=0}^{m-1} \frac{(-1)^lu^{(l)}(L)}{l!}(L-x)^l. \end{align}\] {#eq: sublabel=eq:eq:ILDL,eq:eq:IRDR}
Lemma 4. [22] The solution of the Carleman equation \[\label{eq:Carleman} \int_a^b \frac{\varphi(t)}{|x-t|^{1-\beta}}\, dt = g(x), \qquad a \le x \le b,\; 0<\beta<1,\qquad{(2)}\] is given by \[\begin{gather} \label{eq:sol95Carleman} \varphi(x) = \frac{\sin(\beta \pi)}{\tilde{A}\pi} \frac{d}{dx} \int_a^x \frac{g(t) dt}{(x-t)^\beta} \\ - \frac{1}{\tilde{A}} \left( \frac{\sin(\beta \pi)}{\pi}\right)^2 \frac{d}{dx} \int_a^x \frac{Z(t) dt}{(x-t)^\beta} dt \frac{d}{dt} \int_a^t \frac{d\tau}{(t-\tau)^{1-\beta}} \int_\tau^b \frac{g(s)ds}{Z(s) (s-\tau)^\beta}, \end{gather}\qquad{(3)}\] where \[\tilde{A}=2+2\cos(\beta \pi)=4\cos^2(\frac{\beta \pi}{2}),\]
\[Z(t)=\left( \frac{b-t}{t-a} \right)^{\beta/2},\] and \(g(x)\in H_\beta^*\), where \(H_\beta^*\) denotes the weighted Hölder space of order more than \(\beta\).
We will now proceed following a similar approach to [9], and extending such analysis to the case of the Riesz-Caputo operator with the boundary conditions specified in 1 . In order to do that, we will transform the original fractional differential equation into an integral equation, and use the Schauder’s and the Banach’s fixed point theorems to show that the solution exists and is unique.
Lemma 5. Let \(1<\alpha<2\), \(a_i,b_i,c_i \in \mathbb{R}\) with \(a_i+b_i \neq 0\) for \(i=0,1\), and assume \(f:[0,L]\times \mathbb{R}\rightarrow \mathbb{R}\) and \(f(x,u(x))\in H_{2-\alpha}^*.\) Then \(u \in
C^1[0,L]\) is a solution of 1 if and only if it satisfies the integral equation \[\begin{align}
\label{eq:reform12}
u(x)
&= \frac{A}{\alpha - 1}
\int_0^x f(s,u(s))\,(x - s)^{\alpha - 1}\,ds \notag \\
&+ \frac{B}{\alpha - 1}
\int_0^x h'(\tau)
\left(\frac{L - \tau}{\tau}\right)^{\frac{2 - \alpha}{2}}
(x - \tau)^{\alpha - 1}\,d\tau
+ C_1 x + C_0
\end{align}\qquad{(4)}\] where \(A = - \frac{1}{\pi} \tan \left( \frac{\alpha \pi}{2} \right)\), \(B=2 \frac{\sin^2 \left(\frac{\alpha \pi}{2} \right)}{\pi^2}\), \[h(t)= \int_0^t \frac{1}{(t-r)^{\alpha-1}} F(r) dr,\] with \(F(r)=\int_r^L \left( \frac{L-s}{s} \right)^{-(2-\alpha)/2} \frac{f(s,u(s))}{(s-r)^{2-\alpha}} ds\).
The constants \(C_0\) and \(C_1\) are determined by the boundary data and are given by: \[\label{eq:A}
C_1
= \frac{c_1 - b_1\left(A J_1 + B J_2\right)}{a_1 + b_1},\qquad{(5)}\] and \[\label{eq:B}
C_0
= \frac{
c_0
- b_0 \left[
\frac{A}{\alpha - 1} I_1
+ \frac{B}{\alpha - 1} I_2
\right]
- b_0 L C_1
}{
a_0 + b_0
}.\qquad{(6)}\] with \[I_1=\int_0^L f(s,u(s)) (L-s)^{\alpha-1}ds,\] \[I_2=\int_0^L h'(\tau)(L-\tau)^{\alpha/2}(\tau)^{-(2-\alpha)/2} d\tau,\] \[J_1=\int_0^L f(s,u(s)) (L-s)^{\alpha-2}ds,\] \[J_2=\int_0^L h'(\tau)(L-\tau)^{(\alpha-2)/2}(\tau)^{-(2-\alpha)/2} d\tau.\]
Proof. Without loss of generality, throughout the proof we consider a generic interval \(I=[a,b]\). We will proceed by rewriting eq. 1 in the form ?? for which we have an explicit
solution, and then compute the corresponding constants that satisfy the boundary conditions.
Since \(1<\alpha<2\), we have \(m=2\) and \(\beta=2-\alpha\).
The Carleman equation ?? can be written as \[\label{eq:rewrite95carleman}
\Gamma(\beta)\, I^{(\beta)}_{a,b}\varphi(x) = g(x),
\qquad \beta\in(0,1),\tag{7}\]
The Caputo fractional derivatives can also be written as [10]: \[D^\alpha_{a,x}u(x) =\frac{1}{\Gamma(2-\alpha)}\int_a^x \frac{u''(t)}{(x-t)^{\alpha-1}}\,dt = \bigl(I^{(2-\alpha)}_{a,x}D^2u\bigr)(x),\] \[D^\alpha_{x,b}u(x) =\frac{1}{\Gamma(2-\alpha)}\int_x^b \frac{u''(t)}{(t-x)^{\alpha-1}}\,dt = \bigl(I^{(2-\alpha)}_{x,b}D^2u\bigr)(x).\]
The Rietz–Caputo derivative is therefore \[D^\alpha_{\mathrm{RC}}u(x) =\frac{1}{2}\Bigl(D^\alpha_{a,x}u(x)+D^\alpha_{x,b}u(x)\Bigr) =\frac{1}{2}\bigl(I^{(2-\alpha)}_{a,b}D^2u\bigr)(x).\]
Multiplying both sides of 1 by \(2\Gamma(2-\alpha)\) yields \[\Gamma(2-\alpha)\, I^{(2-\alpha)}_{a,b}D^2u(x) =2f(x,u(x)).\]
Now if we set \[\varphi(t)=u''(t), \qquad g(x)=2f(x,u(x)),\] we recognize the Carleman equation in the form 7 and we can write the unique solution by making use of Lemma 4. For semplicity, we split the solution in two parts4: \[\varphi(x)=A\,\varphi_1(x)+B\,\varphi_2(x),\] where \[\varphi_1(x) =\frac{d}{dx}\int_a^x \frac{f(t,u(t))}{(x-t)^{2-\alpha}}\,dt,\] \[\varphi_2(x) =\frac{d}{dx}\int_a^x \left(\frac{b-t}{t-a}\right)^{\frac{2-\alpha}{2}} \frac{h'(t)}{(x-t)^{2-\alpha}}\,dt,\] with \[A=-\frac{1}{\pi}\tan\!\left(\frac{\alpha\pi}{2}\right), \qquad B=2\frac{\sin^2(\alpha\pi/2)}{\pi^2}.\] and \[h(t)=\int_a^t \frac{1}{(t-r)^{\alpha-1}}\,F(r)\,dr, \qquad F(r)=\int_r^b \left(\frac{b-s}{s-a}\right)^{-\frac{2-\alpha}{2}} \frac{f(s,u(s))}{(s-r)^{2-\alpha}}\,ds.\]
Now, since \(\varphi(x)=u''(x)\), the solution of the fractional differential problem is obtained by integrating twice.
\[\begin{align} u'(x)&=\int_a^x \varphi(t)\,dt + C_1 \nonumber \\ &= A \int_a^x \frac{f(s,u(s))}{(x-s)^{2-\alpha}}\,ds + B \int_a^x \frac{ \left(\dfrac{b-\tau}{\tau-a}\right)^{\frac{2-\alpha}{2}} }{ (x-\tau)^{2-\alpha} } \,h'(\tau)\,d\tau + C_1, \end{align}\] with \(C_1=u'(a)\).
\[\begin{align} u(x)&=\int_a^x u'(t)\,dt + C_0 \nonumber \\ &=\frac{A}{\alpha-1}\int_a^x f(s,u(s))(x-s)^{\alpha-1}\,ds \nonumber \\ &+\frac{B}{\alpha-1}\int_a^x h'(s)\left(\frac{b-s}{s-a}\right)^{\frac{2-\alpha}{2}} (x-s)^{\alpha-1}\,ds +C_1(x-a)+C_0, \end{align}\] with \(C_0=u(a)\).
The boundary conditions given by eq. 1 are: \[a_0u(a)+b_0u(b)=c_0,
\qquad
a_1u'(a)+b_1u'(b)=c_1.\] where, by hypothesis \(a_0+b_0 \neq 0\) and \(a_1+b_1 \neq 0\).
Evaluating the solution and its derivative at \(x=a\) and \(x=b\), one obtains a linear system for \(C_0\) and \(C_1\) that
gives: \[C_1=\frac{c_1-b_1(AJ_1+BJ_2)}{a_1+b_1},\] \[C_0=\frac{c_0-b_0\bigl[\frac{A}{\alpha-1}I_1+\frac{B}{\alpha-1}I_2\bigr]
-b_0(b-a)C_1}{a_0+b_0},\] where \[I_1=\int_a^b f(s,u(s))(b-s)^{\alpha-1}\,ds,
\quad
I_2=\int_a^b h'(s)(b-s)^{\alpha/2}(s-a)^{-(2-\alpha)/2}\,ds,\] \[J_1=\int_a^b \frac{f(s,u(s))}{(b-s)^{2-\alpha}}\,ds,
\quad
J_2=\int_a^b
h'(s)(b-s)^{(\alpha-2)/2}(s-a)^{-(2-\alpha)/2}\,ds.\] By choosing \(a=0\) and \(b=L\), all the claims follow. ◻
Theorem 1. Assume the hypotheses of Lemma 5 and suppose that \(f:[0,L]\times\mathbb{R}\to\mathbb{R}\) is Lipschitz continuous with respect to the second variable, i.e. \[|f(x,y_1)-f(x,y_2)|\le L_f |y_1-y_2|, \qquad \forall x\in[0,L],\; y_1,y_2\in\mathbb{R}.\] If \[L_f < \frac{1}{K},\] where \[\begin{align} K &:= \frac{1}{\alpha-1} \Bigl( |A|L^{\alpha-1} + |B|L^{\alpha-1} \Bigr)\\ &+ \frac{|b_1|}{|a_1+b_1|} \Bigl( |A|J_1^\ast+|B|J_2^\ast \Bigr)\\ &+ \frac{|b_0|}{|a_0+b_0|} \Bigl( |A|I_1^\ast+|B|I_2^\ast+L \Bigr), \end{align}\] and \(I_i^\ast,J_i^\ast\) are suitable positive constants defined below, then problem 1 admits a unique solution \(u\in C^1[0,L]\).
Proof. By Lemma 5, \(u\in C^1[0,L]\) is a solution of problem 1 if and only if it satisfies the integral equation \[u(x) = (\mathcal{T}u)(x),\] where \[\begin{align} (\mathcal{T}u)(x) &= \frac{A}{\alpha - 1} \int_0^x g_u(s)\,(x - s)^{\alpha - 1}\,ds \notag \\ &\quad + \frac{B}{\alpha - 1} \int_0^x \Phi_u(\tau)\,(x - \tau)^{\alpha - 1}\,d\tau + C_1(u) x+C_0(u), \end{align}\] and \[\Phi_u(\tau) = h_u'(\tau) \left(\frac{L - \tau}{\tau}\right)^{\frac{2 - \alpha}{2}}, \qquad g_u(x)=f(x,u(x)).\]
The constants \(C_0(u)\) and \(C_1(u)\) depend linearly on \(g_u\) and are given by \[C_1(u) = \frac{c_1 - b_1\left(A J_1(g_u) + B J_2(g_u)\right)}{a_1 + b_1},\] \[C_0(u) = \frac{ c_0 - b_0 \left[ \frac{A}{\alpha - 1} I_1(g_u) + \frac{B}{\alpha - 1} I_2(g_u) \right] - b_0 L C_1(u) }{ a_0 + b_0 }.\]
Let \(u,v\in C^1[0,L]\). Then \[(\mathcal{T}u)(x)-(\mathcal{T}v)(x) = T_1(x)+T_2(x)+T_3(x),\] where \[\begin{align} T_1(x) &= \frac{A}{\alpha-1} \int_0^x \bigl[g_u(s)-g_v(s)\bigr](x-s)^{\alpha-1}\,ds, \\ T_2(x) &= \frac{B}{\alpha-1} \int_0^x \bigl[\Phi_u(\tau)-\Phi_v(\tau)\bigr](x-\tau)^{\alpha-1}\,d\tau, \\ T_3(x) &= \bigl(C_1(u)-C_1(v)\bigr)x + \bigl(C_0(u)-C_0(v)\bigr). \end{align}\]
By the Lipschitz continuity of \(f\), \[|g_u(s)-g_v(s)| \le L_f |u(s)-v(s)| \le L_f \|u-v\|_\infty.\] Hence \[|T_1(x)| \le \frac{|A|}{\alpha-1} L_f \|u-v\|_\infty \int_0^x (x-s)^{\alpha-1}\,ds \le \frac{|A|}{\alpha-1} L_f \|u-v\|_\infty \frac{L^{\alpha}}{\alpha}.\]
Since the mapping \(g\mapsto \Phi\) defines a linear and bounded operator from \(L^\infty(a,b)\) into itself, there exists a constant \(C>0\) such that \[\|\Phi_u-\Phi_v\|_\infty \le C\,\|g_u-g_v\|_\infty.\] Thus \[|T_2(x)| \le \frac{|B|}{\alpha-1} C\,L_f \|u-v\|_\infty \int_0^x (x-\tau)^{\alpha-1}\,d\tau \le \frac{|B|}{\alpha-1} L_f \|u-v\|_\infty \frac{L^{\alpha}}{\alpha},\] where the constant \(C\) has been absorbed into \(|B|\).
Using the linear dependence of \(C_1\) and \(C_0\) on \(g_u\), we obtain \[|C_1(u)-C_1(v)| \le \frac{|b_1|}{|a_1+b_1|} \Bigl( |A|\,|J_1(g_u-g_v)| + |B|\,|J_2(g_u-g_v)| \Bigr),\] and similarly \[|C_0(u)-C_0(v)| \le \frac{|b_0|}{|a_0+b_0|} \Bigl( |A|\,|I_1(g_u-g_v)| + |B|\,|I_2(g_u-g_v)| + L\,|C_1(u)-C_1(v)| \Bigr).\]
By boundedness of the integral operators \(I_i\) and \(J_i\) in \(L^\infty(0,L)\), there exist constants \(I_i^\ast,J_i^\ast>0\) such that \[|I_i(g_u-g_v)| \le I_i^\ast L_f\|u-v\|_\infty, \qquad |J_i(g_u-g_v)| \le J_i^\ast L_f\|u-v\|_\infty.\]
Therefore \[|T_3(x)| \le L_f K_3 \|u-v\|_\infty,\] for a suitable constant \(K_3>0\).
Collecting the estimates for \(T_1\), \(T_2\) and \(T_3\), we obtain \[\|\mathcal{T}u-\mathcal{T}v\|_\infty \le L_f K \|u-v\|_\infty,\] where \(K\) is defined in the statement of the theorem. If \(L_f<1/K\), then \(\mathcal{T}\) is a contraction on \(C^1[0,L]\). Hence, by the Banach fixed point theorem, \(\mathcal{T}\) admits a unique fixed point, which is the unique solution of problem 1 . ◻
A spline is a piecewise polynomial function widely used in computational mathematics for tasks such as curve fitting, data smoothing, and numerical approximation of integral and differential equations. On the real line splines can be written as linear
combination of the translate of cardinal B-splines, which are compactly supported piecewise polynomials having breakpoints at the integers. The cardinal B-spline \(N_n\) of integer degree \(n>0\) can be defined through the truncated power function \(x^n_+ := \max(0,x)^n\) as \[N_{n}(x) = (n+1) [0,1,\ldots,n+1] (y-x)^n_+,\] where \([x_0,x_1,\ldots,x_{n+1}] f\) denotes the divided difference of the function \(f\) on the sequence of knots \(\{x_0,x_1,\ldots,x_{n+1}\}\). This definition works
also for the case \(n=0\), but we will consider only the case of (at least) continuous polynomials.
On a finite interval \([0,L]\), \(L \in \mathbb{N}\), splines can be represented as linear combination of the optimal basis, which have integer knots of multiplicity one inside the interval
and knots of multiplicity \(n+1\) at the endpoints. Optimality refers to a set of properties, such as compact support, high smoothness, endpoint conditions, partition of unity and numerical stability, that make these bases
ideal for approximating functions. For more details on splines see the basic monographs [23], [24].
The optimal basis \({\cal N}_n = \left\{N_{j,n}(x), 0\le j \le L+n-1 \right\}\) on the interval \([0,L]\) has \(L+n\) basis functions. The \(2n\) functions \(N_{j,n}\) and \(N_{L+n-1-j,n}\), \(0 \le j \le n-1\), are respectively the left and right boundary functions. The \(L-n\) functions \(N_{j,n}\), \(n \le j \le L-1\), are the internal functions (see [25] for more details).
We observe that the optimal basis is centrally symmetric and fulfills the endpoint conditions \[\label{eq:endpoint} N_{j,n}(0) = \delta_{j0}\,, \quad N_{L+n-1-j,n}(L) = \delta_{j0}\,, \qquad 0 \le j \le L+n-1,\tag{8}\] where \(\delta_{j0}\) denotes the Kronecker symbol. The optimal basis \({\cal N}_n\) can be refined by dilation [23], [24]. Let \(h\) be the refinement step, that is, the distance between the refined knot sequence: \[\begin{cases} x_{0,h} = x_{1,h}=\cdots =x_{n-1,h} = 0; \\ x_{j,h} = (j-n+1)h\,, \text{ for } n\le j \le M_h; \\ x_{M_h+1,h}=x_{M_h+2,h}=\cdots=x_{M_h+n-1,h}=L\,, \end{cases}\] where \(M_h = L/h+n\). The refined basis \({\cal N}_{h,n} = \{N_{j,h,n}(x), \;0 \le j \le M_h\}\), \(x \in [0,L]\), has \(L/h-n\) internal functions and \(2n\) boundary functions. We notice that refining the knots only leads to an increase in the number of internal functions, while the number of boundary functions remains unchanged.
We will use these basis functions to represent the approximate solution of the fractional differential problem 1 . Another great advantage of using such a representation is that fractional derivatives of the optimal B-spline basis functions can be computed analytically using the explicit formula provided in [7], [25] and by leveraging the symmetry properties of the Riesz-Caputo operator and the basis functions together. Specifically,
Proposition 2 ([7]). Let \(\Phi = \{\phi_j(x),0\le j \le M\}\) be a centrally symmetric basis function in the interval \([0,L]\), i.e., \(\phi_j(x) = \phi_{M-j}(L-x),\) for \(0 \le j \le M\) and \(x \in [0,L]\). Then, the left and right Caputo derivatives of the basis function \(\Phi\) satisfy the following symmetry properties: \[\label{eq:symm95derphi} D^\alpha_{x,L}\phi_{j}(x)=D^\alpha_{0,L-x} \phi_{j}(L-x), \qquad 0 \le j \le M.\qquad{(7)}\]
In [8] we discussed the advantages of using Greville Abscissae as collocation points, and compared this choice with other possible distributions of points. Motivated by the effectiveness of choosing Greville Abscissae, also shown by numerical experiments, in this paper we will adopt the same option. Greville Abscissae are defined as [23], [24], [26] : \[\label{eq:Grevpoint} \eta_{p,h,n}=\frac{x_{p,h}+...+x_{p+n-1,h}}{n}, \qquad 0\le p \le M_h.\tag{9}\] We solve the fractional boundary value problem 1 by the spline collocation method introduced in [27]. First, we approximate the solution to the differential problem by a spline function: \[\label{eq:spline95approx} u(x) \approx u_{h,n}(x) = \sum_{k=0}^{M_h} \ell_{k,h,n} N_{k,h,n}(x), \qquad x \in [0,L]\,.\tag{10}\] Then, by collocating the differential problem 1 at the Greville Abscissae 9 , we obtain the square nonlinear system \[\label{eq:coll95eq} \left \{ \begin{array}{ll} \displaystyle \sum_{k=0}^{M_h} \ell_{k,h,n} D_{RC}^\alpha N_{k,h,n}(\eta_{p,h,n})=f\bigl(\eta_{p,h,n}, \sum_{k=0}^{M_h} \ell_{k,h,n} N_{k,h,n}(\eta_{p,h,n})\bigr), \;\;\;0<p<M_h, \\ \\ \displaystyle a_i \sum_{k=0}^{M_h} \ell_{k,h,n} N_{k,h,n}^{(i)}(\eta_{0,h,n})+b_i\sum_{k=0}^{M_h} \ell_{k,h,n} N_{k,h,n}^{(i)}(\eta_{M_h,h,n})=c_i, \;\;i=0,..,{\lfloor \alpha \rfloor}, \end{array} \right.\tag{11}\] where the coefficients \(\{\ell_{k,h,n}\}\) are the \(M_h+1\) unknowns.
In this section we aim to prove the convergence for the proposed method.
Theorem 3. Assume the hypotheses of Theorem 1, i.e., that \(f:[0,L]\times \mathbb{R}\to \mathbb{R}\) is continuous and Lipschitz in the second variable with constant \(L_f\) sufficiently small. Let \(u(x)\) be the exact solution of problem 1 and \(u_{h}(x)\) its spline collocation approximation defined by 11 . Then \[\label{eq:conv95norm} \lim_{h\to 0} \| u - u_h \|_\infty = 0.\qquad{(8)}\]
Proof. For simplicity, we omit the subscript \(n\) in the spline basis.
Let \(\bar{\eta}_h\) denote the Greville abscissa nearest to \(x\). Then \[\label{eq:error95split} \| u(x) - u_h(x) \|_\infty \le \underbrace{\| u(x) - u_h(\bar{\eta}_h) \|_\infty}_{\text{discrete error}} + \underbrace{\| u_h(\bar{\eta}_h) - u_h(x) \|_\infty}_{\text{interpolation error}}.\tag{12}\]
By continuity of \(u_h\), the second term vanishes as \(h \to 0\): \[\lim_{h\to 0} \| u_h(\bar{\eta}_h) - u_h(x) \|_\infty = 0.\]
By Lemma 5, both \(u\) and \(u_h\) satisfy the integral equation \[\begin{align} u(x) &= \frac{A}{\alpha-1} \int_0^x f(t,u(t)) (x-t)^{\alpha-1} dt \\ & + \frac{B}{\alpha-1} \int_0^x h'(t)\left(\frac{L - \tau}{\tau}\right)^{\frac{2 - \alpha}{2}} (x-t)^{\alpha-1} dt + C_1 x + C_0, \end{align}\] with analogous expression for \(u_h\) (with \(f(t,u_h(t))\) and \(h'(t,u_h)\)). Subtracting the two equations and taking absolute values: \[\begin{align} |u(x)-u_h(\bar{\eta}_h)| &\le \frac{|A|}{\alpha-1} \int_0^x |f(t,u(t)) - f(t,u_h(t))| (x-t)^{\alpha-1} dt \\ &\quad + \frac{|B|}{\alpha-1} \int_0^x |h'(t,u) - h'(t,u_h)|\left(\frac{L - \tau}{\tau}\right)^{\frac{2 - \alpha}{2}} (x-t)^{\alpha-1} dt \\ &\quad + |C_1(u)-C_1(u_h)| \, |x| + |C_0(u)-C_0(u_h)|. \end{align}\] By the Lipschitz continuity of \(f\), \[\| f(\cdot,u) - f(\cdot,u_h) \|_\infty \le L_f \| u - u_h \|_\infty.\]
Moreover, the mappings \(g\mapsto h'\) and the linear combinations appearing in \(C_0,C_1\) are linear and bounded operators on \(L^\infty\) (see discussion in Section 2.1). Hence, there exists a constant \(K>0\), independent of \(h\), such that \[|u(x)-u_h(\bar{\eta}_h)| \le K L_f \| u - u_h \|_\infty + K_2 |x-\bar{\eta}_h|.\] Taking supremum over \(x \in [0,L]\): \[\| u - u_h(\bar{\eta}_h) \|_\infty \le K L_f \| u - u_h \|_\infty + K_2 \max_x |x-\bar{\eta}_h|.\]
By Theorem 1, \(L_f\) is sufficiently small so that \(K L_f < 1\). Then \[(1 - K L_f) \| u - u_h(\bar{\eta}_h) \|_\infty \le K_2 \max_x |x-\bar{\eta}_h|.\] As \(h\to 0\), \(\max_x |x-\bar{\eta}_h| \to 0\) (Greville points become dense), so the right-hand side vanishes. Together with Step 1, we conclude \[\lim_{h\to 0} \| u - u_h \|_\infty = 0.\] ◻
In this section, we validate the theoretical results through numerical experiments. Without loss of generality, we set \(L=1\) . For the numerical tests, we utilize quadratic and cubic optimal B-splines, which offer
flexibility in representing the solution with different levels of accuracy and computational cost. The nonlinear square system (11 ) is solved using MATLAB’s solve function, which is specifically designed for
efficiently handling nonlinear systems.
In the first example, we solve the nonlinear fractional boundary value problem
\[\left \{ \begin{array}{lll} D_{RC}^\alpha u(x)=u(x)^2+g_1(x), \;\;\;x \in [0,1], \\ \\ u(0)+u(1)=1, \\ \\ u'(0)+u'(1)=2.5, \end{array} \right.\] where \(g_1(x)\) is given in 7. The analytical solution of this problem is \(u(x)=x^2\sqrt{x}\). The values of the infinity norm of the error \(E_{h,n}\) for \(n=2\) and \(n=3\) are listed in Table ¿tbl:table1951? and Table ¿tbl:table1952?, respectively.
1.5cm0cm
| \(h / \alpha\) | 1.25 | 1.5 | 1.75 | 1.95 |
|---|---|---|---|---|
| \(\frac 1 { 8}\) | \(3.98e-03\) | \(2.93e-03\) | \(3.33e-03\) | \(3.95e-03\) |
| \(\frac 1 {16}\) | \(1.09e-03\) | \(8.00e-04\) | \(9.93e-04\) | \(1.25e-03\) |
| \(\frac 1 {32}\) | \(2.93e-04\) | \(2.20e-04\) | \(3.01e-04\) | \(4.01e-04\) |
| \(\frac 1 {64}\) | \(7.87e-05\) | \(6.12e-05\) | \(9.25e-05\) | \(1.31e-04\) |
1.5cm0cm
| \(h / \alpha\) | 1.25 | 1.5 | 1.75 | 1.95 |
|---|---|---|---|---|
| \(\frac 1 { 8}\) | \(6.33e-04\) | \(4.28e-04\) | \(4.34e-04\) | \(5.93e-04\) |
| \(\frac 1 {16}\) | \(1.60e-04\) | \(1.06e-04\) | \(9.23e-05\) | \(1.65e-04\) |
| \(\frac 1 {32}\) | \(4.08e-05\) | \(2.69e-05\) | \(1.84e-05\) | \(5.46e-05\) |
| \(\frac 1 {64}\) | \(1.06e-05\) | \(7.00e-06\) | \(4.12e-06\) | \(2.05e-05\) |
In this example, we approximate the solution of the following fractional problem \[\label{eq:ex2} \left \{ \begin{array}{lll} D_{RC}^\alpha u(x)=\sin \bigl(\frac{\pi}{4} u(x) \bigr)+g_2(x), \;\;\;x \in [0,1], \\ \\ u(0)+u(1)= \frac{\sqrt{2}}{2}, \\ \\ u'(0)+u'(1)=\frac{ \pi}{4}\left(1+\frac{\sqrt{2}}{2} \right), \end{array} \right.\tag{13}\] where \(g_2(x)\) is given in 7. In this case the analytical solution is the periodic function \(u(x)=\sin\bigl(\frac{\pi}{4} x \bigl)\). The values of the infinity norm of the error \(E_{h,n}\) are reported in Tables ¿tbl:table2951?-¿tbl:table2952?.
1.5cm0cm
| \(h / \alpha\) | 1.25 | 1.5 | 1.75 | 1.95 |
|---|---|---|---|---|
| \(\frac 1 { 8}\) | \(3.47e-04\) | \(2.14e-04\) | \(2.45e-04\) | \(2.87e-04\) |
| \(\frac 1 {16}\) | \(7.06e-05\) | \(4.54e-05\) | \(5.64e-05\) | \(7.06e-05\) |
| \(\frac 1 {32}\) | \(1.40e-05\) | \(9.58e-06\) | \(1.29e-05\) | \(1.73e-05\) |
| \(\frac 1 {64}\) | \(2.73e-06\) | \(2.01e-06\) | \(2.96e-06\) | \(4.26e-06\) |
1.5cm0cm
| \(h / \alpha\) | 1.25 | 1.5 | 1.75 | 1.95 |
|---|---|---|---|---|
| \(\frac 1 { 8}\) | \(2.84e-06\) | \(4.60e-06\) | \(1.05e-05\) | \(2.18e-05\) |
| \(\frac 1 {16}\) | \(5.81e-07\) | \(8.67e-07\) | \(2.37e-06\) | \(5.69e-06\) |
| \(\frac 1 {32}\) | \(1.09e-07\) | \(1.64e-07\) | \(5.20e-07\) | \(1.43e-06\) |
| \(\frac 1 {64}\) | \(1.89e-08\) | \(3.04e-08\) | \(1.12e-07\) | \(3.52e-07\) |
In all the presented numerical example we have verified the convergence theoretically proved in section 3.
In this paper we study a nonlinear fractional boundary value problem involving the Riesz-Caputo operator. We prove existence and uniqueness of the solution by showing the equivalence with a Fredholm integral equation, and then using the Schauder’s and the Banach fixed point theorems. We also present a collocation method that project the solution in a spline space. The approximate solution is represented as a linear combination of optimal B-spline bases, whose coefficients are determined by solving a nonlinear square system. The convergence analysis of such method is provided with theoretical insights. The numerical experiments confirm the theoretical findings. As outlined in Section 1, the available literature involving the Riesz-Caputo operator is very limited, especially when dealing with nonlinear problems. In this work we have rigorously proved the convergence of the numerical method. A detailed analysis of the order of convergence is still missing and will be the subject of future studies.
We derive explicitly the formulas for \(C_0\) and \(C_1\) appearing in Lemma 5. Recall the representation
\[\begin{align}
\label{eq:A1} u(x)&=\frac{A}{\alpha-1}\int_a^x f(s,u(s))(x-s)^{\alpha-1}\,ds \nonumber \\ &+\frac{B}{\alpha-1}\int_a^x h'(s)\left(\frac{b-s}{s-a}\right)^{\frac{2-\alpha}{2}} (x-s)^{\alpha-1}\,ds +C_1(x-a)+C_0,
\end{align}\tag{14}\] and \[\begin{align}
\label{eq:A2} u'(x)&= A \int_a^x \frac{f(s,u(s))}{(x-s)^{2-\alpha}}\,ds + B \int_a^x \frac{ \left(\dfrac{b-\tau}{\tau-a}\right)^{\frac{2-\alpha}{2}} }{ (x-\tau)^{2-\alpha} } \,h'(\tau)\,d\tau + C_1,
\end{align}\tag{15}\] where \(A = - \frac{1}{\pi} \tan \left( \frac{\alpha \pi}{2} \right)\), \(B= 2\frac{\sin^2 \left(\frac{\alpha \pi}{2} \right)}{\pi^2}\), \[h(t)= \int_a^t \frac{1}{(t-r)^{\alpha-1}} F(r) dr.\] with \(F(r)=\int_r^b \left( \frac{b-s}{s-a} \right)^{-(2-\alpha)/2} \frac{f(s,u(s))}{(s-r)^{2-\alpha}} ds\).
\[\begin{align} u(a)&=C_0, \end{align}\]
\[\begin{align} u(b)&=\frac{A}{\alpha-1}\int_a^b f(s,u(s))(b-s)^{\alpha-1}\,ds \nonumber \\ &+\frac{B}{\alpha-1}\int_a^b h'(s)\left(\frac{b-s}{s-a}\right)^{\frac{2-\alpha}{2}} (b-s)^{\alpha-1}\,ds +C_1(b-a)+C_0 \\ &=\frac{A}{\alpha-1} I_1+\frac{B}{\alpha-1} I_2 +C_1(b-a)+C_0 \end{align}\] setting \[\begin{align} I_{1}&=\int_{a}^{b}f(s,u(s))(b-s)^{\alpha-1}\,ds, \tag{16} \\ I_{2}&=\int_{a}^{b}h'(\tau)(b-\tau)^{\alpha/2}(\tau-a)^{-(2-\alpha)/2}\,d\tau. \tag{17} \end{align}\]
We use 15 . First,
\[\begin{align} u'(a)&= C_1, \end{align}\]
\[\begin{align} u'(b)&= A \int_a^b \frac{f(s,u(s))}{(b-s)^{2-\alpha}}\,ds + B \int_a^b \frac{ \left(\dfrac{b-\tau}{\tau-a}\right)^{\frac{2-\alpha}{2}} }{ (b-\tau)^{2-\alpha} } \,h'(\tau)\,d\tau + C_1, \end{align}\] and denoting
\[\begin{align} J_{1}&=\int_{a}^{b}\frac{f(s,u(s))}{(b-s)^{2-\alpha}}\,ds, \\ J_{2}&=\int_{a}^{b}h'(\tau)(b-\tau)^{(\alpha-2)/2}(\tau-a)^{-(2-\alpha)/2}\,d\tau. \end{align}\] We have \[\begin{align} u'(b)&=- A J_1+ B J_2+ C_1, \end{align}\]
The boundary conditions are \[\begin{align} a_0 u(a) + b_0 u(b) &= c_0, \\ a_1 u'(a) + b_1 u'(b) &= c_1. \end{align}\] Substituting the expressions above gives \[\begin{align} a_{0}C_{0}+b_{0}\Bigl[\frac{A}{\alpha-1}I_{1}+\frac{B}{\alpha-1}I_{2}+C_{1}(b-a)+C_{0}\Bigr]&=c_{0},\\ a_{1}C_{1}+b_{1}\bigl[A J_{1}+B J_{2}+C_{1}\bigr]&=c_{1} \end{align}\] hence \[\begin{align} (a_{0}+b_{0})C_{0}+b_{0}(b-a)C_{1}+b_{0}\Bigl[\frac{A}{\alpha-1}I_{1}+\frac{B}{\alpha-1}I_{2}\Bigr]&=c_{0}, \tag{18}\\ (a_{1}+b_{1})C_{1}+b_{1}\bigl[A J_{1}+B J_{2}\bigr]&=c_{1}. \tag{19} \end{align}\] Equation 19 becomes \[C_{1}= \frac{c_{1}-b_{1}\bigl(A J_{1}+B J_{2}\bigr)}{a_{1}+b_{1}} .\] Substituting into 18 \[C_{0}= \frac{c_{0}-b_{0}\Bigl[\frac{A}{\alpha-1}I_{1}+\frac{B}{\alpha-1}I_{2}\Bigr]-b_{0}(b-a)C_{1}}{a_{0}+b_{0}} .\] Setting \(a=0\) and \(b=L\), we obtain exactly formulas ?? –?? of Lemma 5.
Here we list the analytical expressions of the known terms \(g_i(x)\), for \(i=1,2\) used in the numerical tests of Section 4.
In the following we denote by \(_pF_q(a;b;z)\) the generalized hypergeometric function.
In the first example, the function \(g_1(x)\) is given by \[\begin{align} g_1(x)& =-x^{5}+\frac{1}{2} \left[\frac{\Gamma(3.5)}{\Gamma(3.5-\alpha)}- \frac{3.75}{2\sqrt\pi} \Gamma(\alpha-2.5)\right] x^{2.5- \alpha}\\ &-\frac{1}{2}\frac{3.75}{\Gamma(2- \alpha)}\frac{L^{2.5-\alpha}}{(\alpha-2.5)} \, _2F_1\left(\alpha-2.5,\alpha-1;\alpha-1.5;\frac{x}{L}\right)\text{for } 1 < \alpha < 2. \end{align}\]
In the second example, we use \[g_2(x) =-\sin(\sin(\omega x))-\frac{\omega^2}{2\Gamma(2-\alpha)} \left(g_{11}(x)+g_{21}(x)\right)\] where \[\begin{align} g_{11}(x)&=\frac{\omega x^{3-\alpha} \, _1F_2\left(1;2-\frac{\alpha}{2},\frac{5}{2}-\frac{\alpha}{2};-\frac{1}{4} \omega^2 x^2\right)}{(\alpha-3)(\alpha-2)},\\ g_{21}(x)&= (L-x)^{2-\alpha} \frac{\omega (x-L) \cos (\omega x) \, _1F_2\left(\frac{3}{2}-\frac{\alpha}{2};\frac{3}{2},\frac{5}{2}-\frac{\alpha}{2};-\frac{1}{4} \omega^2 (L-x)^2\right)}{\alpha-3}\\ &-(L-x)^{2-\alpha}\frac{\sin (\omega x) \, _1F_2\left(1-\frac{\alpha}{2};\frac{1}{2},2-\frac{\alpha}{2};-\frac{1}{4} \omega^2 (L-x)^2\right)}{\alpha-2}, \end{align}\] and \(\omega=\frac{ \pi}{4}\).
The authors are members of the INdAM Research Group GNCS.
The authors declare that they have no conflict of interest.
Department of Basic and Applied Sciences for Engineering, Sapienza University of Rome, Italy. Email: chiara.sorgentone@uniroma1.it↩︎
Department of Industrial and Information Engineering and Economics, University of L’Aquila, Italy. Email: enza.pellegrino@univaq.it↩︎
Department of Basic and Applied Sciences for Engineering, Sapienza University of Rome, Italy. Email: francesca.pitolli@uniroma1.it↩︎
In this step we have made use of the trigonometric identities: \[\sin((2-\alpha)\pi)=-\sin(\alpha \pi), \;\; \sin(\alpha \pi)=2\sin(\frac{\alpha \pi}{2})\cos(\frac{\alpha \pi}{2}), \;\;\cos((2-\alpha)\pi/2)=-\cos(\frac{\alpha \pi}{2}).\]↩︎
Again we split the integration in two parts: \[I_1(x) = \int_a^x \varphi_1(t)\,dt = \left[ \int_a^t \frac{f(s,u(s))}{(t-s)^{2-\alpha}}\,ds \right]_{t=a}^{t=x} = \int_a^x \frac{f(s,u(s))}{(x-s)^{2-\alpha}}\,ds .\] and \[\int_a^x \varphi_2(t)\,dt = \left[ \int_a^t \frac{ \left(\dfrac{b-\tau}{\tau-a}\right)^{\frac{2-\alpha}{2}} }{ (t-\tau)^{2-\alpha} } \,h'(\tau)\,d\tau \right]_{t=a}^{t=x} = \int_a^x \frac{ \left(\dfrac{b-\tau}{\tau-a}\right)^{\frac{2-\alpha}{2}} }{ (x-\tau)^{2-\alpha} } \,h'(\tau)\,d\tau .\]↩︎
We start by computing the first double integral: \[\tilde{I}_1 = A \int_a^x \left[ \int_a^t \frac{f(s,u(s))}{(t-s)^{2-\alpha}} \,ds \right] dt .\] Now if change the integration order we get: \[\begin{align} \tilde{I}_1 &= A \int_{s=a}^{x} f(s,u(s)) \left[ \int_{t=s}^{x} (t-s)^{\alpha-2} \, dt \right] ds =\frac{A}{\alpha - 1} \int_{a}^{x} f(s,u(s)) \, (x - s)^{\alpha - 1} \, ds. \end{align}\] Similarly for \(\tilde{I}_2\) we get: \[\tilde{I}_2 = \frac{B}{\alpha - 1} \int_a^x h'(\tau) \left[ \frac{b - \tau}{\tau - a} \right]^{\frac{2 - \alpha}{2}} (x - \tau)^{\alpha - 1} \, d\tau.\] Finally we get the linear term by integrating the constant \(C_1\): \(\int_a^x C_1 dt=C_1(x-a)\).↩︎