Numerical approximations for the variable coefficient fractional diffusion equations with non-smooth data

TL;DR

This paper develops a spectral approximation scheme for variable coefficient fractional diffusion equations, transforming them into constant coefficient forms to analyze stability and error, with validation through numerical experiments.

Contribution

It introduces a novel spectral method for variable coefficient fractional diffusion equations by transforming them into constant coefficient problems, providing error estimates for smooth and non-smooth data.

Findings

Spectral scheme achieves optimal convergence rates.

Error estimates are valid for non-smooth coefficients.

Numerical experiments confirm theoretical predictions.

Abstract

In this article we study the numerical approximation of a variable coefficient fractional diffusion equation. Using a change of variable, the variable coefficient fractional diffusion equation is transformed into a constant coefficient fractional diffusion equation of the same order. The transformed equation retains the desirable stability property of being an elliptic equation. A spectral approximation scheme is proposed and analyzed for the transformed equation, with error estimates for the approximated solution derived. An approximation to the unknown of the variable coefficient fractional diffusion equation is then obtained by post processing the computed approximation to the transformed equation. Error estimates are also presented for the approximation to the unknown of the variable coefficient equation with both smooth and non-smooth diffusivity coefficient and right-hand side. Three numerical experiments are given whose convergence results are in strong agreement with the theoretically derived estimates.

Tables

Table 1. Table 4.1: Convergence properties of Experiment 1.

$N$	${\Vert u-u_N\Vert }_{{\omega }^{(-\left(\alpha -\beta \right),-\beta )}}$	$\kappa$	${\Vert D\left(u-u_N\right)\Vert }_{{\omega }^{(-\left(\alpha -\beta \right)+1,-\beta +1)}}$	$\kappa$	${\Vert u-u_N\Vert }_{L^{\infty }\left(I\right)}$	$\kappa$
16	4.87E-04		1.15E-02		4.41E-04
20	3.09E-04	2.15	8.93E-03	1.18	2.95E-04	1.89
24	2.12E-04	2.16	7.26E-03	1.18	2.00E-04	2.23
28	1.54E-04	2.16	6.09E-03	1.19	1.54E-04	1.78
32	1.16E-04	2.16	5.22E-03	1.19	1.21E-04	1.86
36	9.08E-05	2.16	4.56E-03	1.19	9.46E-05	2.14
Pred.		2.20		1.20		1.20

Table 2. Table 4.2: Convergence properties of Experiment 2.

$N$	${\Vert u-u_N\Vert }_{{\omega }^{(-\left(\alpha -\beta \right),-\beta )}}$	$\kappa$	${\Vert D\left(u-u_N\right)\Vert }_{{\omega }^{(-\left(\alpha -\beta \right)+1,-\beta +1)}}$	$\kappa$	${\Vert u-u_N\Vert }_{L^{\infty }\left(I\right)}$	$\kappa$
16	4.57E-04		1.07E-02		4.41E-04
20	2.88E-04	2.20	8.29E-03	1.21	2.95E-04	1.90
24	1.96E-04	2.19	6.71E-03	1.21	2.00E-04	2.24
28	1.42E-04	2.19	5.61E-03	1.21	1.53E-04	1.78
32	1.07E-04	2.19	4.80E-03	1.21	1.20E-04	1.87
36	8.33E-05	2.19	4.18E-03	1.21	9.42E-05	2.15
Pred.		2.20		1.20		1.20

Table 3. Table 4.3: Convergence properties of Experiment 3.

$N$	${\Vert u-u_N\Vert }_{{\omega }^{(-\left(\alpha -\beta \right),-\beta )}}$	$\kappa$	${\Vert D\left(u-u_N\right)\Vert }_{{\omega }^{(-\left(\alpha -\beta \right)+1,-\beta +1)}}$	$\kappa$	${\Vert u-u_N\Vert }_{L^{\infty }\left(I\right)}$	$\kappa$
16	4.99E-04		1.14E-02		5.73E-04
20	3.11E-04	2.24	8.74E-03	1.25	3.71E-04	2.05
24	2.10E-04	2.24	7.03E-03	1.24	2.63E-04	1.98
28	1.51E-04	2.23	5.85E-03	1.24	1.88E-04	2.25
32	1.13E-04	2.22	4.99E-03	1.24	1.35E-04	2.59
36	8.80E-05	2.22	4.33E-03	1.23	1.07E-04	2.01
Pred.		1.20		1.20		1.20

Equations

$$\displaystyle-D\big{(}\big{(}r\,{}_{0}I_{x}^{2-\alpha}+(1-r)\,{}_{x}I_{1}^{2-\alpha}\big{)}\,K(x)\,Du(x)\big{)}$$

$$u(0)=u(1)=0.$$

$${}_0{I}_x^{\sigma }w(x):=\frac{1}{\Gamma (\sigma )}{\int }_0^{x}w(s)(x-s{)}^{\sigma -1}ds,{}_x{I}_1^{\sigma }w(x):=\frac{1}{\Gamma (\sigma )}{\int }_x^{1}w(s)(s-x{)}^{\sigma -1}ds,$$

