An efficient numerical method for a Riemann–Liouville two-point boundary value problem

Abstract

In this paper a numerical method is considered for a two-point boundary value problem with a Riemann–Liouville fractional derivative, where the exact solution may have weak singularity. The linear interpolation is used to approximate the functions in the fractional integral transformed from the Riemann–Liouville boundary value problem. In order to capture the singular phenomena of the exact solution, an adaptive mesh is developed by equidistributing a monitor function. The stability is derived by a modified Grönwall inequality. It is shown that the scheme is second-order convergent. Numerical experiments are provided to demonstrate the theoretical results.

Introduction

This article is prompted by recent publications [1], [2] where the authors consider the following Riemann–Liouville two-point boundary value problem

$$-D^{2-\delta }u\left(x\right)+{\left(bu\right)}^{\prime }\left(x\right)+c\left(x\right)u\left(x\right)=f\left(x\right),x\in (0,1),$$$$u\left(0\right)=0,\alpha u\left(1\right)+\beta u^{\prime }\left(1\right)=\gamma ,$$

where $0<\delta <1$, $b^{\prime },c,f\in C^{q,\delta }(0,1]$ for some $q\in \mathbb{N}$, $c\ge 0$, $b^{\prime }+c\ge 0$, $\alpha ,\beta ,\gamma$ are given constants, and $D^{2-\delta }$ denotes a Riemann–Liouville fractional derivative defined by

$$D^{2-\delta }u\left(x\right)=\frac{1}{\Gamma \left(\delta \right)}\frac{d^2}{d{x}^2}{\int }_0^x{\left(x-t\right)}^{\delta -1}u\left(t\right)\mathrm{d}t.$$

Here $C^{q,\delta }(0,1]$ denotes the space of functions $v\in C[0,1]\cap C^q(0,1]$ such that

$$\left|v^{\left(k\right)}\left(x\right)\right|\le C\left(1+x^{1-\delta -k}\right)\mathrm{for}k=0,1,\dots ,q\mathrm{and}x\in (0,1]$$

with a positive constant $C$. It is proved in [1, Theorem 3] that the problem (1.1)–(1.2) has a unique solution $u\in C^{q,\delta }(0,1]$, which implies

$$\left|u^{\left(k\right)}\left(x\right)\right|\le C\left(1+x^{1-k-\delta }\right)\mathrm{for}k=0,1,\dots ,q\mathrm{and}x\in (0,1].$$

The bounds of the derivatives exhibit a typical singularity at $x=0$. Such a problem arises in modeling anomalous diffusion processes [3], [4].

In order to obtain a reliable numerical solution for problem (1.1)–(1.2), it is advantageous to use a mesh that concentrates nodes near the singularity. One approach is the use of graded meshes [1], [5]. Another approach is the use of adaptive meshes generated by equidistributing monitor functions over the domain of the problem [2]. Nowadays, the adaptive mesh approaches based on the equidistribution principle have attracted considerable attention from researchers since they are known to produce optimal rates of convergence [6], [7], [8], [9], [10]. In [2] an equidistributing mesh based on a posteriori error analysis has been used to solve problem (1.1)–(1.2) and it is shown that the scheme is first-order convergent. In this article, an integral discrete scheme based on the linear interpolation approximation and an adaptive mesh is used to solve problem (1.1)–(1.2). The adaptive mesh follows from equidistribution principle, which is based on a priori error estimate for capturing the singular behavior of the exact solution. It is shown that the scheme is second-order convergent, which improves the numerical results given in [2].

Notation. Throughout the paper, $C$ will denote a generic positive constant that is independent of the discretization parameter $N$. Note that $C$ is not necessarily the same at each occurrence. To simplify the notation we set $g_i=g\left(x_i\right)$ for any function $g$ defined on the interval $[0,1]$. We use the (pointwise) maximum norm on the interval $[0,1]$ by $\Vert \cdot \Vert$.