$$L_{\omega}^{2}(\Omega)\,:=\,\bigg{\{}f(x)\,:\,\|f\|_{\omega}^{2}:=\int_{\Omega}\omega(x)\,f(x)^{2}\,dx\ <\ \infty\bigg{\}},\quad(f\,,\,g)_{\omega}\,:=\,\int_{\Omega}\omega(x)\,f(x)\,g(x)\,dx\,.$$

$$\displaystyle H^{m}_{\omega^{(\alpha,\beta)}}(\Omega):=\bigg{\{}v\,:\|v\|_{m,\omega^{(\alpha,\beta)}}^{2}:=\sum_{j=0}^{m}\big{|}v\big{|}_{j,\omega^{(\alpha,\beta)}}^{2}=\sum_{j=0}^{m}\big{\|}D^{j}v\big{\|}_{\omega^{(\alpha+j,\beta+j)}}^{2}<\infty\bigg{\}}.$$

$$P_n^{(\alpha ,\beta )}(x):=\sum _{k=0}^n{p}_{n,k}(x-1{)}^{n-k}(x+1{)}^k,x\in (-1,1),p_{n,k}:=\frac{1}{2^n}\left(\genfrac{}{}{0ex}{}{n+\alpha }{k}\right)\left(\genfrac{}{}{0ex}{}{n+\beta }{n-k}\right).$$

$$G_n^{(\alpha ,\beta )}(x):=P_n^{(\alpha ,\beta )}(2x-1),x\in [0,1],$$

$${\int }_0^1{\omega }^{(\alpha ,\beta )}(x)G_j^{(\alpha ,\beta )}(x)G_k^{(\alpha ,\beta )}(x)dx={\delta }_{j,k}|\Vert G_j^{(\alpha ,\beta )}|{\Vert }^2,{\omega }^{(\alpha ,\beta )}(x):=(1-x{)}^{\alpha }{x}^{\beta },$$

$$|\|G_{j}^{(\alpha,\beta)}|\|:=\Big{(}\frac{1}{(2j\,+\,\alpha\,+\,\beta\,+1)}\frac{\Gamma(j+\alpha+1)\,\Gamma(j+\beta+1)}{\Gamma(j+1)\,\Gamma(j+\alpha+\beta+1)}\Bigr{)}^{1/2}=|\|G_{j}^{(\beta,\alpha)}|\|.$$

$$\begin{array}[]{l}\displaystyle D^{k}G_{n}^{(\alpha,\beta)}(x)=\frac{\Gamma(n+k+\alpha+\beta+1)}{\Gamma(n+\alpha+\beta+1)}G_{n-k}^{(\alpha+k,\beta+k)}(x),\quad 0\leq k\leq n;\\[10.84006pt] \displaystyle D^{k}\Big{(}\omega^{(\alpha+k,\beta+k)}(x)G_{n-k}^{(\alpha+k,\beta+k)}(x)\Big{)}=\frac{(-1)^{k}\,n!}{(n-k)!}\,\omega^{(\alpha,\beta)}(x)\,G_{n}^{(\alpha,\beta)}(x),\quad 0\leq k\leq n.\end{array}$$

$$\frac{1}{2}\le \frac{|\Vert G_j^{(\alpha -\beta ,\beta )}|{\Vert }^2}{|\Vert G_{j+1}^{(\beta -1,\alpha -\beta -1)}|{\Vert }^2}=\frac{j+1}{j+\alpha }\le 1,j\ge 0.$$

$$\begin{array}[]{c}\displaystyle\int_{-1}^{1}\tilde{\omega}^{(\alpha,\beta)}(x)P_{j}^{(\alpha,\beta)}(x)P_{k}^{(\alpha,\beta)}(x)dx=\delta_{j,k}\big{\|}P_{j}^{(\alpha,\beta)}\big{\|}^{2}_{\tilde{\omega}^{(\alpha,\beta)}},\ \ \mbox{ where }\ \ \tilde{\omega}^{(\alpha,\beta)}(x):=(1-x)^{\alpha}(1+x)^{\beta},\\[7.22743pt] \mbox{ Using }\ \ \big{\|}P_{j}^{(\alpha,\beta)}\big{\|}_{\tilde{\omega}^{(\alpha,\beta)}}=\Big{(}\frac{2^{(\alpha+\beta+1)}}{(2j\,+\,\alpha\,+\,\beta\,+1)}\frac{\Gamma(j+\alpha+1)\,\Gamma(j+\beta+1)}{\Gamma(j+1)\,\Gamma(j+\alpha+\beta+1)}\Big{)}^{1/2}=\big{\|}P_{j}^{(\beta,\alpha)}\big{\|}_{\tilde{\omega}^{(\beta,\alpha)}}\end{array}$$

$${\int }_{-1}^1{\tilde{\omega }}^{(\alpha ,\beta )}(x)P_j^{(\alpha ,\beta )}(x)P_k^{(\alpha ,\beta )}(x)dx=2^{\alpha +\beta +1}{\int }_0^1{\omega }^{(\alpha ,\beta )}(x)G_j^{(\alpha ,\beta )}(x)G_k^{(\alpha ,\beta )}(x)dx.$$

$$\begin{array}[]{l}\displaystyle D^{k}P_{n}^{(\alpha,\beta)}(x)=\frac{\Gamma(n+k+\alpha+\beta+1)}{2^{k}\,\Gamma(n+\alpha+\beta+1)}P_{n-k}^{(\alpha+k,\beta+k)}(x),\\[10.84006pt] \displaystyle D^{k}\Big{(}\tilde{\omega}^{\alpha+k,\beta+k}(x)P_{n-k}^{(\alpha+k,\beta+k)}(x)\Big{)}=\frac{(-1)^{k}\,2^{k}\,n!}{(n-k)!}\,\tilde{\omega}^{\alpha,\beta}(x)P_{n}^{(\alpha,\beta)}(x).\end{array}$$

$$\big{(}v\,-\,P_{N,a,b}v\ ,\ \phi_{N}\big{)}_{\omega^{(a,b)}}\ =\ 0\,,\ \ \forall\phi_{N}\in\mathcal{S}_{N}.$$

$$\big{\|}v-P_{N,a,b}v\|_{\mu,\omega^{(a,b)}}\ \leq\ C\,\big{(}N\,(N+a+b)\big{)}^{\frac{\mu-t}{2}}\,|v|_{t,\omega^{(a,b)}}.$$

$$-D\big{(}\big{(}r\,{}_{0}I_{x}^{2-\alpha}+(1-r)\,{}_{x}I_{1}^{2-\alpha}\big{)}\,D\underbrace{D^{-1}\,K(x)\,Du(x)}_{:=\,\widetilde{w}(x)}\big{)}\ =\ f(x)\,.$$

$$D\tilde{w}(x)$$

$$⟹u(x)$$

$$\text{mbox}Now,u(1)=0⟹{\int }_0^1\frac{D\tilde{w}(s)}{K(s)}ds$$

$$\displaystyle\mathcal{L}_{r}^{\alpha}\widetilde{w}(x)\ :=\ -D\big{(}\big{(}r\,{}_{0}I_{x}^{2-\alpha}+(1-r)\,{}_{x}I_{1}^{2-\alpha}\big{)}\,D\widetilde{w}(x)\big{)}$$

$$\text{mbox}subjectto\tilde{w}(0)=0\text{mbox}and{\int }_0^x\frac{D\tilde{w}(s)}{K(s)}ds$$

$${\mathcal{L}}_r^{\alpha }w(x)$$

$$\text{mbox}subjecttow(0)$$

$$\Vert w{\Vert }_{{\omega }^{(-(\alpha -\beta ),-\beta )}}+\Vert Dw{\Vert }_{{\omega }^{(-(\alpha -\beta )+1,-\beta +1)}}\le C\Vert f{\Vert }_{{\omega }^{(\beta ,\alpha -\beta )}}$$

$$\text{mbox}Let\tilde{w}(x)=C_1{k}_1(x)+w(x).$$

$$\displaystyle\int_{0}^{1}\frac{1}{K(s)}\big{(}C_{1}\,(1-s)^{\alpha-\beta-1}\,s^{\beta-1}$$

$$⟹C_1$$

$$D(\frac{1}{K})\in L_{{\omega }^{(\alpha -\beta ,\beta )}}^2$$

$$=\frac{{\int }_0^{1}w(s)\frac{K^{\prime }(s)}{K(s{)}^2}ds}{{\int }_0^1\frac{(1-s{)}^{\alpha -\beta -1}{s}^{\beta -1}}{K(s)}ds}.$$

$$den:={\int }_0^1\frac{(1-s{)}^{\alpha -\beta -1}{s}^{\beta -1}}{K(s)}ds,c_1:={\int }_0^1\frac{-Dw(s)}{K(s)}ds.$$

$$|C_1|$$

$$\displaystyle=\frac{1}{den}\ \int_{0}^{1}\omega^{((-(\alpha-\beta)+1)/2\,,\,(-\beta+1)/2)}(s)\,\big{|}Dw(s)\big{|}\,\omega^{((\alpha-\beta-1)/2\,,\,(\beta-1)/2)}(s)\,\frac{1}{K(s)}\,ds$$

$$\displaystyle\leq\frac{1}{den}\ \|Dw\|_{\omega^{(-(\alpha-\beta)+1\,,\,-\beta+1)}}\,\bigg{\|}\frac{1}{K(s)}\bigg{\|}_{\omega^{((\alpha-\beta-1)\,,\,(\beta-1))}}$$

$$\le C\Vert f{\Vert }_{{\omega }^{(\beta ,\alpha -\beta )}}\text{mbox}(using\text{mbox}Theorem\text{ref}thmWellposedw).$$

$$u(x)=C_1{\int }_0^x\frac{(1-s{)}^{\alpha -\beta -1}{s}^{\beta -1}}{K(s)}ds+{\int }_0^x\frac{Dw(s)}{K(s)}ds,$$

Full text

Key words. Fractional diffusion equation, Jacobi polynomials, spectral method

AMS Mathematics subject classifications. 65N30, 35B65, 41A10, 33C45

1 Introduction

It has been shown that fractional partial differential equations (PDEs) can accurately model challenging phenomena including anomalous transport, long-range time memory and spatial interactions [1, 14]. Extensive research has been conducted on fractional PDEs in terms of their modeling, analysis, numerical approximations and applications. In a representative piece of work Ervin and Roop [3] proved the wellposedness of the Galerkin weak formulation of linear elliptic space fractional diffusion equations (FDEs) of order $1<\alpha<2$ on the Sobolev space $H^{\alpha/2}_{0}$. They also proved optimal-order error estimates of its finite element approximations in the energy and $L^{2}$ norms, assuming that the true solution and the solution to the dual problem for an $L^{2}$ right-hand side have full regularity. However, it was later realized that the smoothness of the coefficients and source term for these space fractional differential equations cannot ensure the smoothness of their solution [4, 9, 19, 20]. This is in sharp contrast to integer-order linear elliptic PDEs [5, 6]. For this reason, the usual smoothness assumptions on the true solutions to fractional PDEs in the analysis of the numerical approximations are inappropriate.

It turns out that the spectral methods are particularly well suited for the accurate approximation of FDEs, as they provide a clean expression of the true solution to FDEs for the convenience of analysis [4, 12], and ordinarily lead to a diagonal stiffness matrices (at least for constant-coefficient FDEs). This is in contrast to the dense stiffness matrices generated from the finite element, finite difference, or finite volume approximations. Mao et al. [11] analyzed the regularity of the solution to a symmetric case of the FDE and developed corresponding spectral methods. The solution structure to the general case was resolved completely in [4], in which a spectral method utilizing the weighted Jacobi polynomial was studied and a priori error estimates derived. The two-sided FDE with constant coefficient and Riemann-Liouville fractional derivative was investigated in [12], by employing a Petrov-Galerkin projection in a properly weighted Sobolev space using two-sided Jacobi polyfractonomials as test and trial functions. In [23] and [24], the regularity of the two-sided fractional reaction-diffusion and advection-diffusion-reaction equations are analyzed in the weighted Sobolev spaces, based on which the optimal (or sub-optimal) convergence rates of the spectral Galerkin or Petrov-Galerkin method are proved.

The variable diffusivity $K$ presents another bottleneck of FDEs. It was shown in [18] that the Galerkin weak formulation may lose its coercivity for a smooth $K(x)$ with positive lower and upper bounds, which increases the difficulties for the stability and convergence analysis and accurate simulations. To circumvent this issue, an indirect Legendre spectral Galerkin method was developed for the FDE in [20], in which the high-order convergence rates of numerical approximations were proved only under the regularity assumptions of coefficients and right-hand side term. In [10], with the introduction of an auxiliary variable, a mixed approximation was developed for an FDE and the corresponding error estimates were proved. In [13], a spectral Galerkin method for a different variable coefficient FDE was analyzed, in which the outside and inside fractional derivatives are chosen carefully so that the corresponding Galerkin weak formulation are self-adjoint and coercive.

Recently the wellposedness of the variable coefficient FDE (2.1) was investigated in [22], in which the existence and uniqueness of the solution to the proposed model was proven for any $f\in L^{2}_{\omega^{(\beta-1,\alpha-\beta-1)}}$, with the space defined by (2.3). A spectral approximation method was then studied and several error estimates were derived based on the regularity of the right-hand side term. In this paper we continue to investigate model (2.1) using a different approach than that used in [22]. We prove in this paper that the model is wellposed for $f$ belonging to a larger space $L^{2}_{\omega^{(\beta,\alpha-\beta)}}$, which extends the wellposedness results in [22]. A spectral approximation scheme is proposed and the error estimates are proved to be dependent on the weaker norms of $f$ without loss of accuracy. In addition, we follow the idea of the K-method of interpolation to determine the range of the index of the weighted Sobolev space that the power function belongs to, which provides the theoretical support for estimating the convergence rate of the proposed method.

This paper is organized as follows. In Section 2 we present the formulation of the model and introduce notation and key lemmas used in the analysis. The wellposedness of the model and the regularity of its solution are studied in Section 3, based on which the spectral approximation method is formulated and a detailed analysis of its convergence is proved. Three numerical experiments are presented in Section 4 whose results demonstrate the sharpness of the derived error estimates.

2 Model problem and preliminaries

In this paper we consider the following homogeneous Dirichlet boundary-value problem of a two-sided Caputo flux FDE, which is obtained by incorporating a fractional Fick’s law into a conventional local mass balance law with a variable fractional diffusivity [2, 21]:

[TABLE]

Here $1<\alpha<2$, $D$ is the first-order differential operator, $K(x)$ is the fractional diffusivity with $0<K_{m}\leq K(x)\leq K_{M}<\infty$, $0\leq r\leq 1$ indicates the relative weight of forward versus backward transition probability and $f(x)$ the source or sink term. The left and right fractional integrals of order $0<\sigma<1$ are defined as [15]

[TABLE]

where $\Gamma(\cdot)$ is the Gamma function.

We introduce notation and properties used subsequently in our discussion of the approximation scheme and in its error analysis.

Let $\Omega$ be a bounded open interval and $\omega(x)>0,~{}x\in\Omega$ be a smooth function. We define the weighted $L^{2}$ space, $L_{\omega}^{2}(\Omega)$, and $L^{2}$ weighted inner product as

[TABLE]

In addition, let $\mathbb{N}_{0}:=\mathbb{N}\cup{0}$, and $\omega^{(\alpha,\beta)}$ be a weighting function defined on $\Omega$ and indexed by $\alpha$ and $\beta$. For any $m\in\mathbb{N}$, we introduce the following weighted Sobolev spaces [7, 16]

[TABLE]

For $t\in\mathbb{R}^{+}\backslash\mathbb{N}_{0}$, $H^{t}_{\omega^{(a\,,\,b)}}(0,1)$ is defined by the $K$-method of interpolation, and for $t\in\mathbb{R}^{-}$, $H^{t}_{\omega^{(a\,,\,b)}}(0,1)$ is defined by duality.

The Jacobi polynomials $P_{n}^{(\alpha,\beta)}(x)$ are defined by [16, 17]

[TABLE]

Let $G_{n}^{(\alpha,\beta)}(x)$ denote the translated and dilated Jacobi polynomials to the interval $[0,1]$:

[TABLE]

and summarize the properties of $G_{n}^{\alpha,\beta}$ in the following lemma

Lemma 2.1

For $\omega^{(\alpha,\beta)}(x)=(1-x)^{\alpha}x^{\beta}$, the polynomials $G_{n}^{(\alpha,\beta)}$ have the following orthogonality and norm properties

[TABLE]

where $\delta_{j,k}=1$ if $j\neq k$ and 0 otherwise, and

[TABLE]

In addition, $G_{n}^{(\alpha,\beta)}$ satisfies

[TABLE]

Finally, $G_{n}^{(\alpha,\beta)}$ have the following norm relation

[TABLE]

Proof. The orthogonality property (2.6) of $G_{n}^{(\alpha,\beta)}$ is a direct consequence of the orthogonality relation of $P_{n}^{(\alpha,\beta)}$ [16, 17]

[TABLE]

and the following relation between $G_{n}^{(\alpha,\beta)}$ and $P_{n}^{(\alpha,\beta)}$ leads to (2.7).

[TABLE]

The two equations in (2.8) are direct consequences of the following equations for $P_{n}^{(\alpha,\beta)}(x)$ [11, equations (2.15) and (2.19)]

[TABLE]

The norm relation (2.9) is derived in [22].

Let $\mathcal{S}_{N}$ denote the space of polynomials of degree $\leq N$. We define the weighted $L^{2}$ orthogonal projection $P_{N,a,b}:\,L^{2}_{\omega^{(a,b)}}(0,1)\rightarrow\mathcal{S}_{N}$ by the condition

[TABLE]

Lemma 2.2

[7, Theorem 2.1]** For $\mu\in\mathbb{N}_{0}$ and $v\in H^{t}_{\omega^{(a\,,\,b)}}(0,1)$, with $0\leq\mu\leq t$, there exists a constant $C$, independent of $N,\,\alpha$ and $\beta$ such that

[TABLE]

Remark: In [7] (2.11) is stated for $t\in\mathbb{N}_{0}$. The result extends to $t\in\mathbb{R}^{+}$ using interpolation.

3 Approximation scheme

3.1 Motivation for the approximation scheme

Introduce $D^{-1}\,:\,L^{1}(\mathrm{I})\longrightarrow H^{1}(\mathrm{I})$, defined by $D^{-1}g(x)\,:=\,\int_{0}^{x}\,g(s)\,ds$.

Rewrite (2.1) as

[TABLE]

Consider $\widetilde{w}(x)\ =\ D^{-1}\,K(x)\,Du(x)\ =\ \int_{0}^{x}K(s)\,D\,u(s)\,ds$. Note that $\widetilde{w}(0)=0$. Then,

[TABLE]

Hence if we could determine $\widetilde{w}(x)$ such that

[TABLE]

then our solution to (2.1),(2.2) would be given by (3.1). With this in mind, consider the problem: Determine $w(x)$ satisfying

[TABLE]

Let $\beta\in[\alpha-1,1]$ be determined by $(1-r)\sin(\pi\beta)=r\sin(\pi(\alpha-\beta))$. The following theorem ensures the well-posedness of this problem.

Theorem 3.1** ([4, 8])**

Let $f(x)\in L^{2}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$. Then, there exists a unique solution $w(x)\in L^{2}_{\omega^{(-(\alpha-\beta),-\beta)}}(\mathrm{I})$ satisfying (3.3),(3.4). In addition, there exists $C>0$ such that

[TABLE]

From [4] we have that $k_{1}(x)\,:=\,\int_{0}^{x}(1-s)^{\alpha-\beta-1}\,s^{\beta-1}\,ds\,\in\,ker(\mathcal{L}_{r}^{\alpha}).$

[TABLE]

We have that $\mathcal{L}_{r}^{\alpha}\widetilde{w}(x)\ =\ f(x)$, and $\widetilde{w}(0)\,=\,0$.

The condition (3.2) combined with (3.6) implies

[TABLE]

Let

[TABLE]

Then $C_{1}=c_{1}/den$ can be bounded by

[TABLE]

Combining (3.1),(3.6), (3.7), (3.11) and Theorem 3.1 we have the following.

Theorem 3.2

For $f(x)\in L^{2}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$ and there exists a unique solution $u(x)\in L^{\infty}(\mathrm{I})$ to (2.1),(2.2), given by

[TABLE]

where $w(x)$ is determined by (3.3),(3.4) and $C_{1}$ by (3.7).

Additionally, for $\epsilon_{1},\epsilon_{2}>0$ there exists $C>0$ such that

[TABLE]

Proof:

It is straightforward to show that $u(x)$ given by (3.12) satisfies (2.1),(2.2). Next we show that there exists a unique solution to (2.1),(2.2).

Assume that $u_{1}(x)$ and $u_{2}(x)$ are solutions of (2.1),(2.2). Let $z_{1}(x)$ and $z_{2}(x)$ be defined by:

[TABLE]

Note that $\mathcal{L}_{r}^{\alpha}(z_{1}\,-\,z_{2})\ =\ 0$. Hence, $(z_{1}\,-\,z_{2})\in ker(\mathcal{L}_{r}^{\alpha})$. Thus from [4], for constants $A$ and $B$,

[TABLE]

Using (3.10) and (3.11)

[TABLE]

Consequently, we obtain

[TABLE]

Estimate (3.13) then follows from (3.15) and (3.16).

3.2 Approximation scheme

To compute an approximation to $u(x)$, $u_{N}(x)$, we firstly compute an approximation to $w(x)$, $w_{N}(x)$, satisfying (3.3),(3.4), and then use $w_{N}(x)$ in place of $w(x)$ in (3.7) and (3.12) to obtain $u_{N}(x)$.

3.2.1 Approximation of $w(x)$ satisfying (3.3),(3.4)

Proceeding as in [4], $f(x)\in L^{2}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$ may be expressed as $f(x)\ =\ \sum_{i=0}^{\infty}\frac{f_{i}}{|\|G_{i}^{(\beta,\alpha-\beta)}|\|^{2}}\,G_{i}^{(\beta,\alpha-\beta)}(x)$, where $f_{i}$ is given by

[TABLE]

With $f_{i}$ defined in (3.17), let

[TABLE]

[TABLE]

Using Stirling’s formula we have that

[TABLE]

Theorem 3.3** ([4, 8])**

Let $f(x)\in L^{2}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$ and $w_{N}(x)$ be as defined in (3.18). Then, $w(x)\ :=\ \lim_{N\rightarrow\infty}w_{N}(x)\ =\ \omega^{(\alpha-\beta,\beta)}(x)\,\sum_{j=0}^{\infty}c_{j}\,G_{j}^{(\alpha-\beta,\beta)}(x)\,\in L^{2}_{\omega^{(-(\alpha-\beta),-\beta)}}(\mathrm{I})$ and satisfies (3.3),(3.4).

Theorem 3.4

For $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}}(\mathrm{I})$, $t\geq 0$, and $w_{N}(x)$ given by (3.18), there exists $C>0$ such that

[TABLE]

Proof: Using the definition of the $\|\cdot\|_{\omega^{(-(\alpha-\beta)\,,-\beta)}}$ norm,

[TABLE]

Similarly, using (2.8)

[TABLE]

Using (3.20),

[TABLE]

Combining (3.23) and (3.24),

[TABLE]

3.2.2 Approximation of $u(x)$ satisfying (3.12)

The approximation $u_{N}(x)$ of $u(x)$ is obtained by substituting $w_{N}(x)$ in place of $w(x)$ in (3.7) and (3.12). With $den$ defined in (3.9), let

[TABLE]

Note that $|C_{1}\,-\,C_{1,N}|\ =\ |c_{1}\,-\,c_{1,N}|/den$. Hence the rate of convergence as $N\rightarrow\infty$ of $|C_{1}\,-\,C_{1,N}|$ is equal to the rate of convergence of $|c_{1}\,-\,c_{1,N}|$. Now,

[TABLE]

In case $D(\frac{1}{K})\in L^{2}_{\omega^{(\alpha-\beta\,,\,\beta)}}(\mathrm{I})$ ,

[TABLE]

We have the following error estimates for $u\,-\,u_{N}$.

Theorem 3.5

For $f\in H^{t}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$, $t\geq 0$, then for $\epsilon_{1},\epsilon_{2}>0$ there exists $C>0$ (independent of $N$ and $\alpha$) such that

[TABLE]

Proof: From (3.12), we have using (3.25) and (3.26)

[TABLE]

Then, from (3.32) and

[TABLE]

we obtain (3.29).

For $D(\frac{1}{K})\in L^{2}_{\omega^{(\alpha-\beta\,,\,\beta)}}(\mathrm{I})$, we apply integration by parts to (3.31) to obtain

[TABLE]

Therefore,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Combining (3.33)-(3.36) with (3.28) and (3.21) we obtain (3.30).

We conclude this section with an error bound for $D(u-u_{N})$.

Lemma 3.1

For $f\in H^{t}_{\omega^{(\beta,\alpha-\beta)}}(\mathrm{I})$, $t\geq 0$, then there exists $C>0$ (independent of $N$ and $\alpha$) such that

[TABLE]

Proof: From (3.31) it follows that

[TABLE]

Thus,

[TABLE]

where in the last step we have used (3.26) and (3.22)), and the fact that $(\alpha-\beta-1)$ and $(\beta-1)\,>\,-1$.

4 Numerical experiments

In this section we present three numerical experiments to demonstrate our approximation scheme, and to compare the experimental rate of convergence of the approximation with the theoretically predicated rate.

Numerical example. Let $K(x)=1/(1+x^{\gamma})$ and

[TABLE]

Then the solution $u(x)$ is given by (3.12) where

[TABLE]

and ${}_{2}F_{1}(a,b;c,x)$ denotes the Gaussian three parameter hypergeometric function.

In order to determine the theoretical rate of convergence for $\|u-u_{N}\|_{\omega^{(-(\alpha-\beta),-\beta)}}$, $\|u-u_{N}\|_{L^{\infty}}$ and $\|D(u-u_{N})\|_{\omega^{(-(\alpha-\beta)+1\,,\,-\beta+1)}}$ from (3.29), (3.30), and (3.37), respectively, we need to determine the largest value for $t$ such that $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}}(\mathrm{I})$. The most singular terms for $f(x)$ in (4.1) are $x^{1-\alpha}$ and $(1-x)^{1-\alpha}$. Using Lemma A.1 (in the Appendix) we have that $x^{1-\alpha}\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}}(\mathrm{I})$, for $t<\ 3-\alpha-\beta$, and $(1-x)^{1-\alpha}\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}}(\mathrm{I})$, for $t<\ 3-\alpha-(\alpha-\beta)$.

Then, for Experiment 1 ($\alpha=1.60$, $\beta=0.80$) $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}\,}(\mathrm{I})$ for $t\,<\,3-\alpha-\max\{\alpha-\beta\,,\beta\}=0.90$, which leads to theoretical asymptotic convergence rates of $\|u-u_{N}\|_{\omega^{(-(\alpha-\beta),-\beta)}}\sim N^{-2.20}$ (using (3.30)), $\|u-u_{N}\|_{L^{\infty}}\sim N^{-1.20}$ (using (3.29)) and $\|D(u-u_{N})\|_{\omega^{(-(\alpha-\beta)+1\,,\,-\beta+1)}}\sim N^{-1.20}$ (using (3.37)).

Assuming that $\|\xi-\xi_{N}\|_{L_{\rho}}\sim\,N^{-\kappa}$, the experimental convergence rate is calculated using

[TABLE]

Experiment 1. In this experiment we select $\alpha=1.60$, $r=0.50$, $\beta=0.80$ and $\gamma=0.80$, which leads to $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}\,}(\mathrm{I})$ for $t\,<\,0.60$.

Experiment 2. In this experiment, we take $\alpha=1.30$, $r=0.63$, $\beta=0.50$ and $\gamma=0.80$. The above analysis gives that $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}\,}(\mathrm{I})$ for $t\,<\,0.90$. The corresponding theoretical asymptotic convergence rates are $\|u-u_{N}\|_{\omega^{(-(\alpha-\beta),-\beta)}}\sim N^{-2.20}$, $\|u-u_{N}\|_{L^{\infty}}\sim N^{-1.20}$ and $\|D(u-u_{N})\|_{\omega^{(-(\alpha-\beta)+1\,,\,-\beta+1)}}\sim N^{-1.20}$.

Experiment 3. In this experiment we select $\alpha=1.30$, $r=0.63$, $\beta=0.50$ and $\gamma=0.10$, which leads to $f(x)\in H^{t}_{\omega^{(\beta\,,\,\alpha-\beta)}\,}(\mathrm{I})$ for $t\,<\,0.90$. However, in this case, $D(\frac{1}{K})\not\in L^{2}_{\omega^{(\alpha-\beta\,,\,\beta)}}(\mathrm{I})$ due to the relatively strong singularity of $K(x)$ at $x=0$, which means that (3.28) is not applicable. Hence we can only apply the bound (3.26) of $c_{1}-c_{1,N}$, which leads to the estimate (3.29) of $u-u_{N}$ instead of (3.30), and consequently an estimate for the convergence rate of $\|u-u_{N}\|_{\omega^{(-(\alpha-\beta),-\beta)}}$ of $1.20$ by (3.29), instead of $2.20$ if using (3.30).

The experimental convergence rates for $\|u-u_{N}\|_{L^{2}_{\omega^{(-(\alpha-\beta),-\beta)}}}$ and $\|D(u-u_{N})\|_{L^{2}_{\omega^{(-(\alpha-\beta)+1,-\beta+1)}}}$ are in strong agreement with the theoretically predicted rates for the first two experiments. For Experiment 3, we note that the numerical convergence rate of $\|u-u_{N}\|_{L^{2}_{\omega^{(-(\alpha-\beta),-\beta)}}}$ is 2.20, which corresponds to the case for $D(\frac{1}{K})\in L^{2}_{\omega^{(\alpha-\beta\,,\,\beta)}}(\mathrm{I})$ (even though this is not the case in this experiment). Additionally, we remark that the error in the $L^{\infty}$ norm is difficult to measure accurately due to the oscillatory nature of the error function caused by the polynomial approximation, as illustrated by Figure 1.

Appendix A Appendix

In this section we investigate which $H^{s}_{\omega^{(\alpha,\beta)}}(\mathrm{I})$ space $u(x)=x^{\mu}$ lies in. For brevity of notation, in this section we use $H^{s}_{(\alpha,\beta)}(\mathrm{I})\equiv H^{s}_{\omega^{(\alpha,\beta)}}(\mathrm{I})$.

Lemma A.1

*Let $u(x)\,=\,x^{\mu}$. Then, $u\in H^{s}_{(\alpha,\beta)}(\mathrm{I})$ for $s\,\geq\,0$ satisfying $s\,<\,2\mu+\beta+1$.

Proof: Let $\chi(x)\in C^{\infty}[0,\infty)$ denote the cutoff function satisfying

[TABLE]

and let $\chi_{\delta}(x)\,:=\,\chi(\frac{x}{\delta})$, for $\delta>0$. Note that

[TABLE]

For $\delta$ to be determined, let $u\,=\,v+w$ where $v(x)\ =\ \chi_{\delta}(x)\,u(x)$ and w(x)\ =\ \big{(}1-\chi_{\delta}(x)\big{)}\,u(x). We have that

[TABLE]

Thus,

[TABLE]

which implies that, for $m\,<\,2\mu\,+\,\beta+1$, $v\in H^{m}_{(\alpha,\beta)}(\mathrm{I})$ and

[TABLE]

Next, consider $w(x)$.

[TABLE]

The first term on the RHS of (A.2) vanishes for $x<\delta/4$, and the second term vanishes for $x\,<\,\delta/4$ and $x\,>\,3\delta/4$. Using this,

[TABLE]

Hence, for $n\,>\,2\mu\,+\,\beta\,+\,1$

[TABLE]

(Remark: For $n\,>\,2\mu\,+\,\beta\,+\,1$ the exponent of $\delta$ in (A.6) is negative, so the ${}^{\prime}1^{\prime}$ term in (A.5) is bounded by the $\delta$ term.)

For $0<t<1$ we have from (A.1) and (A.6)

[TABLE]

Setting $\delta\ =\ t^{2/(n-m)}$ leads to $K(t,u)\,\leq\,C\,t^{(2\mu\,+\,\beta+1-m)/(n-m)}$.

Recall that

[TABLE]

The larger the value of $\theta$ $(0<\theta<1)$ in (A.9) such that the integral is finite, the “nicer” (i.e., more regular) is the function $u$. Hence from (A.9), we are interested in the integrand about $t=0$. We have trivially that for $u_{1}=u$, $u_{2}=0$ in (A.7) that $K(t,u)\,\leq\,\|u\|_{H^{m}_{(\alpha,\beta)}}\,\leq\,C$. Hence it follows that

[TABLE]

Using (A.9),

[TABLE]

For $s\ =\ (1-\theta)m\ +\ \theta n\ =\ m\ +\ \theta(n-m)$, then $s\ <\ 2\mu\,+\,\beta+1$ using (A.10)   .

Hence we can conclude that $u(x)\ =\ x^{\mu}\,\in\,H^{s}_{(\alpha,\beta)}(\mathrm{I})$ for $s<2\mu\,+\,\beta+1$.