text:
\begin{document}
\bibliographystyle{plain}

\baselineskip 13pt

\title[Approximation by Legendre expansions] {Optimal  error estimates for
 Legendre  expansions  of singular functions with fractional derivatives of bounded variation}
\author[
	W. Liu, \,   L. Wang\,  $\&$\, B. Wu
	]{\;\; Wenjie Liu$^\dagger$,   \;\;  Li-Lian Wang$^\ddagger$ \;\; and\;\; Boying Wu$^\dagger$
		}

	\thanks{${}^\dagger$Department of Mathematics, Harbin Institute of Technology, 150001, China.  The research of the first author was supported by the China Postdoctoral Science Foundation Funded Project (No. 2017M620113), the National Natural Science Foundation of China (Nos. 11801120, 71773024 and 11271100), the Fundamental Research Funds for the Central Universities (Grant No.HIT.NSRIF.2020081) and the Natural Science Foundation of Heilongjiang Province of China (Nos. A2016003 and G2018006). Emails: liuwenjie@hit.edu.cn (Wenjie Liu) and mathwby@hit.edu.cn (Boying Wu).\\
		\indent ${}^\ddagger${\em Corresponding author}.  Division of Mathematical Sciences, School of Physical
		and Mathematical Sciences, Nanyang Technological University,
		637371, Singapore. The research of the second author is partially supported by Singapore MOE AcRF Tier 2 Grants: MOE2018-T2-1-059 and MOE2017-T2-2-144. Email: lilian@ntu.edu.sg.

	}

\begin{abstract} We present a new fractional Taylor formula for
singular functions whose  Caputo fractional derivatives are of bounded variation.  It bridges and ``interpolates" the usual Taylor formulas with two consecutive integer orders.  This enables us to obtain an analogous formula for the Legendre expansion coefficient of this type of singular functions, and further derive the optimal (weighted) $L^\infty$-estimates and $L^2$-estimates of the Legendre polynomial approximations.  This set of results can enrich the existing theory for $p$ and $hp$ methods for singular problems, and
 answer some open questions posed in some recent literature.

\end{abstract}
\keywords{Approximation by Legendre polynomials,  functions with interior and endpoint singularities, optimal estimates, fractional Taylor formula}
 \subjclass[2000]{41A10, 41A25,  41A50,   65N35, 65M60}
\maketitle

\section{Introduction}

\begin{comment}
Among the vast family of orthogonal polynomials  in the name of  German mathematician Carl Gustav Jacob Jacobi (1804-1851)-- Jacobi polynomials,  the Chebyshev and Legendre polynomials are perhaps the most  widely used  in numerical analysis (particularly, in spectral methods).  Historically,  the former has been advocated and promoted  to a prevailing ``Chebyshev faith",  which asserts  that  the Chebyshev polynomials are better than any other set of Jacobi polynomials
(cf. \cite{Rivlin1969JAT-An,Fox1972Book,Mason2003Book}).   After all,  (i) the Chebyshev approximation is optimum in the maximum norm; (ii) the Chebyshev polynomials:
$T_n(x)=\cos (n\, \arccos x))$ and related Gauss-quadrature nodes/weights are explicitly known;  (iii)  more importantly,  they enjoy the fast Fourier transform (FFT)  among others.  In fact,  Boyd and Pestschek \cite{Boyd2014JSC}
provided  a fair comparison with  evidences to show that Legendre and other Jacobi polynomials  can win in several exceptional situations.  For example,
the ``spectral elements" \cite{Patera1984JCP-Spectral}  were initially Chebyshev-based,  but shortly swapped  to the Legendre polynomials for at least two reasons:   the unpleasant weight function: $1/\sqrt{1-x^2} $ brings about troublesome extra terms in the use of integration by parts;  and  FFT is overkill for low polynomial degree. As such, spectral elements using Legendre polynomials are predominant.
On the other hand, it is  noteworthy that in a single domain spectral method,  the use of {\em compact combination of Legendre polynomials} as basis functions
\cite{Shen1994SISC}, developed into  the generalized Jacobi polynomials/functions \cite{Guo06JSC-Optimal,Guo2009ANM,Chen2016MC},  led to optimal spectral algorithms for boundary value problems.

\end{comment}

The study of  Legendre approximation to singular functions  is a subject of fundamental importance in the theory and applications of $hp$ finite element methods. We refer to   the seminal series  of papers by
Gui and Babu\v{s}ka \cite{Gui1986NM-Part1,Gui1986NM-Part2,Gui1986NM-Part3}  and many other developments in e.g.,  \cite{Schwab1998Book,Babuska2000NM,Babuska2001SINUM}.
In particular, the very recent work of  Babu\v{s}ka and Hakula \cite{Babuska19CMAME-Pointwise}  provided a review of known/unknown results and posed
 a few open questions  on the pointwise error estimates of Legendre expansion of a typical  singular function discussed in
 \cite{Gui1986NM-Part1}:
 \begin{equation}\label{uxy1}
 u(x)=(x-\theta)^\mu_+=\begin{cases}
 0,  & -1< x\le \theta,\\[2pt]
 (x-\theta)^\mu,    & \theta< x<1,
 \end{cases}
 \quad  |\theta|\le 1,\;\;\; \mu>-1.
 \end{equation}
One significant development along this line is the $hp$ approximation theory in the framework of  Jacobi-weighted Besov spaces \cite{Babuska2000NM,Babuska2001SINUM,Babuska2002MMMAS,Guo2013JAT-Direct}.
Such Besov spaces are defined through space interpolation of Jacobi-weighted Sobolev spaces with integer regularity indices using the $K$-method.  It is important to point out  that the non-uniformly Jacobi-weighted Sobolev spaces has been employed in spectral approximation theory    \cite{Funaro1992Book,Schwab1998Book,Guo2004JAT,Guo06JSC-Optimal,Shen2011Book}.

\subsection{Related works} Different from  the Sobolev-Besov framework,  Trefethen  \cite{Trefethen2008SIREV,Trefethen2013Book} characterised
the regularity of singular functions by using the space of absolute continuity  and bounded variation (AC-BV),
in the study of Chebyshev expansions of such functions.

One motivative example therein is $u(x)=|x|$ in $\Omega=(-1,1)$ which has the regularity:   $u,u'\in {\rm AC}(\bar \Omega)$ and
$u''\in {\rm BV}(\bar \Omega)$ (where the integration of the norm  is in the Riemann-Stieltjes (RS) sense).
As a result,  the maximum error of its Chebyshev expansion can attain  optimal order  (but can only be suboptimal in a usual Sobolev framework). There have been many follow-up  works on the improved error estimates of  Chebyshev approximation or more general Jacobi polynomial approximation  under this AC-BV framework (see, e.g., \cite{Majidian2017ANM,Wang18JCAM-convergence,Wang18AML-new,Xiang20NM-Optimal}).  However, the regularity index and the involved derivatives are of integer order, so it is not suitable to best  characterise the regularity of many singular functions, e.g., \eqref{uxy1} and $u(x)=|x|^\mu$ with non-integer  $\mu.$ In other words,  if one naively applies the estimates,  then the loss of order might occur.
Nevertheless,  the solutions of many singular problems (in irregular domains or with singular coefficients/operators among others)
typically exhibit this kind of singularities.

 To fill this gap, we introduced  for the first time in \cite{Liu19MC-Optimal} certain fractional Sobolev-type spaces
and derived optimal Chebyshev polynomial approximation to  functions with interior and endpoint singularities within this new framework.  This study also inspired the discovery of generalised Gegenbauer functions of fractional degree, as an analysis tool  and a class of special functions with rich properties  \cite{Liu20JAT-Asymptotics}.

\subsection{Our contributions} Undoubtedly,  the Taylor formula  plays a foundational role in numerical  analysis and algorithm development.  We present  a new fractional Taylor formula for  AC-BV functions with fractional regularity index (see Theorem \ref{NewFTFV1}) that  ``interpolates"  and seamlessly  bridges the Taylor formulas of two consecutive integer orders.  From this tool, we can derive  an analogous formula for the Legendre expansion coefficient of the same class of functions, which
turns out  the cornerstone of all the analysis.
Then we obtain a set of optimal Legendre approximation results in $L^\infty$- and $L^2$-norms for functions with both interior and endpoint singularities.
We highlight that the use of  function space involving fractional integrals/derivatives to characterise regularity follows  that in  \cite{Liu19MC-Optimal}, but we further refine this framework by introducing the Caputo derivative. When the fractional regularity index takes  integer value, it reduces to the AC-BV space in Trefethen  \cite{Trefethen2008SIREV,Trefethen2013Book} (with adaption to the  Legendre approximation).  We  point out that  the argument for the Legendre approximation herein  is different from that for the Chebyshev approximation in   \cite{Liu19MC-Optimal}.   It is also noteworthy that  Babu\v{s}ka and Hakula \cite{Babuska19CMAME-Pointwise} discussed the point-wise error estimates
of  the Legendre expansion for the specific  function \eqref{uxy1} (which is also the subject of \cite{Gui1986NM-Part1}) including known and unknown results. In fact, it appears necessary to study the point-wise error in the Legendre or other Jacobi cases. For example, the estimating the $L^\infty$-error like the Chebyshev expansion can only lead to suboptimal results for
functions with the interior singularity, e.g., $u(x)=|x|,$ as  a loss of half order occurs.  It was observed numerically, but how to obtain optimal estimate appears open (see, e.g., \cite{Wang18AML-new}). Here, we shall provide an answer to this, and also to some conjectures in  \cite{Babuska19CMAME-Pointwise}.
We remark that we aim at deriving sharp and optimal estimates valid for all polynomial orders.  According to \cite{Babuska19CMAME-Pointwise},  in most applications the polynomial orders are relatively small compared to those in the asymptotic range, while  the existing theory does not address the behaviour of the pre-asymptotic error.  As a result, our arguments and results are  different from those in  \cite{Xiang20NM-Optimal},  where  some  asymptotic formulas were employed to derive Jacobi approximation of specific singular functions for large polynomial orders.
As a final remark,   this paper will be largely devoted to the $L^\infty$- and $L^2$-estimates of the  finite Legendre expansions, which lay the groundwork for establishing the approximation theory of other orthogonal projections, interpolations and quadrature for singular functions.   Indeed, these results can  enrich the theoretical foundation of   $p$ and $hp$ methods (cf. \cite{Funaro1992Book,Schwab1998Book,Bernardi1997Book,Canuto2006Book,Hesthaven2007Book,Shen2011Book}).

In a nutshell, the present study together with \cite{Liu20JAT-Asymptotics,Liu19MC-Optimal}  is far from being the last word on this subject.

The rest of the paper is organised as follows. In section \ref{sect:fintderG}, we derive the fractional Taylor formula for the AC-BV functions and present some  preliminaries to pave the way for all forthcoming  discussions. In section \ref{sect3main}, we obtain the main results on Legendre approximation of functions with interior singularities and extend the tools to study the
endpoint singularities in  section \ref{sect4main}.

\begin{comment}

In fact, we can  interpret the functions in those spaces being  of AC$_{m}$-BV$_{\!s,\theta}$-regularity  for some $s\in (0,1)$ and given $\theta\in (-1,1)$ corresponding to the location of the interior singularity (cf. Definition \ref{functionA}).
  For example,  the regularity indices $m,s$ of  \eqref{uxy1}  are $m+s=\alpha+1$ with $m=[\alpha+1].$ The ``fractional regularity"
  is characterised by the Riemann-Liouville (RL)  fractional integration of $u^{(m)}$ at $\theta$ of order $1-s$,  which can capture the additional regularity with the index  ``$s$",  as the RL  fractional integration can increase  the regularity of such singular functions. When $s=1,$ it reduces to the space  in \cite{Trefethen2008SIREV,Trefethen2013Book}.   When $\theta=\pm 1,$ it
leads to the framework for endpoint singularities (cf. Definition \ref{functionA-1}).
We  point out that the analytical study of Legendre expansion of \eqref{uxy1} in Gui and Babu\v{s}ka \cite{Gui1986NM-Part1}
actually sheds light on our analysis tools to  derive optimal error estimates.  More precisely, for $\theta\in (-1,1),$  the following exact formula for the Legendre coefficient was derived in \cite{Gui1986NM-Part1}:
\begin{equation}\label{unL-Gui}
\hat u_n^L=\frac{2n+1} 2 \int_{-1}^1 (x-\theta)^\alpha_+\, P_n(x)\,{\rm d}x=C_n^{(\alpha)} (1-\theta)^{\alpha+1} P_n^{(\alpha+1,-\alpha-1)}(\theta),
\end{equation}
where $C_n^{(\alpha)}$ is an explicit constant.

 The analysis therein heavily  relied on the  formula (cf. \cite[p. 833]{Gradshteyn1975Book}): {\em for $\mu>0$ and $\gamma\in (0,1),$
\begin{equation}\label{unL-Gui2}
\int_0^1 (1-y)^{\mu-1} P_n(1-\gamma y)\, {\rm d} y=\frac{\Gamma(\mu)\, n!} {\Gamma(\mu+n+1)} P_n^{(\mu, -\mu)}(1-\gamma),
\end{equation}
where $P_n^{(\mu, -\mu)}(x)$ is the {\rm(}generalised{\rm)} Jacobi polynomial with one parameter possibly $<-1.$} Interestingly,
with the substitution:  $t=1-\gamma y$ and $x=1-\gamma,$ the formula \eqref{unL-Gui2} turns out to be
\begin{equation}\label{unL-Gui3}
\frac{1} {\Gamma(\mu)}\int_x^1 \frac{P_n(t)} {(t-x)^{1-\mu}}\,{\rm d}t=\frac{n!} {\Gamma(\mu+n+1)} {(1-x)^{\mu} P_n^{(\mu, -\mu)}(x)},
\end{equation}
where the left-hand side is the RL fractional integral of $P_n(x),$ i.e.,  $I_{1-}^\mu P_n(x)$ (cf. \eqref{leftintRL}).

The above snapshots clearly indicate  that  the RL fractional calculus
and  (generalised) Jacobi polynomials  with the parameters in  non-classical  ranges
played an important part in the analysis of \cite{Gui1986NM-Part1}.   It is not a coincidence that  these threads can be
found  in the
optimal $L^\infty$-error estimates of Chebyshev approximation \cite{Liu19MC-Optimal}. In fact,
they can be seen more evidently  in this context, as
 the generalised  Gegenbauer  functions of fractional degree  (GGF-Fs) \cite{Liu19MC-Optimal,Liu20JAT-Asymptotics} (counterpart of generalised Jacobi polynomials) and  the RL fractional integration by parts are  indispensable for
 the derivation of the exact formula for  the Legendre expansion coefficient (see Theorems \ref{HatUnCaseC} and \ref{Thm:52}).
 As already mentioned, the analysis in \cite{Gui1986NM-Part1,Gui1986NM-Part2,Gui1986NM-Part3}  had deep impact upon
 the $hp$-approximation theory in Jacobi-weighted Besov  spaces \cite{Babuska2000NM,Babuska2001SINUM,Babuska2002MMMAS,Guo2013JAT-Direct}.  To some extent, we however drive this  to a different
 path with the aid of tools from RL fractional calculus, leading to the space of functions, e.g., with  AC$_{m}$-BV$_{\!s,\theta}$-regularity that can seamlessly bridge the gaps  in \cite{Trefethen2008SIREV,Trefethen2013Book}.
 \smallskip

We summarise and highlight that this work can stand out from  many existing works mentioned earlier,  in particular,
the very recent papers  \cite{Wang18AML-new,Liu19MC-Optimal,Babuska19CMAME-Pointwise,Xiang20NM-Optimal} attributed to  the following features.
\smallskip
\begin{itemize}
\item[(i)] Different from \cite{Wang18AML-new,Xiang20NM-Optimal} and the related works cited therein, we consider the Legendre approximation of singular functions with AC$_{m}$-BV$_{\!s,\theta}$-regularity, rather than the usual AC$_m$-BV-regularity \cite{Trefethen2008SIREV,Trefethen2013Book}.    Most existing results corresponding to  the special case with $s=1.$
Moreover, our argument leads to  sharper  results as we shall illustrate in the text.
\smallskip
\item[(ii)] \smallskip
\item[(iii)] We aim at deriving sharp and optimal estimates valid for all polynomial orders.  As pointed out in  \cite{Babuska19CMAME-Pointwise},  in most applications the polynomial orders are relatively small compared to those in the asymptotic range, while  the existing theory does not address the behaviour of the pre-asymptotic error.  As a result, our arguments and results are  different from those in  \cite{Xiang20NM-Optimal},  where  some  asymptotic formulas were employed to derive Jacobi approximation of specific singular functions for large polynomial orders.

\smallskip
\item[(iv)]  Some results in \cite{Liu19MC-Optimal,Liu20JAT-Asymptotics} pave  the way for this study,  but much more delicate analysis and additional tools are needed for the Legendre approximations.
In fact, it is  believed that the Legendre case has a more extendability of the arguments and results to the other Jacobi cases.
\end{itemize}
\smallskip

  \newpage
It is known that polynomial approximation theory is of fundamental importance in numerical analysis and  algorithm development of many computational methods, e.g.,     $p$/$hp$ finite elements  or spectral/spectral-element methods (see, e.g., \cite{Funaro1992Book,Schwab1998Book,Bernardi1997Book,Canuto2006Book,Hesthaven2007Book,Shen2011Book} and references therein).

Typically, the  documented approximation results   take  the form
\begin{equation}\label{existestm}
\| Q_N u-u \|_{\mathcal S_{\rm l}}\le cN^{-\sigma} |u|_{\mathcal B_{\rm r}},\quad \sigma\ge 0,
\end{equation}
where $Q_N$ is an orthogonal projection (or interpolation operator) upon the set of all polynomials of degree at most $N,$ and $c$ is a positive constant independent of $N$ and $u$. In \eqref{existestm},    ${\mathcal S}_{\rm l}$ is  a certain Sobolev space,
${\mathcal B}_{\rm r}$ is a related  Sobolev  or Besov space, and $\sigma$ depends on the regularity exponentials of both ${\mathcal B}_{\rm r}$ and ${\mathcal S}_{\rm l}$. In practice,  one would expect  {\em {\rm \bf (a)} the space ${\mathcal B}_{\rm r}$ should contain the classes of functions as broad as possible; {\rm and  {\bf  (b)}}  the space ${\mathcal B}_{\rm r}$ can best characterise their regularity   leading to  optimal order
of  convergence.}  In general,  the  space ${\mathcal B}_{\rm r}$ is of the following types.
\begin{itemize}
\item[(i)]  ${\mathcal B}_{\rm r}$ is the standard weighted Sobolev space $H^m_\omega(\Omega)$ with integer $m\ge 0$ and certain
weight function $\omega(x)$ on $\Omega=(-1,1)$ (see, e.g.,  \cite{Bernardi1997Book,Canuto2006Book,Hesthaven2007Book}).  However,  it  could not  lead to optimal order  for functions with endpoint singularities  (see, e.g., \cite{Bernardi1997Book,Guo2004JAT}) or with interior singularities, e.g., $|x|$ (see \cite{Trefethen2008SIREV}).

\vskip 4pt
\item[(ii)]  ${\mathcal B}_{\rm r}$ is the  non-uniformly Jacobi-weighted Sobolev space (see, e.g., \cite{Funaro1992Book,Babuska2000NM,Babuska2001SINUM,Guo2004JAT,Guo06JSC-Optimal,Shen2011Book}).  For example,  ${\mathcal B}_{\rm r} =H^{m,\beta}(\Omega)$ with integer $m\ge 0$ and $\beta>-1,$  is defined as  a closure of $C^\infty$-functions endowed with the weighted norm
\begin{equation}\label{weightednorm}
\|u\|_{H^{m,\beta}(\Omega)}=\bigg\{\sum_{k=0}^m \int_{-1}^1 |u^{(k)}(x)|^2(1-x^2)^{\beta+k} dx  \bigg\}^{1/2}.
\end{equation}
Compared with the standard Sobolev space in (i), such spaces can better  describe the  endpoint singularities, but
 still produce suboptimal estimates for $(1+x)^\alpha$-type  singular functions with non-integer $\alpha>0$   (cf. \cite[P. 474]{Castillo2002MC}).
Indeed, for the Chebyshev approximation, we
find that  $u=(1+x)^\alpha \in  {H^{m,-1/2}(\Omega)}$ with integer  $m<2\alpha+1/2,$   and
\begin{equation}\label{Jacoweighted00}
   \|\pi_N^C u- u\|_{L^2_\omega(\Omega)}\le c N^{-m}|u|_{H^{m,-1/2}(\Omega)},
   \end{equation}
    where $\pi_N^C u$ is the $L^2_\omega$-orthogonal projection of $u$ (with $\omega=(1-x^2)^{-1/2}$). However, the expected optimal order is  $O(N^{-2\alpha-1/2}),$ so the loss of an order of the fractional part of $2\alpha+1/2$ or one order (when $2\alpha=k+1/2$ with  integer $k\ge 0$),  is inevitable under this framework. This is due to the space $H^{m,\beta}(\Omega)$ is only defined for integer $m\ge 0.$

\vskip 4pt

\item[(iii)] In a series of works \cite{Babuska2000NM,Babuska2001SINUM,Babuska2002MMMAS},   Babu\v{s}ka and Guo  introduced the Jacobi-weighted Besov space
defined by space interpolation based on  the so-called K-method.
One commonly used  Besov space for $(1+x)^\alpha$-type corner singularities  is  ${\mathcal B}^{s,\beta}_{2,2}(\Omega)=(H^{l,\beta}(\Omega),H^{m,\beta}(\Omega))_{\vartheta, 2}$ with integers $l<m$ and  $s=(1-\vartheta) l+\vartheta m, \vartheta\in (0,1),$  equipped with the norm
\begin{equation}\label{weightednorm22}
\begin{split}
 & \|u\|_{{\mathcal B}^{s,\beta}_{2,2}(\Omega)}=\bigg(\int_{0}^\infty t^{-2\vartheta} |K(t,u)|^2 \frac{dt} t  \bigg)^{ 1/2}\!,\quad {\rm where}\\
& K(t,u)=\inf_{u=v+w} \big(\|v\|_{H^{l,\beta}(\Omega)}+t\|w\|_{H^{m,\beta}(\Omega)} \big).
\end{split}
\end{equation}

However, to deal with $(1+x)^\alpha \log^\nu (1+x)$-type corner singularities,  Babu\v{s}ka and Guo had to further modify  the K-method by incorporating  a log-factor into  the norm.

\end{itemize}

The aforementioned framework    might lead to  suboptimal estimates for  functions with interior singularities.
For example, we consider  $u(x)=|x|$ for $ x\in (-1,1).$ Note that   $u'(x)=2H(x)-1$  and $ u''(x)=2\delta(x)$ (where $H,\delta$ are respectively the Heaviside function and the Dirac delta function).   Since $u''\not \in L^2(\Omega)$,
the Chebyshev approximation of $|x|$ has a convergence:
   \begin{equation}\label{Jacoweighted}
   \|\pi_N^C u- u\|_{L^2_\omega(\Omega)}\le c N^{-1}|u|_{H^{1,-1/2}(\Omega)},
   \end{equation}
but   the expected  optimal order is  $O(N^{-3/2})$ (cf. \cite{Trefethen2008SIREV,Trefethen2013Book}). In fact,
as shown in  \cite[Thms 4.2-4.3]{Trefethen2008SIREV} and \cite[Thms 7.1-7.2]{Trefethen2013Book} (also see Lemma \ref{trenstha} below),  one should choose ${\mathcal B}_{\rm r}\subseteq {\rm BV}(\bar\Omega)$ (the space of functions of bounded variation)
to achieve optimality (see Section \ref{sect:existingest}, and  refer to  \cite{Trefethen2008SIREV,Xiang2010NM,Trefethen2013Book,Majidian2017ANM}  for more details).
Unfortunately, the Sobolev spaces therein  were defined through integer-order derivatives, so they could not best
characterise  the regularity of e.g., $u(x)=|x|^{\alpha}$ with non-integer $\alpha>0.$ In other words,  the order of convergence is  suboptimal.

In this paper, we intend  to introduce a new framework of fractional Sobolev-type spaces that can  meet the two  requirements {\bf (a)-(b)} and
overcome the deficiencies mentioned above.

We focus on the Chebyshev approximation but the  analysis techniques
are extendable to general Jacobi approximations.    Here, we  put the emphasises  on
estimating  the decay rate of  expansion coefficients for the reason that  the errors of spectral expansions in various norms, and the related interpolation and quadratures
  can be estimated directly from the sums of the coefficients (cf. \cite{Trefethen2008SIREV,Majidian2017ANM}).

The essential ideas and main contributions of this study  are summarised as follows.
 \begin{itemize}
\item[(i)]   We derive the exact representation of the Chebyshev expansion coefficients  (see Theorem \ref{IdentityForUn}) by using  the fractional calculus properties of  GGF-Fs and  fractional integration by parts (under the weakest possible conditions).
  This allows us to naturally define the fractional Sobolev spaces to characterise the regularity
 of a large class of singular functions, leading to optimal order of convergence.

\vskip 4pt
   \item [(ii)]  As a by-product,
   our estimates for the functions in this framework with integer regularity exponential  can  improve  the existing bounds in usual Sobolev spaces (see, e.g., \cite{Trefethen2008SIREV,Xiang2010NM,Trefethen2013Book,Majidian2017ANM}).

    \item[(iii)]  We present some useful analytical formulas on fractional calculus of GGF-Fs, and  the Chebyshev expansions of some specific singular functions.    Some of them  are new or difficult to be derived by other means  (cf. \cite{Gui1986NM,Boyd1989AMC,Wang18JCAM-convergence}).  They are  also  useful for the design of spectral algorithms.
 \end{itemize}
\vskip 4pt

 The paper is organised as follows.  In Sections \ref{sect:gegenbaure}-\ref{sect:fintderG}, we introduce the GGF-Fs, and present their important  properties, including the uniform bounds and Riemann-Liouville fractional integral/derivative formulas.  We derive the main results in Section \ref{mainsect:ms}, and    improve the existing estimates in Sobolev spaces with integer-order derivatives in Section \ref{sect:existingest}.   We discuss in Section \ref{Sect6Analysis} the extension of the main results to the analysis of interpolation, quadrature and endpoint singularities.
 \end{comment}

\begin{comment}
\section{Preliminaries and illustrative examples}

\subsection{Illustrative examples} For clarity of presentations, we start with the some relatively simple but illustrative examples to motivate the introduction of the analysis tools.

Given  $\theta\in [-1,1),$ we consider
\begin{equation}\label{modelexample}
u(x)=g(x)|x-\theta|^\alpha,\quad x\in (-1,1),
\end{equation}
where $g(x)\in C^\infty(-1,1).$

 More recently, Wang \cite{Wang18AML-new} showed that if $f,f',\ldots,f^{(m-1)}$ are absolutely continuous, $f^{(m)}$ has a bounded variation and $\int^1_{-1}|f^{(m+1)}(x)|(1-x^2)^{-1/4}dx<\infty$
\begin{equation}\label{ancoefWang}
|\hat u_n^L|\le \frac{2}{\sqrt{\pi(2n-2m-1)}}\frac{\Gamma(n-m+1/2)}{\Gamma(n+1/2)}\int^1_{-1}|f^{(m+1)}(x)|(1-x^2)^{-1/4}dx.
\end{equation}
The $L^2$ orthogonal projection is defined by
$$\pi_N^L  u(x):=\sum_{n=0}^{N}\hat u_n^L\, P_n(x).$$

The decay of the

\section{Preliminaries}\label{sect:fintderG}
\setcounter{equation}{0}
\setcounter{lmm}{0}
\setcounter{thm}{0}

of fractional  calculus, related  spaces  of functions, and   with some their properties.

\subsection{Fractional  calculus and related  spaces  of functions}
\end{comment}

\section{Fractional integral/derivative formulas of GGF-Fs}\label{sect:fintderG}
\setcounter{equation}{0}
\setcounter{lmm}{0}
\setcounter{thm}{0}

In this section, we make necessary preparations for the forthcoming   discussions. More precisely, we first introduce several spaces of functions that will be used to characterise the regularity of the class of functions of interest.  We then recall the definition of  the Riemann-Liouville (RL) fractional integrals, and  present   a useful RL fractional integration parts formula.  Finally, we collect some relevant properties of generalised Gegenbauer functions of fractional degree (GGF-Fs), which were first introduced and studied in
\cite{Liu19MC-Optimal,Liu20JAT-Asymptotics}.

\subsection{Spaces of functions}\label{Subsection2.1}

 Let $\Omega=(a,b)\subset {\mathbb R}$ be a finite open interval. For  real $p\in [1, \infty],$ let $L^p(\Omega)$ (resp. $W^{m,p}(\Omega)$ with $m\in {\mathbb N},$ the set of all positive integers)
 be the usual
$p$-Lebesgue space (resp. Sobolev space), equipped  with the norm $\|\cdot\|_{L^p(\Omega)}$ (resp. $\|\cdot\|_{W^{m,p}(\Omega)}$), as  in Adams \cite{Adams1975Book}.

 Let $C(\bar \Omega)$ be the classical space of continuous functions, and  ${\rm AC}(\bar \Omega)$
the space of  absolutely continuous functions on $\bar \Omega.$  It is known that
every absolutely continuous function is uniformly continuous (but the converse is not true), and hence continuous (cf. \!\cite[p.\! 483]{Ponnusamy2012Book}). It is known  that {\em a real function  $f(x)\in {\rm AC}(\bar \Omega)$ if and only if $f(x)\in L^1(\Omega)$, $f(x)$ has a derivative $f'(x)$
almost everywhere on $[a,b]$ such that   $f'(x)\in L^1(\Omega), $ and $f(x)$ has the integral representation:}
\begin{equation}\label{AVint}
f(x)=f(a)+\int_a^x f'(t)\,\rd t,\quad \forall\,  x\in [a,b],
\end{equation}
(cf.\! \cite[Chap.\! 1]{Samko1993Book} and  \cite[p.\! 285]{Lang1993Book}).

\begin{comment}
According to
\cite[Ch. 3]{Leoni2009Book} (also see \cite[P. 206]{Brezis2011Book} or  \cite[Ch. 1]{Samko1993Book}), a real function $f(x)\in {\rm AC}(\bar \Omega),$ if for any $\varepsilon>0,$ there exists $\delta>0,$ such that for every finite sequence of disjoint intervals $(a_k, b_k)\subset \Omega$ such that $\sum_k |b_k-a_k|<\delta,$ we have $\sum_k |f(b_k)-f(a_k)|<\varepsilon.$

Recall that (cf. \cite[Ch. 1]{Samko1993Book} or  \cite[P. 285]{Lang1993Book}):  a real function  $f(x)\in {\rm AC}(\bar \Omega)$ if and only if $f(x)\in L^1(\Omega)$, $f(x)$ has a derivative $f'(x)$
almost everywhere on $[a,b]$ such that   $f'(x)\in L^1(\Omega), $ and $f(x)$ has the integral representation:
\begin{equation*}\label{AVint}
f(x)=f(a)+\int_a^x f'(t)\,dt,\quad \forall x\in [a,b].
\end{equation*}
\end{comment}

Let ${\rm BV}(\bar \Omega)$ be the space of functions of bounded  variation on $[a,b].$  We say that  {\em a real function $f(x)\in {\rm BV}(\bar \Omega), $ if there exists a constant $C>0$ such that
$$V({\mathbb P};f):=\sum\limits_{i=0}^{k-1}|f(x_{i+1})-f(x_i)|\le C$$
 for every finite partition
${\mathbb P}=\{x_0, x_1,\cdots, x_k\}$ {\rm(}satisfying $x_i<x_{i+1}$ for all $0\le i\le k-1${\rm)} of $[a,b].$}
Then  the total variation of $f$ on $[a,b]$ is defined as $V_{\bar \Omega}[f]:=\sup\{V({\mathbb P};f)\},$ where the supreme is taken over all the partitions of $\bar \Omega$ (cf.
\cite[p.\! 207]{Brezis2011Book} or \cite[Chap.\! X]{Lang1993Book}).

An important characterisation of a BV-function  is the Jordan decomposition
(cf. \!\cite[Thm.\! 11.19]{Ponnusamy2012Book}):  {\em a function is of bounded variation if and only if it can be expressed as the difference of two increasing functions on $[a,b].$}
As a result,   every function in  ${\rm BV}(\bar\Omega)$ has at most a countable number of discontinuities, which are  either jump or removable discontinuities,
so it is differentiable almost everywhere.  Indeed, according to
 \cite[p.\! 223]{Appell2014Book},  {\em  if $f \in {\rm AC}(\bar \Omega),$ then}
\begin{equation*}\label{Vnorm}
V_{\bar \Omega}[f]= \int_{\Omega} |f'(x)|\, {\rm d}x.
\end{equation*}

In fact, we have  ${\rm BV}(\bar \Omega) \subset {\rm AC}(\bar \Omega)={\rm W}^{1,1}( \Omega)$  in the sense that   every $f(x) \in  {\rm AC}(\bar \Omega)$ has an almost everywhere classical derivative $f'\in L^{1}(\Omega)$ (cf.  \eqref{AVint}) and  $f'(x)$ is the weak derivative
of $f(x).$ Conversely, even $f \in W^{1.1}(\Omega),$ modulo a modification on a set of measure
zero, is an absolutely continuous function (cf.\! \cite[p.\! 206]{Brezis2011Book} and   \cite[p.\! 84; p.\! 96]{Buttazzo1998Book}).

For BV-functions,  we can define  the Riemann-Stieltjes (RS) integral  (cf.\! \cite[Chap.\!\! X]{Lang1993Book}).  A function $f(x)$ is said to be  RS($g$)-integrable, if $\int_\Omega f \rd g<\infty$ for  $g\in {\rm BV} (\bar \Omega).$
From  \cite[Prop.\! 1.3]{Lang1993Book}, we have  that
{\em  if $f$ is  {\rm RS($g$)-integrable}, then
\begin{equation}\label{integralProp}
\Big|\int_\Omega f(x)\, \rd g(x)\Big|\le \|f\|_\infty\, V_{\bar \Omega}[f],\quad  \int_\Omega  |\rd g(x)|=V_{\bar \Omega}[g],
\end{equation}
where $\|f\|_\infty$ is the $L^\infty$-norm of $f$ on $[a,b].$}

In the analysis, we shall also use the splitting rule of a RS integral, which is different from the usual integral.
\begin{lmm}[{see Carter and Brunt \cite[Thm 6.1.1 \& Thm 6.1.6]{Carter2000Book}}]\label{LmmCarter2000Book-1}
  If the interval $\Omega$ is a union of a finite number of pairwise disjoint intervals
$
\Omega=\Omega_{1} \cup \Omega_{2} \cup \cdots \cup \Omega_{m},
$
then
\begin{equation*}
\int_{\Omega} f \rd g=\sum_{j=1}^{m} \int_{\Omega_{j}} f \rd g
\end{equation*}
in the sense that if one side exists, then so does the other, and the two are equal. Moreover,   for any function $f$ defined at $\theta,$ then
 \begin{equation*}\label{BVtheta}
  \int_{[\theta, \theta]} f \rd  g =f(\theta)(g(\theta+)-g(\theta-)).
  \end{equation*}
\end{lmm}

\subsection{Formula of fractional integration by parts}
Recall the formula of integration by parts  involving
the Riemann-Stieltjes integrals  (cf. \!\cite[(1.20)]{Klebaner2005Book}):  {\em if  $f,g\in {\rm BV}(\bar\Omega)$, we have}
	\begin{align} \label{IPPW1122}
	\int_{a}^b f(x)  \, {\rm d} g(x) =\{f(x)g(x)\}\big|_{a^+}^{b^-} - \int_{a}^b  g(x)\,  {\rm d} f(x),
	\end{align}
where we denote
\begin{equation*}\label{uvpmlim}
	f(x)\big|_{a^+}^{b^-}=\lim\limits_{x\to b^-} f(x)- \lim\limits_{x\to a^+} f(x)=f(b-)-f(a+).
\end{equation*}
In particular, if $f, g \in \mathrm{AC}(\bar{\Omega}),$ we have
\begin{equation*}\label{uvpmlim00}
\int_{a}^{b} f(x) g^{\prime}(x)\, \rd x+\int_{a}^{b} f^{\prime}(x) g(x)\, \rd x=\{f(x) g(x)\}\big|_{a} ^{b}.
\end{equation*}
In what follows, we shall derive a formula of fractional integration parts from \eqref{IPPW1122} in a weaker  sense  than   the existing counterparts (cf.\! \cite{Samko1993Book,Bourdin2015ADE}).  For this purpose, we recap on  the definition of  the  Riemann-Liouville  fractional integral (cf.\!  \cite[p.\! 33, p.\! 44]{Samko1993Book}):  {\em  for any   $f\in L^1(\Omega),$  the left-sided and right-sided  Riemann-Liouville  fractional integrals of  real order $\rho \ge 0$ are defined  by}
	\begin{equation}\label{leftintRL}
	\begin{split}
	& (I_{a+}^\rho f) (x)=\frac 1 {\Gamma(\rho)}\int_{a}^x \frac{f(y)}{(x-y)^{1-\rho}} \rd y; \quad

	(I_{b-}^\rho f) (x)=\frac 1 {\Gamma(\rho)}\int_{x}^b \frac{f(y)}{(y-x)^{1-\rho}} \rd y,
	\end{split}
	\end{equation}
{\em for $x\in \Omega,$ where $\Gamma(\cdot)$ is the usual Gamma function. For $\mu \in (k-1, k]$ with $k \in \mathbb{N},$ the left-sided  and right-side Caputo fractional derivatives of order $\mu$ are respectively  defined by}
\begin{equation}\label{CaputoDev}
({ }^{C}\! D_{a+}^{\mu} f)(x)= (I_{a+}^{k-\mu} f^{(k)})(x); \quad ({ }^{C}\! D_{b-}^{\mu} f)(x)=(-1)^k (I_{b-}^{k-\mu} f^{(k)})(x).
\end{equation}

The following formula of fractional integration by parts plays an important role in the analysis, which can be derived from \eqref{IPPW1122}  (see  Appendix \ref{AppendixA}).
\begin{lmm}\label{FracIntPart-0} Let $\rho\ge 0, f(x)\in L^1(\Omega)$ and $ g(x)\in {\rm AC}(\bar \Omega).$

	\begin{itemize}
		\item[(i)] If $ I_{b-}^{\rho} f(x)\in {\rm BV}(\bar\Omega),$ then
	\begin{equation}\label{FracIntPart-1}
\begin{split}
\int^{b}_{a}   f(x) & \,I_{a+}^{\rho} g'(x)\,{\rm d}x  = \big\{g(x)\,  I_{b-}^{\rho} f(x)\big\}\big|_{a^+}^{b^-} - \int^{b}_{a} g(x)
\,{{\rm d}\big\{I_{b-}^{\rho} f(x)\big\}}.
\end{split}
\end{equation}
	\item[(ii)] If $ I_{a+}^{\rho} f(x)\in {\rm BV}(\bar\Omega),$ then
	\begin{equation}\label{FracIntPart-2}
\begin{split}
\int^{b}_{a}   f(x) & \,I_{b-}^{\rho} g'(x)\,{\rm d}x  = \big\{g(x)\,  I_{a+}^{\rho} f(x)\big\}\big|_{a^+}^{b^-} - \int^{b}_{a} g(x)
\,{{\rm d}\big\{I_{a+}^{\rho} f(x)\big\}}.
\end{split}
\end{equation}
	\end{itemize}
\end{lmm}
\begin{rem} {\em  If $\rho=0,$ then they reduce  \eqref{IPPW1122}.  It is known that the fractional integral can improve the regularity.  Indeed,  for $0<\rho<1$ and $u\in L^1(\Omega),$ we have $I_{a+}^{\rho} u, I_{b-}^{\rho} u\in L^p(\Omega)$ with $p\in [1, \rho^{-1})$ {\rm(cf.\! \cite[Prop.\! 2.1]{Bourdin2015ADE}).} } \qed
\end{rem}

Compared with those in \cite{Samko1993Book,Bourdin2015ADE},  a weaker condition is imposed on $f(x)$ in  \eqref{FracIntPart-1}-\eqref{FracIntPart-2}, which turns out essential in dealing with the singular functions.  Moreover, for such functions,
 the limit values   $\lim_{x\to a^+}  I_{a+}^\rho f(x)$ in \eqref{FracIntPart-1},  and   $\lim_{x\to b^-}  I_{b-}^\rho f(x)$ in \eqref{FracIntPart-2}  might  be nonzero,  in contrast to a usual integral with $\rho=1.$ For example,  for $\rho\in (0,1),$ we have
\begin{equation*}\label{intformu007}
 I_{a+}^{1-\rho} \, (x-a)^{\rho-1}=
 I_{b-}^{1-\rho} \, (b-x)^{\rho-1}=  \Gamma(\rho),
\end{equation*}
which follow from the explicit formulas  $($cf.  \cite{Samko1993Book}$)$:   for  real $\eta>-1$ and $\rho\ge 0,$
\begin{equation}\label{intformu}
 I_{a+}^\rho (x-a)^\eta=  \dfrac{\Gamma(\eta+1)}{\Gamma(\eta+\rho+1)} (x-a)^{\eta+\rho};
 \quad I_{b-}^\rho  (b-x)^\eta=  \dfrac{\Gamma(\eta+1)}{\Gamma(\eta+\rho+1)} (b-x)^{\eta+\rho}.

\end{equation}

In fact, we have the following more general formula,  which finds useful in  exemplifying some estimates to be presented later.
We  sketch the derivation in Appendix \ref{AppendixA0}.
\begin{prop}\label{boundaryvalue}  Let $f(x)=(x-a)^{\gamma}g(x)$ with real $\gamma>-1$, where $g(x)$  is bounded and Riemann integrable on $[a,a+\delta)$ for some $\delta>0.$  Then  for real  $\rho >0$,  we have
	\begin{equation}\label{integbA}
	\lim_{x\to a^+}(I_{a+}^\rho\, f) (x)=
	\begin{cases}
	0,  & {\rm if}\;\;  \rho >-\gamma, \\
	g(a)\Gamma(\gamma+1),&{\rm if}\;\; \rho=-\gamma,\\
	\infty, &{\rm if}\;\; \rho<-\gamma.
	\end{cases}
	\end{equation}
Let  $f(x)=(b-x)^{\gamma}g(x)$, $\gamma>-1$, and  $g(x)$ be bounded and Riemann integrable on $(b-\delta,b].$ Then the same result holds for the limit $\lim\limits_{x\to b^-}(I_{b-}^\rho\, f) (x)$ but with $g(b)$ in place of $g(a).$
 \end{prop}

 \subsection{Fractional Taylor formula} Needless to say, the Taylor formula plays a fundamental role in many branches of mathematics.
For comparison purpose, we recall this well-known formula:  {\em Let $k \geq 1$ be an integer and let $f(x)$ be a real function that is  $k$ times differentiable at the point $x=\theta$.  Further, let $f^{(k)}(x)$ be absolutely continuous on the closed interval between $\theta$ and $x$.  Then we have}
 \begin{equation}\label{TaylorA}
 f(x)= \sum_{j=0}^{k} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j}
 +
 \int_{\theta}^{x} \frac{f^{(k+1)}(t)}{k !}(x-t)^{k}\, \rd t.
 \end{equation}

 Note that since $f^{(k)}(x)$ is an AC-function, $f^{(k+1)}(x)$ exists as an $L^1$-function.

As a second building block for the analysis, we  derive a fractional  Taylor  formula from Lemma  \ref{FracIntPart-0} and
 \eqref{TaylorA}.
\begin{thm}[{\bf Fractional Taylor formula}]\label{NewFTFV1}  Let  $\mu \in (k-1, k]$ with $k \in \mathbb{N},$ and  let $f(x)$ be a real function that is  $(k-1)$times differentiable at the point $x=\theta$.
\begin{itemize}
 \item[(i)] If $f^{(k-1)}\in {\rm AC}([\theta, x])$ and    ${ }^{C}\! D_{\theta+}^{\mu} f \in  {\rm BV}([\theta, x]),$ then we have the left-sided fractional Taylor formula
\begin{equation}\begin{split}\label{LetfTaylor-00}
f(x)=\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j} + \frac{{ }^{C}\! D_{\theta+}^{\mu} f(\theta+)}{\Gamma(\mu+1)}{(x-\theta)^{\mu}}+  \frac{1}{\Gamma(\mu+1)}
 \int_{\theta}^{x} (x-t)^{\mu}\, \rd\{{ }^{C}\! D_{\theta+}^{\mu} f(t)\}.
\end{split}\end{equation}
\item[(ii)] If $f^{(k-1)}\in {\rm AC}([x,\theta])$ and    ${ }^{C}\! D_{\theta-}^{\mu} f \in  {\rm BV}([x,\theta]),$ then we have the right-sided fractional Taylor formula
\begin{equation}\begin{split}\label{LetfTaylor-000}
f(x)=\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j} + \frac{{ }^{C}\! D_{\theta-}^{\mu} f(\theta-)}{\Gamma(\mu+1)}{(\theta-x)^{\mu}}
-\frac{1}{\Gamma(\mu+1)} \int^{\theta}_{x} (t-x)^{\mu}\, \rd \{{ }^{C}\! D_{\theta-}^{\mu} f(t)\}.
\end{split}\end{equation}
\end{itemize}
\end{thm}
\begin{proof} By \eqref{TaylorA}  (with  $k\to k-1$), we have
\begin{equation}\begin{split}\label{LetfTaylor-1}
f(x)&=\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j}  +\frac{1}{(k-1)!} \int_{\theta}^{x}(x-t)^{k-1}f^{(k)}(t)\, \rd t.
\end{split}\end{equation}
From \eqref{intformu}, we find readily that  for $x>t,$
$$(x-t)^{k-1} =- \frac{(k-1)!}{\Gamma(\mu+1)}I_{x-}^{k-\mu}\Big\{\frac{\rd}{\rd t} (x-t)^{\mu}\Big\}.$$
Thus, we  can rewrite \eqref{LetfTaylor-1} as
\begin{equation}\begin{split}\label{LetfTaylor-10}
f(x)&=\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j}   - \frac{1}{\Gamma(\mu+1)} \int_{\theta}^{x}I_{x-}^{k-\mu}\Big\{\frac{\rd}{\rd t} (x-t)^{\mu}\Big\}f^{(k)}(t)\, \rd t.

\end{split}\end{equation}

  Substituting $a,$ $b,$ $\rho,$  $f$ and $g$ in \eqref{FracIntPart-2} of Lemma \ref{FracIntPart-0}  by $\theta,$ $x,$ $k-\mu$, $f^{(k)}(t)$ and $(x-t)^\mu,$ respectively, we obtain that  for $x>\theta,$
   \begin{equation}\begin{split}\label{LetfTaylor-2+1}
\int_{\theta}^{x}I_{x-}^{k-\mu}\Big\{\frac{\rd}{\rd t} (x-t)^{\mu}\Big\}f^{(k)}(t)\, \rd t
&=-{I_{\theta+}^{k-\mu} f^{(k)}(\theta+)}(x-\theta)^{\mu}-\int_{\theta}^{x} (x-t)^{\mu}\, \rd \{{I_{\theta+}^{k-\mu} f^{(k)}(t)}\} \\
&=-{ }^{C}\! D_{\theta+}^{\mu} f(\theta+){(x-\theta)^{\mu}}-\int_{\theta}^{x} (x-t)^{\mu}\, \rd \{{ }^{C}\! D_{\theta+}^{\mu} f(t)\},
  \end{split}\end{equation}
  where in the last step, we used the definition  \eqref{CaputoDev}.   Thus, we obtain  \eqref{LetfTaylor-00} from   \eqref{LetfTaylor-10}-\eqref{LetfTaylor-2+1} immediately.

  The right-sided formula \eqref{LetfTaylor-000} can be obtained in a very similar fashion.
\end{proof}

\begin{rem}\label{existing} {\em When $\mu=k,$ the fractional Taylor formulas  \eqref{LetfTaylor-00} and \eqref{LetfTaylor-000}
lead to  \eqref{TaylorA}.
The fractional formula can be viewed as the ``interpolation" of  the integer-order Taylor  formulas with the regularity indexes  $k-1$ and $k.$ Apparently,  the integer-order Taylor formula \eqref{TaylorA} is exact for all $f\in {\mathcal P_k}={\rm span}
\{(x-\theta)^j\,:\, 0\le j\le k\}.$ In the fractional case, the exactness of \eqref{LetfTaylor-00}  is for all $f\in {\mathcal P_{k-1}}\cup\{(x-\theta)^\mu\}$ {\rm(}i.e., the remainder vanishes{\rm).} We can verify this readily from \eqref{CaputoDev} and the fundamental formula:
$ { }^{C}\! D_{\theta+}^{\mu} \{(x-\theta)^\mu\}=\Gamma({\mu+1}).$ Note that the right-sided formula \eqref{LetfTaylor-000} is exact for all  $f\in {\mathcal P_{k-1}}\cup\{(\theta-x)^\mu\}.$} \qed
\end{rem}

We remark that there are several versions of fractional Taylor formulas for functions with different regularities.  For example,
Anastassiou \cite[(21)]{Anastassiou2009Chaos-On} stated the right-sided fractional Taylor formula:  for real $\mu\ge 1,$  let $k=[\mu]$ be its integer part, and assume that $f,f',\ldots,f^{(k-1)}\in {\rm AC}([x,\theta]).$ Then
\begin{equation*}\label{Anast}
f(x)=\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j !}(x-\theta)^{j}+\frac{1}{\Gamma(\mu)} \int_{x}^{\theta}(t-x)^{\mu-1} \, { }^{C}\!D_{\theta{-}}^{\mu} f(t)\, \rd t.
\end{equation*}
Kolwankar and Gangal   \cite{Kolwankar1998PRL-Local}
presented  some  local fractional Taylor
expansion with a  different fractional derivative in the remainder.

\begin{comment}
Denote the ordinary derivatives by  ${\mathcal D}=d/dx$ and  ${\mathcal D}^k=d^k/dx^k $ with integer  $k\ge 2$.
\begin{defn}\label{RLFintder} {\em
 For any   $u\in L^1(\Omega),$  the left-sided and right-sided  Riemann-Liouville  fractional integrals of order $s \in {\mathbb R}^+$ are defined  by
	\begin{equation}\label{leftintRL}
	\begin{split}
	& (I_{a+}^s u) (x)=\frac 1 {\Gamma(s)}\int_{a}^x \frac{u(y)}{(x-y)^{1-s}} dy; \quad

	(I_{b-}^s u) (x)=\frac 1 {\Gamma(s)}\int_{x}^b \frac{u(y)}{(y-x)^{1-s}} dy, \;\;\; x\in \Omega.
	\end{split}
	\end{equation}
	A function $u\in L^1(\Omega)$ is said to possess  a   left-sided {\rm(}resp. right-sided\,{\rm)} Riemann-Liouville fractional derivative
 ${\mathcal D}_{\!a+}^s u$ {\rm (}resp. $ {\mathcal D}_{b-}^s u${\rm)}  of
order $s \in (0,1),$ if  $I_{a+}^{1-s}\, u\in {\rm AC}(\bar\Omega) $ {\rm(}resp. $ I_{b-}^{1-s}\, u\in {\rm AC}(\bar\Omega)${\rm)}.  Moreover,   we have
\begin{equation}\label{reprezentacja20s01}
({\mathcal D}_{\!a+}^{s} u)(x)= {\mathcal D}\big\{I_{a+}^{1-s}\, u\big\}(x),\quad ({\mathcal D}_{b-}^{s} u)(x)= -{\mathcal D}\big\{ I_{b-}^{1-s}\, u\big\}(x),\quad x\in \Omega.
\end{equation}
   Similarly, for $s\in [k-1,k)$ with $k\in {\mathbb N},$ the  higher order  left-sided and right-sided   Riemann-Liouville fractional derivatives for $u\in L^1(\Omega)$
   satisfying  $I_{a+}^{1-s}\, u,  I_{b-}^{1-s}\, u \in {\rm AC}^k(\bar\Omega)$ {\rm(}i.e., the space of all $f(x)$ having continuous derivatives up to order $k-1$ on $\bar \Omega,$ and $f^{(k-1)}\in {\rm AC}(\bar\Omega)${\rm)} are defined  by
	\begin{equation}\label{left-fra-der-rl}
({\mathcal D}_{\!a+}^s u)(x)={\mathcal D}^k\big\{\, I_{a+}^{k-s}\, u\big\}(x);\quad
({\mathcal D}_{b-}^s u)(x)=(-1)^k{\mathcal D}^k\big\{\,I_{b-}^{k-s}\, u\big\}(x).
	\end{equation} }
	\end{defn}

\vskip 4pt

As a generalisation of  \eqref{AVint},  we have the following fractional integral representation,   which can also be regarded as the   definition of   Riemann-Liouville fractional derivatives alternative to Definition \ref{RLFintder} (see \cite[Prop.\!  5]{Bourdin2015ADE} and    \cite[P. 45]{Samko1993Book}).
\begin{prop}\label{fractrep}
 A function $u\in L^1(\Omega)$  possesses  a left-sided Riemann-Liouville fractional derivative
 ${\mathcal D}_{\!a+}^s u$  of
order $s \in (0,1),$ if and only if there exist  $C_a\in {\mathbb R}$ and  $\phi \in L^1(\Omega)$ such that
\begin{equation}\label{reprezentacja}
u(x)= \frac{C_a}{\Gamma(s)}(x-a)^{s-1} +  (I_{a+}^s\,\phi)(x) \;\;  \text{a.e. on}\;\;  [a,b],
\end{equation}
where $C_a=( I_{a+}^{1-s}u)(a)$ and $\phi(x) = ({\mathcal D}_{\!a+}^s u)(x)$ a.e. on $ [a,b].$

Similarly, a function
$u\in L^1(\Omega)$  has a right-sided Riemann-Liouville fractional derivative
 ${\mathcal D}_{\!b-}^s u$  of
order $s \in (0,1),$ if and only if there exist  $C_b\in {\mathbb R}$ and  $\psi\in L^1(\Omega)$ such that
\begin{equation}\label{reprezentacja2}
u(x)= \frac{C_b}{\Gamma(s)}(b-x)^{s-1} +  (I_{b-}^s\,\psi)(x) \;\;   \text{a.e. on}\;\;  [a,b],
\end{equation}
where $C_b=( I_{b-}^{1-s}u)(b)$ and $\psi(x) = ({\mathcal D}_{b-}^s u)(x) $ a.e. on $ [a,b].$
\end{prop}
	\end{comment}

\begin{comment}
\subsection{Fractional Sobolev-type spaces}

We define the fractional Sobolev space to characterize functions with interior or endpoint singularities (cf. \cite{Liu19MC-Optimal}).
For  a fixed   $\theta\in \Omega=(-1,1), $ we denote $\Omega_\theta^-:=(-1,\theta)$ and  $\Omega_\theta^+:=(\theta,1).$  For   $m\in {\mathbb N}_0$ and {$s\in (0,1),$}
we define the fractional Sobolev-type space:
	\begin{equation}\label{fracA}
\begin{split}
{\mathbb  W}^{m+s}_{\theta}(\Omega):= \big\{ &u\in L^1(\Omega)\, :\,   u, u',\cdots, u^{(m-1)}\in {\rm AC}(\bar \Omega)  \;\; {\rm and}\\
&  I_{\theta-}^{1-s} u^{(m)}\in  {\rm BV}(\bar\Omega_\theta^-),\quad  I_{\theta+}^{1-s} u^{(m)}\in  {\rm BV}(\bar\Omega_\theta^+)\big\},
	\end{split}
\end{equation}
equipped with the norm:
\begin{equation}\label{Wmsnorm}
\|u\|_{{\mathbb  W}^{m+s}_{\theta}(\Omega)}=\sum_{k=0}^m \|u^{(k)}\|_{L^1(\Omega)} +   U^{m,s}_\theta,
\end{equation}
where the semi-norm is defined by
\begin{itemize}
\item for $m=1,2,\cdots,$ and  {$s\in (0,1),$}
\begin{equation}\label{seminormF}
  \begin{split}
  U^{m,s}_\theta:=&\int_{-1}^\theta  \big | {\mathcal D}_{\!\theta-}^{s} \, u^{(m)}(x) \big|\,  dx
  +\int^{1}_\theta    |{\mathcal D}_{\!\theta+}^{s}  u^{(m)}(x)| dx\\
  &+\big|\big\{  I_{\theta-}^{1-s} u^{(m)}\big\}(\theta+)\big|+\big|\big\{  I_{\theta+}^{1-s} u^{(m)}\big\}(\theta-)\big|;
  \end{split}
  \end{equation}
  \item for  $m=0$ and {$s\in (1/2,1),$}
      \begin{equation}\label{Pnorm00}
  \begin{split}
 U^{0,s}_\theta:= & \int_{-1}^\theta  \big | {\mathcal D}_{\!\theta-}^{s} \, u(x) \big|\,(1-x^2)^{s/2-1/2} \rd x
  +\int^{1}_\theta    |{\mathcal D}_{\!\theta+}^{s}  u(x)| \, (1-x^2)^{s/2-1/2} \rd x\\
  &\quad\quad\quad+\big|\big\{(1-x^2)^{s/2-1/2}\,  I_{\theta-}^{1-s} u\big\}(\theta+)\big|+\big|\big\{(1-x^2)^{s/2-1/2}\,  I_{\theta+}^{1-s} u\big\}(\theta-)\big|.
  \end{split}
  \end{equation}
 \end{itemize}
{In analogue to the integer derivative  case  (see \cite{Trefethen2008SIREV}),
 the integrals in \eqref{seminormF}-\eqref{Pnorm00} are in the  Riemann-Stieltjes sense.
 For example, for $I_{\theta-}^{1-s} u^{(m)}\in  {\rm BV}(\bar\Omega_\theta^-),$ we understand
 $$ \int_{-1}^\theta  \big | {\mathcal D}_{\!\theta-}^{s} \, u^{(m)}(x) \big|\,  dx=\int_{-1}^\theta  \big | d
(I_{\theta-}^{1-s} u^{(m)}) \big|=V_{g}[-1,\theta],\quad g=I_{\theta-}^{1-s} u^{(m)}, $$
where ${\mathcal D}_{\!\theta-}^{s} \, u^{(m)}(x) $ is a Radon measure.

By virtue of   \eqref{integralProp}, it can be bounded by the total variation of $I_{\theta-}^{1-s} u^{(m)}$ on $\bar \Omega_\theta^-.$}

 \begin{rem}\label{propremk}
{ \em
 The parameter $\theta$ is related to the location of the singular point of $u(x)$. For example, if $u=|x|,$ then $\theta=0.$
For a function of multiple interior  singular points, we partition $(-1,1)$ into multiple subintervals and introduce the same number of parameters accordingly.
  }

  Accordingly, we define the corresponding space as
\begin{equation}\label{intergerASobB}
\begin{split}
{\mathbb  W}^{m+1}(\Omega):= \big\{ &u\in L^1(\Omega)\, :\,    u',\cdots, u^{(m-1)}\in {\rm AC}(\bar\Omega), \; u^{(m)}\in  {\rm BV}(\bar\Omega)\big\},
	\end{split}
\end{equation}
as in \cite{Trefethen2013Book,Majidian2017ANM}.

\end{rem}

{To deal with endpoint singularities,  we define the  fractional Sobolev-type spaces corresponding to $\theta=\pm 1$  by}
\begin{equation}\label{resultpm1}
\begin{split}
& {\mathbb  W}^{m+s}_{1-}(\Omega):= \big\{ u\in L^1(\Omega)\, :\,   u, u',\cdots, u^{(m-1)}\in {\rm AC}(\bar\Omega),\;  I_{1-}^{1-s} u^{(m)}\in  {\rm BV}(\bar\Omega)\big\},\\
& {\mathbb  W}^{m+s}_{-1+}(\Omega):= \big\{ u\in L^1(\Omega)\, :\,   u, u',\cdots, u^{(m-1)}\in {\rm AC}(\bar\Omega),\;  I_{-1+}^{1-s} u^{(m)}\in  {\rm BV}(\bar\Omega)\big\}.
\end{split}
\end{equation}

{Accordingly, the semi-norm
  $U^{m,s}_+$ (resp.    $U^{m,s}_-$)  only involves the right  (resp.  left) Riemann-Liouville fractional integrals/derivatives. For example,  when  $\theta=-1,$  the semi-norm corresponding to \eqref{seminormF} becomes
  \begin{equation}\label{seminormF--}
  \begin{split}
  U^{m,s}_{-1+}:=
  \int^{1}_{-1}    |{\mathcal D}_{\!-1+}^{s}  u^{(m)}(x)| dx+\big|\big\{  I_{-1+}^{1-s} u^{(m)}\big\}(-1+)\big|.
  \end{split}
  \end{equation}
    }
\end{comment}

\begin{comment}
\begin{theorem}\label{BoundGegC}
	For real  $ \lambda >1$, we have the bound: for $x^2\le 1-{\lambda(\lambda-1)}/{(\lambda+\nu)^2}$,
	\begin{equation}\label{BoundGegC-0}
	\max\big\{	(1-x^2)^{\lambda/2} \big|{}^{r\!}G_\nu^{(\lambda)}(x)\big|,	\; (1-x^2)^{\lambda/2}\big| {}^{l}G_\nu^{(\lambda)}(x)\big| \big\}\le
	\widehat	\kappa_\nu^{(\lambda)}
	,\;\;\;  \nu\ge 0,
	\end{equation}
	where
	\begin{equation}\label{kappaC}
	\begin{split}
	\widehat\kappa_\nu^{(\lambda)}=\frac{\Gamma(\lambda+1/2)}{\sqrt \pi}\bigg(
	\frac{ \cos^2(\pi\nu/ 2)\Gamma^2(\nu/2+1/2) }{\Gamma^2( ({\nu}+1)/2+\lambda)}+\frac{4\sin ^2\big(\pi  {\nu} /2\big)}{\nu^2+2\lambda \nu+\lambda}\frac{\Gamma^2(\nu/2+1)}{\Gamma^2({\nu}/2+\lambda)}\bigg)^{\frac 12}.
	\end{split}	\end{equation}
\end{theorem}
\begin{proof} {\color{red} To   be revised in the earlier paper!}
	Thanks to  \eqref{obsvers},  it suffices to  prove the result for  $ {}^{r\!}G_\nu^{(\lambda)}(x).$

We denote
	\begin{equation}\label{notationGMHB}
	\begin{split}
	&  \widehat M(x):= (1-x^2)^{\lambda/2} {}^{r\!}G_\nu^{(\lambda)}(x); \quad \widehat H(x):=\widehat M^2(x)+ \frac 1 {\rho(x)}\big(\widehat M'(x)\big)^2,\\
	& \rho(x):=\big((\nu+\lambda)^2(1-x^2)-\lambda (\lambda-1)\big)(1-x^2)^{-2}.
	\end{split}
	\end{equation}
	Note that by Proposition \ref{nonvember},  $\widehat H(x)$ is continuous on $[-1,1].$

	Taking $a=-\nu, b=\nu+2\lambda,c=\lambda +1/2$ and $z=(1\pm x)/2$ in  \eqref{SLPHF},  we find that
	$G(x)$ satisfies the Sturm-Liouville problem
	\begin{equation}\label{SLPGGF-F}
	\begin{split}
	\big(\omega_{\lambda+1}(x) G'(x)\big)'+\nu(\nu+2\lambda)\omega_{\lambda}(x) G(x)=0,\quad x\in (-1,1).
	\end{split}
	\end{equation}
	A direct calculation from   \eqref{SLPGGF-F} leads to
	\begin{equation}\label{BoundGegB-1}
	(1-x^2)\widehat M''(x)-x \widehat M'(x)+(1-x^2)\rho(x) \widehat M(x)=0,  \quad x\in (-1,1).
	\end{equation}
	Thus,  we have
	\begin{equation}\label{BoundGegB-2}
	\begin{split}
	\widehat H'(x)

	 &=\frac{2\lambda(\lambda-1)x}{(\lambda+\nu)^2(1-x^2)^2-\lambda(\lambda-1)(1-x^2)}\big(\widehat M'(x)\big)^2,\quad x\in (-1,1).
	\end{split}
	\end{equation}
	Hence, for $x^2\le 1-{\lambda(\lambda-1)}/{(\lambda+\nu)^2},$  $	\widehat H(x)$ is increasing for $-\sqrt{1-{\lambda(\lambda-1)}/{(\lambda+\nu)^2}}\le x<0,$ and decreasing  for $0<x\le \sqrt{1-{\lambda(\lambda-1)}/{(\lambda+\nu)^2}},$  so we have
	\begin{equation}\label{BoundGegB-3}
	\widehat M ^2(x)\le\widehat	H(x)\le \widehat	H(0)=\widehat M^2(0)+ {\rho^{-1}(0)}\big(\widehat M'(0)\big)^2,\quad \forall\, x\in [-1,1].
	\end{equation}
		By  \eqref{nonitA} and \eqref{Fatzero},  we obtain that for $\lambda\ge 0,$
	\begin{equation}\label{LBoundGeg-5}
	\begin{split}
	M(0)=&{}^{r\!}G_\nu^{(\lambda)}(0)=
	{}_2F_1\Big(\!\!-\nu, \nu+2\lambda;\lambda+\frac 1 2;\frac{1} 2\Big)=\frac{\sqrt{\pi}\, \Gamma(\lambda+1/2)}{\Gamma(- {\nu}/2+
		1/2)\Gamma({\nu}/2+\lambda+ 1/2)}\\
	=&\sin \big(\pi ( {\nu} +1)/2\big)\frac{\Gamma(\lambda+1/2)\Gamma( {\nu}/2+1/2)}{\sqrt{\pi}\, \Gamma( {\nu}/2+\lambda+
		1/2)}.
	\end{split}\end{equation}
	By \eqref{LBoundGeg-5}, we have
	\begin{equation}\label{BoundGegB-4}
	\begin{split}
	\widehat M(0)
	=&\cos \big(\pi {\nu}/2\big)\frac{\Gamma(\lambda+1/2)\Gamma( {\nu}/2+1/2)}{\sqrt{\pi}\, \Gamma( {\nu}/2+\lambda+
		1/2)}.
	\end{split}\end{equation}
	Recall the identity (cf. \cite[(15.5.1)]{Olver2010Book}):
	\begin{equation}\label{BoundGegB-5}\frac{d}{dx} {}_2F_1(a, b;c;z)=\frac{ab}{c} {}_2F_1(a+1, b+1;c+1;z).
	\end{equation}
	By \eqref{rgjfdef} and  \eqref{BoundGegB-5}, we obtain
	\begin{equation}\label{BoundGegB-6}
	\begin{split}
	\frac{d}{dx}{}^{r\!}G_\nu^{(\lambda)}(x)
	=&\frac{d}{dx} {}_2F_1\Big(\!\!-\nu, \nu+2\lambda;\lambda+\frac 1 2;\frac{1-x} 2\Big)\\
	=&	\frac{\nu(\nu+2\lambda)}{2\lambda+1} {}_2F_1\Big(\!\!-\nu+1, \nu+2\lambda+1;\lambda+\frac 3 2;\frac{1-x} 2\Big).
	\end{split}\end{equation}
	If $\nu> 1/2$, using \eqref{BoundGegB-6}, we have
	\begin{equation}\label{BoundGegB-6+}
	\begin{split}
	\frac{d}{dx}{}^{r\!}G_\nu^{(\lambda)}(x)
	=	\frac{\nu(\nu+2\lambda)}{2\lambda+1}{}^{r\!}G_{\nu-1}^{(\lambda+1)}(x).
	\end{split}\end{equation}
	Thanks to
	$$\widehat M'(x)= -\lambda x(1-x^2)^{\lambda/2-1}\, {}^{r\!}G_\nu^{(\lambda)}(x)+(1-x^2)^{\lambda/2}\frac{d}{dx}{}^{r\!}G_\nu^{(\lambda)}(x),$$
	with \eqref{nonitA}, \eqref{Fatzero} and  \eqref{BoundGegB-6}, we obtain
	\begin{equation}\label{BoundGegB-7}
	\begin{split}
	\{\rho^{-1/2}(x)\widehat M'(x)\}|_{x=0}
	=&\frac{2\sin(\pi\nu/2 )}{\sqrt{\pi}}\frac{\Gamma(\nu/2+1)\,\Gamma(\lambda+1/2)}{\sqrt{\nu^2+2\lambda \nu+\lambda}\,\, \Gamma(\lambda+\nu/2)}.
	\end{split}\end{equation}
	From  \eqref{BoundGegB-3}, 	\eqref{BoundGegB-4} and \eqref{BoundGegB-7}, we derive
	\eqref{BoundGegC-0}.
\end{proof}

\end{comment}

\section{Legendre expansions of functions with interior singularities}\label{sect3main}
\setcounter{equation}{0}
\setcounter{lmm}{0}
\setcounter{thm}{0}

It is known that much of the error analysis for orthogonal polynomial approximation and associated interpolation and quadrature  relies on  the decay rate of the expansion coefficient (cf.\! \cite{Xiang2012SINUM,Majidian2017ANM}).
Remarkably, we find that  the  spirit in  deriving the fractional Taylor formula in Theorem \ref{NewFTFV1}  can be extended to obtain  an analogous  formula for  the Legendre expansion coefficient
\begin{equation}\label{ancoef}
\hat u_n^L=\frac{2n+1}{2}\int_{-1}^{1}{u(x)P_n(x)}\, {\rm d} x,
\end{equation}
where $P_n(x)$ is the Legendre polynomial of degree $n.$
This formula lays the groundwork for all the forthcoming analysis. In fact, the argument  is also different from that for  the Chebyshev expansion coefficient  in \cite{Trefethen2008SIREV,Trefethen2013Book,Majidian2017ANM,Liu19MC-Optimal}.

\subsection{Fractional formula for the Legendre expansion coefficient}\label{sub31}

In what follows,  we assume that $u$ has a limited regularity with an interior singularity at $\theta\in (-1,1),$  e.g.,
$u(x)=|x-\theta|^\alpha$ with $\alpha>-1.$   Note   that
the results   can be extended to multiple interior singularities straightforwardly.
\begin{theorem}\label{IdentityForUn}
Let  $\mu \in (k-1, k]$ with $k \in \mathbb{N}$ and let $\theta\in(-1,1).$  If  $u,u',\ldots,u^{(k-1)}\in {\rm AC}([-1,1]),$  ${ }^{C}\! D_{\theta+}^{\mu} u \in  {\rm BV}([\theta, 1])$ and
  ${ }^{C}\! D_{\theta-}^{\mu} u \in  {\rm BV}([-1,\theta]),$  then  we have the following  representation of the Legendre expansion coefficient for each $n\ge\mu+1,$
\begin{equation}\label{HatUnCaseC}
\begin{split}
\hat u_n^L = \frac{2n+1}{2} &\Big\{ (I_{1-}^{\mu+1} P_n)(\theta) ({{ }^{C}\! D_{\theta+}^{\mu} u)(\theta+)}
+\int_{\theta}^{1} (I_{1-}^{\mu+1} P_n)(x) \,  \rd \big\{{ }^{C}\! D_{\theta+}^{\mu} u(x)\big\}\\
&\quad  +(I_{-1+}^{\mu+1} P_n)(\theta) ({{ }^{C}\! D_{\theta-}^{\mu} u)(\theta-)}  - \int^{\theta}_{-1}(I_{-1+}^{\mu+1} P_n)(x)\, \rd \big\{{ }^{C}\! D_{\theta-}^{\mu} u(x)\big\}\Big\},

\end{split}
\end{equation}
where the fractional integrals of $P_n(x)$ can be evaluated explicitly  by
\begin{equation}\label{twoGGFs}
\begin{split}
&(I_{1-}^{\mu+1} P_n)(x)=
\frac{(1-x)^{\mu+1}}{\Gamma(\mu+2)} \frac{P_{n}^{(\mu+1, -\mu-1)}(x)}{P_{n}^{(\mu+1,-\mu-1)}(1)},\\
& (I_{-1+}^{\mu+1} P_n)(x)
=\frac{(1+x)^{\mu+1}}{\Gamma(\mu+2)} \frac{P_{n}^{( -\mu-1,\mu+1)}(x)}{P_{n}^{(\mu+1,-\mu-1)}(1)}.

\end{split}

\end{equation}
Here $P_{n}^{(\mu+1, -\mu-1)}(x)$ and $P_{n}^{( -\mu-1,\mu+1)}(x)$ are the generalised Jacobi polynomials defined by the hypergeometric functions  as in  Szeg\"o \cite[p.\! 64]{Szego1975Book}.

\end{theorem}
\begin{proof} Given the regularity of $u$, we  obtain from  the fractional Taylor formulas in Theorem \ref{NewFTFV1} that
for $x\in (\theta,1),$
\begin{equation}\begin{split}\label{IdentityForUn-1}
u(x)=\sum_{j=0}^{k-1} \frac{u^{(j)}(\theta)}{j!}(x-\theta)^{j} + \frac{{ }^{C}\! D_{\theta+}^{\mu} u(\theta+)}{\Gamma(\mu+1)}{(x-\theta)^{\mu}}+  \frac{1}{\Gamma(\mu+1)}
 \int_{\theta}^{x} (x-t)^{\mu}\, \rd \big\{{ }^{C}\! D_{\theta+}^{\mu} u(t)\big\},
\end{split}\end{equation}
and for $x\in (-1,\theta),$
\begin{equation}\begin{split}\label{IdentityForUn-2}
u(x)=\sum_{j=0}^{k-1} \frac{u^{(j)}(\theta)}{j!}(x-\theta)^{j} + \frac{{ }^{C}\! D_{\theta-}^{\mu} u(\theta-)}{\Gamma(\mu+1)}{(\theta-x)^{\mu}}
-\frac{1}{\Gamma(\mu+1)} \int^{\theta}_{x} (t-x)^{\mu}\, \rd \big\{{ }^{C}\! D_{\theta-}^{\mu} u(t)\big\}.
\end{split}\end{equation}
Substituting  \eqref{IdentityForUn-1} and  \eqref{IdentityForUn-2} into \eqref{ancoef} leads to
\begin{equation}\label{IdentityForUn-3}
\begin{split}
&\frac {2\,  \hat u_n^L}{2n+1}  =  \int_{-1}^{1}{u(x)P_n(x)}\, {\rm d} x= \sum_{j=0}^{k-1} \frac{u^{(j)}(\theta)}{j!}
\int_{-1}^{1}(x-\theta)^{j} {P_n(x)}\, {\rm d} x \\
&\;\;  +\frac {{ }^{C}\! D_{\theta+}^{\mu} u(\theta+)}{\Gamma(\mu+1)}\int_{\theta}^{1} {(x-\theta)^{\mu}}\,P_n(x)\, {\rm d} x
+\frac{1}{\Gamma(\mu+1)} \int_{\theta}^{1}\Big(\int_{\theta}^{x} (x-t)^{\mu} \rd \{{ }^{C}\! D_{\theta+}^{\mu} u(t)\}\Big)\,P_n(x)\, {\rm d} x\\
&\;\; +\frac{{ }^{C}\! D_{\theta-}^{\mu} u(\theta-)} {\Gamma(\mu+1)}  \int_{-1}^{\theta} {(\theta-x)^{\mu}}\,P_n(x)\, {\rm d} x-\frac{1}{\Gamma(\mu+1)} \int_{-1}^{\theta}\Big(\int_{x}^{\theta} (t-x)^{\mu} \rd \{{ }^{C}\! D_{\theta-}^{\mu} u(t)\}\Big)\,P_n(x)\, {\rm d} x.
\end{split}\end{equation}
From the orthogonality of the Legendre polynomials, we obtain that  for $n\ge \mu+1\ge k,$
\begin{equation}\label{IdentityForUn-4}
\begin{split}
\int_{-1}^{1}(x-\theta)^{j} P_n(x) {\rm d} x=0,\quad 0\le j\le k-1.
\end{split}\end{equation}
We find readily that for a fixed $\theta\in (-1,1),$
\begin{equation}\label{IdentityForUn-60}
\begin{split}
&\int_{\theta}^{1}\Big(\int_{\theta}^{x} {(x-t)^{\mu}}\, \rd \{{ }^{C}\! D_{\theta+}^{\mu} u(t)\}\Big)\,P_n(x)\, {\rm d} x=
\int_{\theta}^{1}\Big(\int_{t}^{1} {(x-t)^{\mu}} \,P_n(x) {\rm d} x\Big) \rd \{{ }^{C}\! D_{\theta+}^{\mu} u(t)\},
\end{split}\end{equation}
and
\begin{equation}\label{IdentityForUn-602}
\begin{split}
& \int_{-1}^{\theta}\Big(\int_{x}^{\theta} (t-x)^{\mu} \rd \{{ }^{C}\! D_{\theta-}^{\mu} u(t)\}\Big)\,P_n(x)\, {\rm d} x=
\int^{\theta}_{-1}\Big(\int_{-1}^{t} {(t-x)^{\mu}} \,P_n(x) \,{\rm d} x\Big)\, \rd \{{ }^{C}\! D_{\theta-}^{\mu} u(t)\}.
\end{split}\end{equation}
In view of the definition of the fractional integral in \eqref{leftintRL} and \eqref{IdentityForUn-4}-\eqref{IdentityForUn-602}, we can rewrite
\eqref{IdentityForUn-3} as
\begin{equation}\label{IdentityForUn-3A}
\begin{split}
\frac {2\,  \hat u_n^L}{2n+1}  &=
(I_{1-}^{\mu+1} P_n)(\theta)  ({{ }^{C}\! D_{\theta+}^{\mu} u)(\theta+)}
+\int_{\theta}^{1} (I_{1-}^{\mu+1} P_n)(t) \,  \rd \{{ }^{C}\! D_{\theta+}^{\mu} u(t)\}\\
&\quad  +(I_{-1+}^{\mu+1} P_n)(\theta) ({ }^{C}\! D_{\theta-}^{\mu} u)(\theta-) - \int^{\theta}_{-1}(I_{-1+}^{\mu+1} P_n)(t)\, \rd \{{ }^{C}\! D_{\theta-}^{\mu} u(t)\},
\end{split}\end{equation}
which yields \eqref{HatUnCaseC}.  The two fractional integral identities of $P_n(x)$
in \eqref{twoGGFs} can be obtained from the  formulas  of the Jacobi polynomials (cf.\! Szeg\"o  \cite[p.\! 96]{Szego1975Book}), due to the Bateman's fractional integral formula (cf.\! \cite{Andrews1999Book}).  This ends the proof.
\end{proof}

We see from the above proof that the identity \eqref{HatUnCaseC}  is rooted in the fractional Taylor formula in Theorem \ref{NewFTFV1}. Also note that
 when $\mu=k$, the formula  \eqref{IdentityForUn} takes a much simpler form. Firstly,  the AC-BV regularity reduces to the setting considered by Trefethen \cite{Trefethen2008SIREV,Trefethen2013Book},  Xiang and Bornemann \cite{Xiang2012SINUM} among others (where one motivative example for the framework therein is to best characterise the regularity of  $u(x)=|x|$).  Secondly,  from
  Szeg\"o  \cite[Chap.\!\! 4]{Szego1975Book}, we find that for $\mu>-2,n\ge 0,$
  \begin{equation}\label{Pnmumu}
  P_{n}^{(\mu+1,-\mu-1)}(1)=\frac{\Gamma(n+\mu+2)} {n!\,\Gamma(\mu+2)},
  \end{equation}
  and for $n\ge k+1,$
  \begin{equation}\label{Pnmumu2}
  P_{n}^{(-k-1,k+1)}(x)= \frac{(n-k-1)! (n+k+1)!} {(n!)^2}\Big(\frac{x-1} 2\Big)^{k+1} P_{n-k-1}^{(k+1,k+1)}(x).
  \end{equation}
  Thus,  we can rewrite the second formula in \eqref{twoGGFs} with $\mu=k$  in terms of the usual Jacobi polynomial as follows
  \begin{equation}\label{Pnmumu3}
(I_{-1+}^{k+1} P_n)(x) = \frac{(-1)^{k+1}\,(n-k-1)!} {2^{k+1}\,n!} (1-x^2)^{k+1} P_{n-k-1}^{(k+1,k+1)}(x).
  \end{equation}
Following the same lines as above and using the parity of Jacobi polynomials, we can  reformulate the first formula  in \eqref{twoGGFs} with $\mu=k$ as
 \begin{equation}\label{Pnmumu4}
(I_{1-}^{k+1} P_n)(x) =  \frac{(n-k-1)!} {2^{k+1}\,n!} (1-x^2)^{k+1} P_{n-k-1}^{(k+1,k+1)}(x)= (-1)^{k+1}\, (I_{-1+}^{k+1} P_n)(x).
  \end{equation}
 In view of this relation, we find from    \eqref{CaputoDev} with $\mu=k$ that

 \eqref{HatUnCaseC} reduces to
  \begin{equation}\label{HatUnCaseC-V2+1}
\begin{split}
\hat u_n^L & = \frac{2n+1}{2} \Big\{
  u^{(k)}(\theta+) (I_{1-}^{k+1} P_n)(\theta)
+\int_{\theta}^{1} (I_{1-}^{k+1} P_n)(x) \,  \rd \big\{u^{(k)}(x)\big\}\\
&\quad -u^{(k)}(\theta-) (I_{1-}^{k+1} P_n)(\theta) + \int^{\theta}_{-1}(I_{1-}^{k+1} P_n)(x)\, \rd \big\{ u^{(k)}(x)\big\}\Big\}\\
&= \frac{2n+1}{2} \bigg\{ \int^{\theta}_{-1}+ \int_{[\theta,\theta]} +\int_{\theta}^{1} \bigg\} (I_{1-}^{k+1} P_n)(x)\, \rd \big\{ u^{(k)}(x)\big\}.
\end{split}
\end{equation}
By virtue of the splitting rule in Lemma \ref{LmmCarter2000Book-1},  we can summarise the formula of  the Legendre expansion coefficient with $\mu=k$ as follows.
\begin{cor}\label{sequal1}
 If $u,u',\ldots,u^{(k-1)}\in {\rm AC}([-1,1])$  and $ u^{(k)} \in  {\rm BV}([-1, 1])$ with $k \in \mathbb{N},$
 then  we have  for  all $n\ge k+1,$
\begin{equation}\label{HatUnCases1}
\begin{split}
\hat u_n^L= \frac{2n+1}{2} &\int^{1}_{-1}(I_{1-}^{k+1} P_n)(x)\, \rd \{ u^{(k)}(x)\},
\end{split}
\end{equation}
where $(I_{1-}^{k+1} P_n)(x)$ can be explicitly evaluated by \eqref{Pnmumu4}.
\end{cor}

It is seen from Theorem \ref{IdentityForUn} that the decay rate of  $\hat u_n^L$ for $u(x)$  with a fixed regularity index $\mu$ is determined by the fractional integrals of $P_n(x).$  Indeed, we have the following bound.
\begin{lmm}\label{BoundIntLeg}
For $\mu>-1/2$ and $n\ge \mu+1,$ we have
\begin{equation}\label{BoundIntLeg-0}
\max_{|x|\le 1}\Big\{\big|(I_{1-}^{\mu+1} P_n)(x)\big|, \big|(I_{-1+}^{\mu+1} P_n)(x)\big|\Big\}\le  \frac{1}{2^{\mu+1} \sqrt \pi }
	\frac{ \Gamma( (n-\mu)/ 2)}{\Gamma( (n+\mu+3)/2)}.
\end{equation}
\end{lmm}
\begin{proof} According to Szeg\"o \cite[p.\! 62]{Szego1975Book},  the generalised Jacobi polynomials  with real parameters $\alpha,\beta$ are defined by the hypergeometric functions as
\begin{equation}\label{BoundIntLeg-1}
P_n^{(\alpha,\beta)}(x) =\frac{(\alpha+1)_{n}}{n !}\, {}_2F_1\Big(\!\!-n, n+\alpha+\beta+1;\alpha+1;\frac{1-x} 2\Big),\quad x\in (-1,1),
\end{equation}
or alternatively,
\begin{equation}\label{BoundIntLeg-2}
P_{n}^{(\alpha, \beta)}(x)=(-1)^{n} \frac{(\beta+1)_{n}}{n !}{ }_{2} F_{1}\Big(\!-n, n+\alpha+\beta+1 ; \beta+1 ; \frac{1+x}{2}\Big),\quad x\in (-1,1).
\end{equation}
Recall  the Euler transform identity (cf. \cite[p.\! 95]{Andrews1999Book}): for $a,b,c \in {\mathbb R}$ and $-c\not \in {\mathbb N}_0,$
\begin{align}\label{Euler}
{}_2F_1(a,b;c;z) =(1-z)^{c-a-b}{}_2F_1(c-a,c-b;c;z),\quad |z|<1.
\end{align}
Taking $a=-n, b=n+\alpha+\beta+1, c=\alpha+1$ and $z=(1-x)/2$ in \eqref{Euler},  we obtain
 \begin{equation}\label{BoundIntLeg-3}
 \begin{split}
{}_2F_1\Big(\!\!-n, n+\alpha+\beta+1;\alpha+1;\frac{1-x} 2\Big)& =\Big(\frac {1+x} 2\Big)^{-\beta}\!
{}_2F_1\Big(\!\!\,n+\alpha+1, -n-\beta;\alpha+1;\frac{1-x} 2\Big)\\
&=\Big(\frac {1+x} 2\Big)^{-\beta}\!
{}_2F_1\Big(\!\! -n-\beta, n+\alpha+1;\alpha+1;\frac{1-x} 2\Big).

   \end{split}
\end{equation}

From \eqref{twoGGFs}, \eqref{BoundIntLeg-1}, \eqref{BoundIntLeg-3} and  ${}_2F_1(a,b;c;0)=1,$ we get
\begin{equation}\label{BoundIntLeg-4}
\begin{split}
(I_{1-}^{\mu+1} P_n)(x)&=
\frac{(1-x)^{\mu+1}}{\Gamma(\mu+2)} \frac{P_{n}^{(\mu+1, -\mu-1)}(x)}{P_{n}^{(\mu+1,-\mu-1)}(1)}=\frac{(1-x)^{\mu+1}}{\Gamma(\mu+2)}\, {}_2F_1\Big(\!\!-n, n+1;\mu+2;\frac{1-x} 2\Big)\\
&=\frac{(1-x^2)^{\mu+1}}{2^{\mu+1}\Gamma(\mu+2)}
{}_2F_1\Big(\!\! -n+\mu+1, n+\mu+2;\mu+2;\frac{1-x} 2\Big).
\end{split}

\end{equation}
Similarly, we can show that
\begin{equation}\label{BoundIntLeg-5}
\begin{split}
(I_{-1+}^{\mu+1} P_n)(x)

&=\frac{(-1)^n (1-x^2)^{\mu+1}}{2^{\mu+1}\Gamma(\mu+2)}
{}_2F_1\Big(\!\! -n+\mu+1, n+\mu+2;\mu+2;\frac{1+x} 2\Big).
\end{split}

\end{equation}
From Liu et al. \cite[Definition 2.1 \& (4.30)]{Liu19MC-Optimal}, we find that for   $\lambda \ge 1$ and $\nu\ge 0,$
	\begin{equation}\label{BoundIntLeg-6}
	\max_{|x|\le 1}\Big\{	(1-x^2)^{\lambda-\frac 12}\Big |\, {}_2F_1\Big(\!\!-\nu, \nu+2\lambda;\lambda+\frac 1 2;\frac{1\pm x} 2\Big) \Big| \Big\}\le
	 \frac{\Gamma(\lambda+ 1/2)}{\sqrt \pi}
	\frac{ \Gamma( ({\nu}+1)/ 2)}{\Gamma( ({\nu}+1)/2+\lambda)}.
	\end{equation}
Thus, taking $\nu\to n-\mu-1$ and  $\lambda\to \mu+3/2$ in \eqref{BoundIntLeg-6}, leads to
	\begin{equation}\label{BoundIntLeg-7}
	\max_{|x|\le 1}\Big\{	(1-x^2)^{\mu+1}\Big | {}_2F_1\Big(\!\!-n+\mu+1, n+\mu+2;\mu+2;\frac{1\pm x} 2\Big) \Big| \Big\}\le
	 \frac{\Gamma(\mu+2)}{\sqrt \pi}
	\frac{ \Gamma( (n-\mu)/ 2)}{\Gamma( (n+\mu+3)/2)}.
	\end{equation}

Finally the bound \eqref{BoundIntLeg-0} follows from  \eqref{BoundIntLeg-4}, \eqref{BoundIntLeg-5} and \eqref{BoundIntLeg-7}.
\end{proof}

\subsection{$L^\infty$-estimates of Legendre orthogonal projections} With the above preparations, we are now ready to analyse  the  $L^\infty$-error estimate of the $L^2$-orthogonal projection:
\begin{equation}\label{Cbexp}
u(x)=\sum_{n=0}^{\infty}\hat u_n^L\, P_n(x),\quad\;\;   (\pi_N^L  u)(x)= \sum_{n=0}^N \hat u_n^L\, P_n(x).
\end{equation}

Below, we present the approximation results on the  $L^\infty$-estimate and the weighted $L^\infty$-estimate.  We shall illustrate  that  the former is suboptimal for functions with interior singularity, but optimal for the endpoint singularity, while the latter is optimal in both cases. Such convergence behaviours  were numerically  observed in    \cite{Wang18AML-new,Wang20arXiv-How}, but lack of  theoretical justifications.
\begin{thm}\label{TruncLeg}     Let $u,u',\ldots,u^{(k-1)}\in {\rm AC}([-1,1])$ with $k \in \mathbb{N}.$
\begin{itemize}
\item[(i)] For  $\mu \in (k-1, k)$ and $\theta\in(-1,1),$ if  ${ }^{C}\! D_{\theta+}^{\mu} u \in  {\rm BV}([\theta, 1])$ and
  ${ }^{C}\! D_{\theta-}^{\mu} u \in  {\rm BV}([-1,\theta]),$ then   for  all $N\ge \mu>1/2,$
		\begin{equation}\label{FracLinfA00}
  \begin{split}
  \big\|u-\pi_N^L  u\big\|_{L^\infty (\Omega)}\le \frac{1}{2^{\mu-1} (\mu-1/2)\sqrt{\pi}}\frac{ \Gamma( ({N-\mu+1})/ 2)}{\Gamma( ({N+\mu})/2)} \, U_\theta^{(\mu)},
  \end{split}
  \end{equation}
 and for all $N\ge \mu,$
  \begin{equation}\label{FracLinfB}
  \begin{split}
  \big\|(1-x^2)^{\frac 1 4 }(u-\pi_N^L  u)\big\|_{L^\infty (\Omega)}\le \frac {1 }{2^{\mu-1} \mu \pi}\frac{ \Gamma( ({N-\mu+1})/2)}{\Gamma( ({N+\mu+1})/2)}\, \, U_\theta^{(\mu)},
  \end{split}
  \end{equation}
		where we denoted
\begin{equation}\label{seminormF0}
U_\theta^{(\mu)}:= V_{[-1,\theta]}[{ }^{C}\! D_{\theta-}^{\mu} u]  +V_{[\theta,1]}[{ }^{C}\!D_{\theta+}^{\mu} u]+ |{ }^{C}\! D_{\theta-}^{\mu} u(\theta-)| + |{ }^{C}\! D_{\theta+}^{\mu} u(\theta+)|.
\end{equation}
\item[(ii)]  If $u^{(k)}\in {\rm BV}(\bar \Omega)$ with $\Omega=(-1,1),$ then the estimates  \eqref{FracLinfA00}-\eqref{FracLinfB} with $\mu=k$ hold,  but the total variation $V_{\bar\Omega}[u^{(k)}]$ is in place of $U_\theta^{(\mu)}.$
\end{itemize}
\end{thm}
\begin{proof}

Using the identity in Theorem \ref{IdentityForUn} and the   bound in Lemma \ref{BoundIntLeg}, we obtain from \eqref{integralProp}  that
\begin{equation}\label{HatUnCaseC0}
\begin{split}
|\hat u_n^L| &= \frac{2n+1}{2}
\max_{x\in \bar\Omega} \big\{\big|(I_{1-}^{\mu+1} P_n)(x)\big|, \big|(I_{-1+}^{\mu+1} P_n)(x)\big|\big\}\,
U_\theta^{(\mu)}\\
& \le \frac{(2n+1)\Gamma( ({n-\mu})/ 2)}{ 2^{\mu+2} \sqrt{\pi}\,\Gamma( ({n+\mu}+3)/2)} \,U_\theta^{(\mu)}.
\end{split}
\end{equation}

We first prove the error bound \eqref{FracLinfA00}.  For simplicity,  we   denote
	\begin{equation}\label{Snsigma}
	{\mathcal S}_n^\mu:=\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+1})/2)},\quad
	{\mathcal T}_n^\mu:=\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu}-1)/2)}.
	\end{equation}
Using  the identity $z\Gamma(z)=\Gamma(z+1)$,  we find readily that
	\begin{equation}\begin{split}\label{FracLinfA-1}
	{\mathcal T}_n^\mu-{\mathcal T}_{n+2}^\mu&=
	\frac{{n+\mu}-1} 2\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+1})/2)}-\frac{{n-\mu}} 2\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+1})/2)}\\
	&=(\mu-1/2)\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+1})/2)}= (\mu-1/2)\,{\mathcal S}_n^\mu.
	\end{split}\end{equation}
  As $|P_n(x)|\le 1,$

	we derive from  \eqref{HatUnCaseC0} that
	\begin{equation}\begin{split}\label{FracLinfA-2}
	\big|(u -&\pi_N^L  u)(x)\big|\le \sum_{n=N+1}^{\infty}|\hat u_n^{L}| \le
\frac{U_\theta^{(\mu)}}{2^{\mu+2} \sqrt{\pi}}\sum_{n=N+1}^{\infty}
	 \frac{ (2n+1)\Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+3})/2)}\\
&\le \frac{U_\theta^{(\mu)}}{2^{\mu} \sqrt{\pi}}\sum_{n=N+1}^{\infty}
	 \frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+1})/2)}
= \frac{U_\theta^{(\mu)}}{2^{\mu} (\mu-1/2)\sqrt{\pi}}  \sum_{n=N+1}^{\infty}\big\{{\mathcal T}_n^\mu-{\mathcal T}_{n+2}^\mu\big\}

\\&=\frac{U_\theta^{(\mu)}}{2^{(\mu)} (\mu-1/2)\sqrt{\pi}} \big\{{\mathcal T}_{N+1}^\mu+{\mathcal T}_{N+2}^\mu\big\}.
	\end{split}\end{equation}
	Since $\mu>1/2,$ we obtain from \eqref{gammratio} immediately that
	\begin{equation}\label{Tsigmabnd}
	{\mathcal T}_{N+2}^\mu={\mathcal R}_{\mu-1/2}^0(1+(N-\mu)/2) \le {\mathcal R}_{\mu-1/2}^0(1+(N-\mu-1)/2)={\mathcal T}_{N+1}^\mu.
	\end{equation}
Therefore, we have from the above that
\begin{equation*}\begin{split}
	\big|(u -&\pi_N^L  u)(x)\big|\le  \frac{ 2  {\mathcal T}_{N+1}^\mu\, U_\theta^{(\mu)}}{2^{\mu} (\mu-1/2)\sqrt{\pi}}=
	\frac{U_\theta^{(\mu)}}{2^{\mu-1} (\mu-1/2)\sqrt{\pi}}\frac{ \Gamma( ({N-\mu}+1)/ 2)}{\Gamma( ({N+\mu})/2)}.
\end{split}\end{equation*}
This leads to  the error bound  \eqref{FracLinfA00}.

\medskip

We now turn to the proof of  \eqref{FracLinfB}.

Recall the Bernstein inequality (cf.\! \cite{Antonov1981VLUM-An}):
		\begin{equation}\label{IneqtyLeg-A}
	\max_{|x|\le 1}\big\{(1-x^2)^{\frac 14} |P_n(x)|\big\}\le \sqrt{ \frac{2}{ \pi}}\Big(n+\frac{1}{2}\Big)^{-\frac 12},\quad n\ge 0.

	\end{equation}
Thus we infer from \eqref{HatUnCaseC0}  and \eqref{IneqtyLeg-A} that
		\begin{equation}\begin{split}\label{FracLinfB-1}
e_N(x)& := \big|(1-x^2)^{\frac 14}  (u -\pi_N^L  u)(x)\big| \le \sum_{n=N+1}^{\infty}  \max_{|x|\le 1}\big\{(1-x^2)^{\frac 14} |P_n(x)|\big\}\, |\hat u_n^L |
	 \\
&\le
\frac{U_\theta^{(\mu)}}{2^{\mu} \pi}\sum_{n=N+1}^{\infty}
	\sqrt{\frac n2+\frac 1 4}\, \frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+3})/2)}.

	\end{split}\end{equation}
Considering   \eqref{IneqtyLeg-1} with $z=z_n=(n+\mu+1)/2$ and  $c=1/4-(\mu+1)/2\, (\le 1/4)$, we find  from  its monotonicity    that
$\widehat {\mathcal R}_c(z_n)\ge  \widehat {\mathcal R}_c(\infty)=1$   (cf. \eqref{NRatioGammaV2-1}).  This immediately implies
	\begin{equation}\begin{split}\label{FracLinfB-2}
\frac{\sqrt{n/2+1/4}}{\Gamma( ({n+\mu+3})/2)}&\le
	 \frac{1}{\Gamma( ({n+\mu})/2+1)},
	\end{split}\end{equation}
so we can bound the summation in 	\eqref{FracLinfB-1} by
		\begin{equation}\begin{split}\label{FracLinfB-3}
	\sum_{n=N+1}^{\infty}
	 \sqrt{\frac n2+\frac 14}\,\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+3})/2+1)}&\le\sum_{n=N+1}^{\infty}
	 \frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu})/2+1)}.
	\end{split}\end{equation}
	Similarly,  denoting
\begin{equation}\label{Snsigma-2}
\widehat{\mathcal S}_n^\mu:=\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu})/2+1)},\quad
\widehat{\mathcal T}_n^\mu:=\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu})/2)} ,
\end{equation}
we  find readily that
\begin{equation}\begin{split}\label{FracLinfB-4}
\widehat{\mathcal T}_n^\mu-\widehat{\mathcal T}_{n+2}^\mu& =
\frac{{n+\mu}} 2\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu})/2+1)}-\frac{{n-\mu}} 2\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu})/2+1)}

= \mu\,\widehat{\mathcal S}_n^\mu.
\end{split}\end{equation}
Following the same lines as in the derivation of \eqref{FracLinfA-2},  we can get
\begin{equation}\begin{split}\label{FracLinfB-5}
e_N(x) & \le

\frac{U_\theta^{(\mu)}}{2^{\mu}  \mu\pi} \big\{\widehat{\mathcal T}_{N+1}^\mu+\widehat{
\mathcal T}_{N+2}^\mu\big\}
\le \frac{ 2 \widehat{\mathcal T}_{N+1}^\mu}{2^{\mu}  (\mu-1)\pi}  \, U_\theta^{(\mu)},
	\end{split}\end{equation}
where we used the property derived from \eqref{gammratio} with $\mu>1,$ that is,
	\begin{equation*}\label{Tsigmabnd-2}
\widehat{\mathcal T}_{N+2}^\mu={\mathcal R}_{\mu-1}^0(1+(N-\mu)/2) \le {\mathcal R}_{\mu-1}^0(1+(N-\mu-1)/2)=\widehat{\mathcal T}_{N+1}^\mu.
	\end{equation*}
Then  the estimate \eqref{FracLinfB}  follows from \eqref{Snsigma-2} and \eqref{FracLinfB-5} straightforwardly.

\smallskip
For $\mu=k,$  using the identity in Corollary \ref{sequal1} to derive the bound in \eqref{HatUnCaseC0}, and then following the same lines, we can obtain the estimates in (ii).
\end{proof}

 Under the regularity assumption in Theorem \ref{TruncLeg},  we infer from \eqref{NRatioGammaV2-1} and the estimates \eqref{FracLinfA00}-\eqref{FracLinfB} the convergence behaviour:
\begin{equation}\label{converorg}
\|u-\pi_N^L  u\|_{L^\infty (\Omega)}=O(N^{-\mu+1/2}),\quad  \|(1-x^2)^{\frac 1 4 }(u-\pi_N^L  u)\|_{L^\infty (\Omega)}=O(N^{-\mu}),
\end{equation}
which exhibit  a half-order convergence difference. Moreover,  the estimate \eqref{FracLinfB-1} implies

\begin{equation}\label{converorg2}
 \big| (u -\pi_N^L  u)(x)\big|= (1-x^2)^{-\frac 14}\, O(N^{-\mu}),\quad\forall\, x\in [-a,a]\subset (-1,1).
\end{equation}
For a  function with an interior singularity  in $[-a,a]$ with $|a|<1$,   one expects the optimal order  $O(N^{-\mu}).$
Note  from  \eqref{FracLinfA-2} and \eqref{FracLinfB-1} that the bounds essentially depend on the maximum of $|P_n(x)|$
and $(1-x^2)^{1/4}|P_n(x)|,$ which behave very differently near the endpoints  as
shown in Figure \ref{PnWithWeight}.  In fact, $|P_n(x)|=O(n^{-1/2})$ for $x\in [-a,a],$ but it is overestimated by the bound $1$ at $x=\pm 1.$ However, from \eqref{IneqtyLeg-A}, we have $(1-x^2)^{1/4}|P_n(x)|\le Cn^{-1/2}$ for all $x\in [-1,1].$ This is actually the cause of  the lost order  in  the (non-weighted) $L^\infty$-estimate in \eqref{converorg}.
\begin{figure}[!ht]
	\begin{center}
		{~}\hspace*{-20pt}
\includegraphics[width=0.42\textwidth]{Pn-Left2.pdf}\qquad\quad
		\includegraphics[width=0.42\textwidth]{Pn-Right2.pdf}
		\caption{$P_n(x)$ and $(1-x^2)^{\frac 14}P_n(x)$ with $x\in [-1,1]$ and $n=100.$}
		\label{PnWithWeight}
	\end{center}
\end{figure}

With these analysis tools  at hand,  we  further examine $u(x)=|x|$ (as a motivative example in Trefethen  \cite{Trefethen2008SIREV}), for which
 Wang \cite{Wang18AML-new,Wang20arXiv-How} observed the order $O(N^{-1})$ numerically,  but  the order is $O(N^{-1/2})$ based on the error estimate of Legendre approximation in $L^\infty$-norm.    From the pointwise error plots in Figure \ref{PoWisePosi} (left),  we see the largest error occurs at the singular point $x=0.$ Indeed, we have the following estimates (with the proof given in Appendix  \ref{AppendixC}), which are sharp as shown in Figure \ref{PoWisePosi} (right).
  \begin{theorem}\label{01estimate} Consider $u(x)=|x|$ for  $x\in [-1,1].$  Then for $N>2,$ we have
 \begin{equation}\label{un|x|-1+1}
\begin{split}
\big|(u -\pi_N^L  u)(0)\big|\le \frac{2}{\pi(N-1)}\,; \; \quad \big|(u-\pi_N^L  u)(\pm 1)\big|&\le
\frac{1 }{2\sqrt{\pi}}\frac{\Gamma( N/2-1)} {\Gamma( N/2+1/2)}.
\end{split}
\end{equation}

 \end{theorem}
 \begin{figure}[!ht]
	\begin{center}
		{~}\hspace*{-20pt}
\includegraphics[width=0.42\textwidth]{PoWisePosi.pdf}\qquad\quad
		\includegraphics[width=0.42\textwidth]{ErrorAbsXZerosPm1.pdf}
		\caption{Left: $|(u -\pi_N^L  u)(x)|$  with $u(x)=|x|$ and different $N.$ Right: Tightness of the bounds in Theorem \ref{01estimate}. }
		\label{PoWisePosi}
	\end{center}
\end{figure}

Finally, we apply the main results to the example $u(x)=|x|^\mu$ with $\mu \in(k-1,k).$  As shown in \cite[Thm. 4.3]{Liu19MC-Optimal},
 $u,u',\ldots,u^{(k-1)}\in {\rm AC}(\bar\Omega),$ ${ }^{C}\! D_{0+}^{\mu} u \in  {\rm BV}([0, 1]),$
 and ${ }^{C}\! D_{0-}^{\mu} u \in  {\rm BV}([-1,0]).$ Thus, we infer from \eqref{converorg} that the expected convergence orders are
 $O(N^{-\mu+1/2})$ in $L^\infty$-norm and  $O(N^{-\mu})$ in $L^\infty_{\varpi}$-norm with $\varpi=(1-x^2)^{1/4}.$ Observe from
 the numerical results in Table \ref{Tab-1}  that the latter is optimal, but the former loses half order.

\begin{table}[!htbp]
 \centering
 \caption{Convergence order of Legendre expansion for $|x|^\mu$.}
\small
 \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
  \hline
  \multicolumn{1}{ |c  }{\multirow{2}{*}{$N$} } &
  \multicolumn{4}{ |c| }{Errors in $L^\infty$-norm} & \multicolumn{4}{ |c| }{Errors in $L^\infty_\varpi$-norm}      \\ \cline{2-9}
  \multicolumn{1}{ |c  }{}                        &
  \multicolumn{1}{ |c| }{$\mu=1.7$}  & order  & $\mu=2.6$  & order  &$ \mu=1.7$  & order & $\mu=2.6$  & order \\ \hline
$ 2^3 $ & 5.81e-03 &--& 2.35e-03 &-- &5.81e-03 &--& 2.22e-03 &--  \\
$ 2^4 $ & 2.03e-03 & 1.52 & 4.38e-04 & 2.42 & 2.03e-03 & 1.52 & 4.38e-04 & 2.34 \\
$ 2^5 $ & 6.72e-04 & 1.60 & 8.01e-05 & 2.45  & 6.72e-04 & 1.60 & 8.01e-05 & 2.45 \\
$ 2^6 $ & 2.15e-04 & 1.65 & 1.40e-05 & 2.52 & 2.15e-04 & 1.65 & 1.40e-05 & 2.52  \\
$ 2^7 $ & 6.74e-05 & 1.67 & 2.37e-06 & 2.56  & 6.74e-05 & 1.67 & 2.37e-06 & 2.56\\
$ 2^8 $ & 2.09e-05 & 1.69 & 3.97e-07 & 2.58 & 2.09e-05 & 1.69 & 3.97e-07 & 2.58\\
  \hline
 \end{tabular} \label{Tab-1}
\end{table}

\subsection{$L^2$-estimates}
As pointed out in \cite[Chap.\!\! 3]{Shen2011Book}, the estimate of the $L^2$-orthogonal projection is the starting point  to derive many other approximation results that  provide fundamental tools for error analysis of spectral and $hp$ methods  (see, e.g., \cite{Bernardi1997Book,Canuto2006Book,Guo1998Book,Hesthaven2007Book,Schwab1998Book,Shen2011Book}).  Most    estimates  therein are for functions in Sobolev or Besov spaces. Here, we consider functions with
 AC-BV regularity, thereby enriching the approximation theory.

We first highlight the fundamental importance of estimating $L^2$-orthogonal projection in \eqref{Cbexp}.
For $u\in H^1(\Omega),$ we define
\begin{equation}\label{H1proj}
(\pi_{N}^1 u)(x)= u(-1)+\int^x_{-1}\pi_{N-1}^L u'(t)\, \rd t\in {\mathcal P}_{\!N},
\end{equation}
where ${\mathcal P}_{\!N}$ denotes the set of  polynomials of degree at most $N.$ Note that for $N\ge 2,$ we have from the orthogonality of Legendre polynomials that
$$(\pi_{N}^1 u)(1)= u(-1)+\int^1_{-1}\pi_{N-1}^L u'(t)\,\rd t=  u(-1)+ \int^1_{-1} u'(t)\,\rd t=u(1).$$
Thus we have $(\pi_{N}^1 u)(\pm 1)=u(\pm 1).$ Moreover, one verifies readily that
$$\int^1_{-1}\big( \pi_{N}^1 u-u\big)'(x) \,v'(x)\,\rd x =0, \;\;\; \forall v\in \mathcal{P}_{\!N}^0:=\{v\in {\mathcal P}_{\!N} : v(\pm 1)=0 \}.$$
Therefore, \eqref{H1proj} defines the $H^{1}_0$-orthogonal projection. Note that
\begin{equation}\label{H1proj2}
\|(u-\pi_{N}^1 u)'\|_{L^{2}(\Omega)}^{2}=\|u'-\pi_{N-1}^Lu' \|_{L^{2}({\Omega})}^{2},
\end{equation}
so the $H^1$-estimate boils down to the estimate of the $L^2$-orthogonal projection. The high-order $H_0^m$-orthogonal projection is treated similarly in a recursive manner (see, e.g.,  \cite{Bernardi1997Book}). On the other hand, the analysis of Gauss-type interpolation and quadrature errors is also based upon the Legendre expansion (see \cite{Shen2011Book}).

With tools in Subsection \ref{sub31}, we can also derive  the following optimal  $L^2$-error bound under  the  AC-BV regularity of $u$, from which we can further establish many other approximation results indispensable for analysis of spectral and $hp$ methods for PDEs. Here, we omit such extensions.
\begin{thm}\label{newL2est} Assume the conditions in  Theorem {\rm \ref{TruncLeg}} hold.
\begin{itemize}
\item[(i)]
 For $\mu\in(k-1,k)$ and $ -1/2<\mu< N,$
		\begin{equation}\label{FracL2B}
		\begin{split}
		\|u-\pi_N^L  u\|_{L^2(\Omega)}\le \sqrt{\frac{2}{(2\mu+1) \pi}
		\frac{ \Gamma( {N-\mu})}{\Gamma( {N+\mu+1})}} \, U_\theta^{(\mu)}.
		\end{split}
		\end{equation}

\item[(ii)] For $\mu=k,$
 the estimates  \eqref{FracL2B} hold,  with the total variation $V_{\bar\Omega}[u^{(k)}]$ is in place of $U_\theta^{(\mu)}.$
\end{itemize}
\end{thm}
\begin{proof} Similar to   \eqref{Tsigmabnd}, we can use \eqref{gammratio} to show that
$$\frac{ \Gamma( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+3})/2)}\le
	\frac{ \Gamma( ({n-\mu-1})/ 2)}{\Gamma( ({n+\mu})/2+1)}.$$
Then from \eqref{DuplicationFormula},	we derive
	\begin{equation}\begin{split}\label{FracL2-4}
	& \frac{ (n/2+1/4)\Gamma^2( ({n-\mu})/ 2)}{\Gamma^2( ({n+\mu+3})/2)} \le \frac{ \Gamma^2( ({n-\mu})/ 2)}{\Gamma( ({n+\mu+3})/2)\Gamma( ({n+\mu+1})/2)}\\
	&\qquad\qquad  \le \frac{ \Gamma( ({n-\mu})/ 2)\Gamma( ({n-\mu-1})/ 2)}{\Gamma( ({n+\mu})/2+1)\Gamma( ({n+\mu+1})/2)}=2^{2(\mu+1)}\frac{ \Gamma( {n-\mu-1})}{\Gamma( {n+\mu+1})}\\
&\qquad\qquad  =\frac{2^{2(\mu+1)}}{2\mu+1}
\bigg(\frac{ \Gamma( {n-\mu-1})}{\Gamma( {n+\mu})}-\frac{ \Gamma( {n-\mu})}{\Gamma( {n+\mu+1})}\bigg).
	\end{split}\end{equation}

Then, by the orthogonality of Legendre polynomials, we derive from
\eqref{HatUnCaseC0} and  \eqref{FracL2-4} that for $\mu> -1/2$,
	\begin{equation}\label{FracL2-5}
	\begin{split}
	\big\|u-\pi_N^{L} u\big\|_{L^2(\Omega)}^2=&\sum_{n=N+1}^{\infty}\frac{2}{2n+1}\big|\hat u_n^{L}\big|^2
	\le   \frac{(U_\theta^{(\mu)})^2}{2^{2\mu+3} \pi}\sum_{n=N+1}^{\infty}
	 \frac{ (2n+1)\Gamma^2( ({n-\mu})/ 2)}{\Gamma^2( ({n+\mu+3})/2)}\\
\le &   \frac{2(U_\theta^{(\mu)})^2}{(2\mu+1) \pi}\frac{ \Gamma( {N-\mu})}{\Gamma( {N+\mu+1})}.
	\end{split}
	\end{equation}

For $\mu=k,$  using the identity in Corollary \ref{sequal1} to derive the bound in \eqref{HatUnCaseC0}, and then following the same lines, we can obtain the estimates in (ii).
		\end{proof}

\begin{comment}
\subsection{Some point-wise error estimates}\label{Subsect:Pt-est}     In particular, the very recent work by  Babu\v{s}ka and  Hakula \cite{Babuska19CMAME-Pointwise} provided  some delicate point-wise estimates for the specific function (as in the seminal work \cite{Gui1986NM-Part1}): $u(x)=(x-\theta)^\beta$ for $\theta<x\le 1$ and $u(x)=0$ for $-1\le x<\theta$  with $|\theta|<1$ and $\beta>0.$

Using Corollary \ref{CorLinfL2} with $k=1$ and $V_{\bar \Omega}[u^{(k)}]=2$,  we obtain from \eqref{FracLinfAV2-1} and \eqref{FracLinfBV2-1} the following estimates.
\begin{cor} Consider $u(x)=|x|$ on $\bar \Omega.$ Then for  $N\ge 1,$ we have
\begin{equation}\label{FracLinfA-V2}
		\begin{split}
		\big\|u-\pi_N^L  u\big\|_{L^\infty (\Omega)}\le \frac{4}{\sqrt{\pi}}\frac{ \Gamma( N/ 2)}{\Gamma(N/2+ 1/2)},
		\end{split}
		\end{equation}
 and
		\begin{equation}\label{FracLinfB-V2}
		\begin{split}
		\big\|(1-x^2)^{\frac 1 4 }(u-\pi_N^L  u)\big\|_{L^\infty (\Omega)}\le \frac{4}{ \pi N}.
		\end{split}
		\end{equation}
\end{cor}
In fact, the estimate \eqref{FracLinfA-V2} is of order $O(N^{-1/2}),$ but the optimal order should be  $O(N^{-1}).$
 How to derive the optimal result appears  open  (cf. \cite{Wang18AML-new}).  To this end, we unfold the mystery behind this.

Some important observations are in order.
\begin{itemize}
\item It is seen from \eqref{P2kvalue} that $|P_{2m}(0)|=O(m^{-1/2}).$ However, it is overestimated in the proof of Theorem \ref{TruncLeg},  which is bounded as $|P_{2m}(0)|\le 1$ (see \eqref{FracLinfA-2} and this type of estimates in e.g., \cite{Wang18AML-new,Xiang20NM-Optimal,Wang20arXiv-How}). This should be blamed for the loss of  order at this point.
Observe from Figure \ref{PoWisePosi} (right) that the estimates are sharp.

\item  As shown in Figure \ref{PnWithWeight} (left),  the Legendre polynomial  $P_n(x)$ attends its maximum/minimum at the endpoints $x=\pm 1,$ but $|P_n(x)|=O(n^{-1/2})$ for $x\in [-a,a]$ where $a=1-\delta$ for some small $\delta.$
In fact, we suppose to have the optimal point-wise error estimate as \eqref{un|x|-1+1} for all $x\in [-a,a].$

\item It is noteworthy that the (weighted) $L^\infty$-estimates including \eqref{FracLinfB} in Theorem  \ref{TruncLeg} 	 and
\eqref{FracLinfBV2}  in Corollary \ref{CorLinfL2}. As shown in Figure \ref{PnWithWeight} (right), the extreme values of   $(1-x^2)^{\frac 1 4} P_n(x)$ behave like $O(n^{-1/2})$ as the global maximum as in \eqref{FracLinfB-1}. The situation is very similar to the Chebyshev polynomial.   In fact, we can also see  from   the property (cf. \cite[Thm. 2.1]{Liu20JAT-Asymptotics}): for  $ n>1/2$  and  $\varphi\in (0,\pi),$
 \begin{equation}\label{IntRep-0}
 \begin{split}
 (\sin \varphi)^{\frac 12}\,P_n(\cos \varphi)=\frac 2 {\sqrt{\pi(2n+1)}}

  {\cos ((n+1/2)\varphi-\pi/4)}
 +{\mathcal  R}_n (\varphi),

 \end{split}\end{equation}
where
\begin{equation}\label{uniformBndR}
 |{\mathcal  R}_n(\varphi)|\le \frac{\sqrt{2}}{4\sqrt{\pi}(n-1/2)^{3/2}} \Big\{|\cot \varphi|+\frac{1} {n-1/2}\Big\}, \quad \forall\, \varphi\in (0,\pi).
\end{equation}
\item In view of the above,   we can obtain from \eqref{FracLinfB} in Theorem  \ref{TruncLeg}  the following point-wise estimate:

\begin{equation}\label{uptestim}
\begin{split}
|(u-\pi_N^L  u)(x)| &\le (1-x^2)^{-\frac 1 4 } \|(1-x^2)^{\frac 1 4 }(u-\pi_N^L  u\|_{L^\infty (\Omega)}\\
& \le  (1-x^2)^{-\frac 1 4 }  \frac {U_\theta^{m,s} }{2^{\sigma-2} (\sigma-1)\pi}\frac{ \Gamma( ({N-\sigma})/2+1)}{\Gamma( ({N+\sigma})/2)},
\end{split}
\end{equation}
for $x\in (-1,1)$. As a result, we can derive optimal estimate for functions with interior singularities from \eqref{uptestim}.

\end{itemize}

\begin{figure}[!ht]
	\begin{center}
		{~}\hspace*{-20pt}
\includegraphics[width=0.38\textwidth]{Pn-Left2}\qquad\quad
		\includegraphics[width=0.38\textwidth]{Pn-Right2}
		\caption{$P_n(x)$ and $(1-x^2)^{\frac 14}P_n(x)$ with $x\in [-1,1]$ and $n=100.$}
		\label{PnWithWeight}
	\end{center}
\end{figure}
\end{comment}

\section{Legendre expansion of functions with endpoint singularities}\label{sect4main}
\setcounter{equation}{0}
\setcounter{lmm}{0}
\setcounter{thm}{0}
\setcounter{cor}{0}

The aforementioned AC-BV framework and  main results can be extended to the study of  the end-point singularities, which typically occur in underlying solutions  of  PDEs  in various situations, for instance, irregular domains, singular coefficients and mismatch of boundary conditions among others. It is known that  the  Legendre expansion of a function with an endpoint singularity  has  a much higher convergence rate than that with an interior singularity of the same type.  We illustrate this through an example which also motivates the seemingly complicated extension. To fix the idea, we focus on the left endpoint singularity but the results can be extended to the right endpoint setting straightforwardly.
\subsection{An illustrative example} We consider $u(x)=(1+x)^\mu g(x)$ with $\mu\in (k-1,k), k\in\mathbb N$ and
 and a sufficiently smooth $g(x)$ on $\Omega.$ Then we can  write
\begin{equation}\label{newA}
u(x)=(1+x)^\mu g(x)=\sum_{m=0}^\infty\frac{g^{(m)}(-1)}{m!} (1+x)^{\mu+m}.

\end{equation}
Then by  \eqref{intformu},
\begin{equation}\label{Dcregul}
\begin{split}
({ }^{C}\! D_{-1+}^{\mu} u)(x)=\big(I_{-1+}^{k-\mu} u^{(k)}\big)(x) & = \sum_{m=0}^\infty\frac{\{\mu+m\}_{k}}{m!} g^{(m)}(-1) \,
I_{-1+}^{k-\mu}\big\{ (1+x)^{m+\mu-k}\big\}\\
& = \sum_{m=0}^\infty\frac{\{\mu+m\}_{k}\, \Gamma(m+\mu-k+1)}{(m!)^2} g^{(m)}(-1) \,
 (1+x)^{m},
\end{split}
\end{equation}
where  $\{a\}_k=a(a-1)\cdots (a-k+1)$ stands for the falling factorial. This implies  $({ }^{C}\! D_{-1+}^{\mu} u)(x)$ is sufficiently smooth. In particular, if $g=1,$ then $({ }^{C}\! D_{-1+}^{\mu} u)(x)$ is equal to a constant.

We deduce from  Theorem  \ref{IdentityForUn} with $\theta\to-1^+$ that  for   $u, u',\cdots, u^{(k-1)}\in {\rm AC}(\bar \Omega),$ and
$  { }^{C}\! D_{-1+}^{\mu}u  \in  {\rm BV}(\bar\Omega)$ with $\mu \in (k-1, k],$ we have
\begin{equation}\label{HatUnCaseC-3Old-1}
\begin{split}
\hat u_n^L = \frac{2n+1}{2} \Big\{
  ({{ }^{C}\! D_{-1+}^{\mu} u)(-1+)} (I_{1-}^{\mu+1} P_n)(-1)
+\int_{-1}^{1} (I_{1-}^{\mu+1} P_n)(x) \,  \rd \big\{{ }^{C}\! D_{-1+}^{\mu} u(x)\big\}\Big\}.

\end{split}
\end{equation}

In fact, we can  show   that $(I_{1-}^{\mu+1} P_n)(-1)\sim n^{-2(\mu+1)},$ so the first term decays like $O(n^{-2\mu-1})$ (which gives  the optimal convergence order for \eqref{newA} (see Table \ref{UnV2}),  and  doubles $O(n^{-\mu-1/2})$ for the interior singularity, e.g., of $|x|^\mu g(x)$).

\begin{lemma}\label{Pnfyi} For $n\ge \mu+1>0,$ we have
\begin{equation}\label{HatUnCaseC-3Old-1+1}
\begin{split}
(I_{1-}^{\mu+1} P_n)(-1)&=  \frac{(-1)^n2^{\mu+1}\Gamma(\mu+1)\sin((\mu+1)\pi)}{\pi} \frac{\Gamma(n-\mu)}{\Gamma(n+\mu+2)}.
\end{split}
\end{equation}
\end{lemma}
\begin{proof}  By \eqref{twoGGFs}, we have
\begin{equation}\label{twoGGFs0}
\begin{split}
&(I_{1-}^{\mu+1} P_n)(-1)=
\frac{2^{\mu+1}}{\Gamma(\mu+2)} \frac{P_{n}^{(\mu+1, -\mu-1)}(-1)}{P_{n}^{(\mu+1,-\mu-1)}(1)}
=
\frac{(-1)^n 2^{\mu+1}}{\Gamma(\mu+2)} \frac{P_{n}^{(-\mu-1, \mu+1)}(1)}{P_{n}^{(\mu+1,-\mu-1)}(1)},
\end{split}
\end{equation}
where we used the parity $P_n^{(\alpha,\beta)}(-x)=(-1)^n P_{n}^{(\beta,\alpha)}(x)$ valid for all real parameters $\alpha,\beta$
(cf.\!   Szeg\"o  \cite[p.\! 64]{Szego1975Book}).  Then we derive \eqref{HatUnCaseC-3Old-1+1} from \eqref{Pnmumu} and \eqref{nonitA} immediately.
\end{proof}

We find from Lemma \ref{BoundIntLeg} that the second term in \eqref{HatUnCaseC-3Old-1} decays at a rate
$O(n^{-\mu-1/2}),$ if one naively works this out with this formula.  However, in view of \eqref{Dcregul},  we can continue to carry out integration by parts upon \eqref{HatUnCaseC-3Old-1} as many as times we want, until the first  boundary term in \eqref{HatUnCaseC-3Old-1} dominates the error.
 This produces the optimal order $O(n^{-2\mu-1})$ (see Table \ref{UnV2} for numerical illustrations).
\begin{table}[!htbp]
	\centering
	\caption{Decay rate of $|\hat u_n^{L}|$ with  $u= (x+1)^\mu \sin x$.}
\small
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
	$n$&{$\mu=0.1$}  & order
		&$ \mu=1.2$  & order  & $\mu=2.6$  & order
		\\
\hline
$2^3$&1.26e-02&--& 6.86e-04&--& 1.42e-04&--\\
$2^4$&5.59e-03&1.17& 6.59e-05&3.38& 1.35e-06&6.72\\
$2^5$&2.47e-03&1.18& 6.41e-06&3.36& 1.86e-08&6.18\\
$2^6$&1.08e-03&1.19& 6.20e-07&3.37& 2.60e-10&6.16\\
$2^7$&4.73e-04&1.19& 5.94e-08&3.38& 3.49e-12&6.22\\
\hline
	\end{tabular} \label{UnV2}
\end{table}

 \begin{comment}
 Indeed, we obtain from \cite[(2.29)]{Liu19MC-Optimal} that the GGF satisfies
\begin{equation}\label{endSingularities-2}
\lim_{x \to -1^+} \big\{  (1-x^2)^\sigma \, {}^{r\!}G_{n-\sigma}^{(\sigma+1/2)}(x)\big\}=(-1)^n
2^{2\sigma}
\frac{\sin(\sigma\pi)}\pi
\frac{\Gamma(\sigma)\Gamma(\sigma+1)\Gamma(n-\sigma+1)}{\Gamma(n+\sigma+1)},
\end{equation}
where $\sigma=m+s=\alpha+1.$ Moreover,  by \eqref{newA} with $m=0,$ we have
\begin{equation}\label{newcaselim}
\lim_{x\to -1^+} I_{-1+}^{1-s} \big\{u^{(m)}(x)\big\}=\Gamma(\alpha-m+1) g(-1).
\end{equation}
Recall that {\rm(}cf.  \cite[(5.11.13)]{Olver2010Book}{\rm):} for $a<b$,
\begin{equation}\label{gamratioA}
\frac{\Gamma(z+a)}{\Gamma(z+b)}=z^{a-b}+\frac{1}{2}(a-b)(a+b-1)z^{a-b-1}+O(z^{a-b-2}),\quad z\gg 1.
\end{equation}
Thus, the second boundary term in  \eqref{HatUnCaseC-3Old-1} for the function \eqref{gxcase} decays like $O(n^{-2\sigma})=O(N^{-(1+2\alpha)}).$ We refer to Table \ref{UnV2} for the numerical illustration.

\begin{table}[!htbp]
	\centering
	\caption{Decay rate of $|\hat u_n^{L}|$ with  $u= (x+1)^\alpha \sin x$.}
\small
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
	$n$&{$\alpha=0.1$}  & order
		&$ \alpha=1.2$  & order  & $\alpha=2.6$  & order
		\\
\hline
$2^3$&1.26e-02&--& 6.86e-04&--& 1.42e-04&--\\
$2^4$&5.59e-03&1.17& 6.59e-05&3.38& 1.35e-06&6.72\\
$2^5$&2.47e-03&1.18& 6.41e-06&3.36& 1.86e-08&6.18\\
$2^6$&1.08e-03&1.19& 6.20e-07&3.37& 2.60e-10&6.16\\
$2^7$&4.73e-04&1.19& 5.94e-08&3.38& 3.49e-12&6.22\\
\hline
	\end{tabular} \label{UnV2}
\end{table}

\begin{thm}[{\bf Fractional Taylor formula}]\label{NewFTFV1-V2}  Let  $\mu \in (k-1, k]$ with $k \in \mathbb{N},$ and  let $f(x)$ be a real function that is  $(k-1)$times differentiable at the point $x=\theta$.
\begin{itemize}
 \item[(i)] If $f^{(k-1)}\in {\rm AC}([\theta, x]),$  let  $v={ }^{C}\! D_{\theta+}^{\mu} f,$ and $v, v',\cdots,v^{(m-1)}\in  {\rm AC}([\theta, x]),$  $v^{(m)}\in {\rm BV}([\theta, x])$  then we have the left-sided fractional Taylor formula
\begin{equation}\begin{split}\label{LetfTaylor-00-V2}
f(x)=&\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j} +(x-\theta)^{\mu} \sum_{i=0}^{m} \frac{v^{(i)}(\theta+)}{\Gamma(\mu+i+1)}{(x-\theta)^{i}}\\
&+  \frac{1}{\Gamma(\mu+m+1)}
 \int_{\theta}^{x} (x-t)^{\mu+m}\, \rd\, \{v^{(m)}(t)\}.
\end{split}\end{equation}
\item[(ii)] If $f^{(k-1)}\in {\rm AC}([x,\theta]),$  let  $v={ }^{C}\! D_{\theta-}^{\mu} f,$ and   $v, v',\cdots,v^{(m-1)}\in  {\rm AC}([x,\theta]),$  $v^{(m)}\in {\rm BV}([x,\theta])$  then we have the right-sided fractional Taylor formula
    \begin{equation}\begin{split}\label{LetfTaylor-000-V2}
f(x)=&\sum_{j=0}^{k-1} \frac{f^{(j)}(\theta)}{j!}(x-\theta)^{j} +(\theta-x)^{\mu} \sum_{i=0}^{m} \frac{v^{(i)}(\theta-)}{\Gamma(\mu+i+1)}{(x-\theta)^{i}}\\
&-  \frac{(-1)^m}{\Gamma(\mu+m+1)}
 \int^{\theta}_{x} (t-x)^{\mu+m}\, \rd\, \{v^{(m)}(t)\}.
\end{split}\end{equation}

\end{itemize}
\end{thm}
\begin{proof}
Using integration by parts for \eqref{LetfTaylor-00} and \eqref{LetfTaylor-000}, we obtain \eqref{LetfTaylor-00-V2} and \eqref{LetfTaylor-000-V2}, respectively.
\end{proof}
\end{comment}

\subsection{Approximation results for functions with endpoint singularities}
With the above understanding, we are now ready to present the identity on
 the Legendre expansion coefficient from \eqref{HatUnCaseC-3Old-1} and integration by parts. Given that  the function has more regularity in this case,  we make the following assumption.
\begin{defn}[{\bf Regularity Assumption}]\label{functionA-1} For $\mu \in (k-1, k]$ with $k\in \mathbb N,$  assume
 $u, \cdots, u^{(k-1)}\in {\rm AC}(\bar \Omega)$ and
$ v_\mu(x):={ }^{C}\! D_{-1+}^{\mu}u  \in  {\rm BV}(\bar\Omega).$  We further assume that $v_\mu, \cdots,v^{(m-1)}_\mu\in  {\rm AC}(\bar \Omega)$ and  $v^{(m)}_\mu\in {\rm BV}(\bar \Omega).$   Accordingly, we denote
 \begin{equation}\label{seminormF0-1}
  \begin{split}
  U^{(\mu,m)}_-:= V_{\bar\Omega}[v^{(m)}_\mu]  + \,|\sin(\mu\pi)| \sum_{l=0}^{m}\big| v^{(l)}_\mu (-1+) \big|.

  \end{split}
  \end{equation}
For simplicity, {\em we say $u$ is  of  {\rm AC-BV$_{\!\!\mu,m}$-regularity}.} \qed
\end{defn}

 Note that for the example \eqref{newA}, $v_\mu$ is sufficiently smooth so we have $m=\infty.$ Under this assumption, we can update
the formula  \eqref{HatUnCaseC-3Old-1} as follows.
\begin{thm}\label{Thm:52} Assume that $u$ is  of  {\rm AC-BV$_{\!\!\mu,m}$-regularity}.
Then  for $n\ge \mu+m+1,$ we have
		\begin{equation}\label{Unend-0}
		\begin{split}
		\hat u_n^{L}& = \frac{2n+1}{2} \bigg\{  \sum_{l=1}^{m}(I_{1-}^{\mu+l+1} P_n)(-1)v^{(l)}_\mu(-1+)+\int_{-1}^{1}(I_{1-}^{\mu+m+1} P_n)(x) \, {\rm d}\{ v^{(m)}_\mu(x)\}\bigg\},
		\end{split}
		\end{equation}
where $(I_{1-}^{\mu+l+1} P_n)(-1)$ has the explicit value given by \eqref{HatUnCaseC-3Old-1+1}.

	\end{thm}
	\begin{proof} Since  $v_\mu, \cdots,v^{(m-1)}_\mu\in {\rm AC}(\bar \Omega),$ we can conduct  integration by parts upon \eqref{HatUnCaseC-3Old-1}:
	\begin{equation*}\label{HatUnCaseC-3Old-10}
\begin{split}
\frac{2\,  \hat u_n^L} {2n+1} & =
 (I_{1-}^{\mu+1} P_n)(-1)\, v_\mu(-1)
+\int_{-1}^{1} (I_{1-}^{\mu+1} P_n)(x) \, v_\mu'(x)\,  \rd x\\
& =\cdots=\sum_{l=1}^{m-1}(I_{1-}^{\mu+l+1} P_n)(-1)v^{(l)}_\mu(-1)+\int_{-1}^{1}(I_{1-}^{\mu+m} P_n)(x) v^{(m)}_\mu(x)  \,  \rd x\\
& = \sum_{l=1}^{m}(I_{1-}^{\mu+l+1} P_n)(-1)v^{(l)}_\mu(-1+)+\int_{-1}^{1}(I_{1-}^{\mu+m+1} P_n)(x) \, {\rm d}\{ v^{(m)}_\mu(x)\},
\end{split}
\end{equation*}
where the boundary values at $x=1$ vanish in view of \eqref{twoGGFs}, and in the last step, we used the factor
$v^{(m)}_\mu\in {\rm BV}(\bar \Omega)$ and \eqref{IPPW1122}.
			\end{proof}

Comparing the formulas of $\hat u_n^L$ in Theorem \ref{IdentityForUn} (with $\theta\to -1^+$, i.e., \eqref{HatUnCaseC-3Old-1})  and Theorem \ref{Thm:52},
we find they  largely differ from the regularity index.  We can use Lemmas \ref{BoundIntLeg} and \ref{Pnfyi} to deal with the fractional integrals of the Legendre polynomial. Accordling,  we can follow the same lines as in the proofs of Theorems \ref{TruncLeg} and \ref{newL2est} to derive the following estimates. To avoid the repetition, we skill the proof, though there is  subtlety in some derivations.
\begin{thm}\label{TruncLegEnd} Assume that $u$ is  of   {\rm AC-BV$_{\!\!\mu,m}$-regularity}. Then  we have the following estimates.

	\begin{itemize}
\item[(i)] For $\mu>1/2$ and $N\ge \mu+m,$
		\begin{equation}\label{FracEndLinfA}
		\begin{split}
		\|u-\pi_N^L  u \|_{L^\infty (\Omega)}&\le \bigg\{\frac{1}{2^{\mu+m-1} (\mu+m-1/2)\sqrt{\pi}}\frac{ \Gamma( ({N-\mu-m+1})/ 2)}{\Gamma( ({N+\mu+m})/2)}\\
&\qquad  +\sum_{j=0}^{m}\frac{2^{\mu+j}\Gamma(\mu+j+1)}{\pi(\mu+j-1)}\frac{\Gamma(N-\mu-j+1)}{\Gamma(N+\mu+j+1)}
\bigg\}\, U^{(\mu,m)}_-.
		\end{split}
		\end{equation}

		\item[(ii)] For $ \mu>-1/2$ and $N> \mu+m,$
		\begin{equation}\label{FracEndL2B}
		\begin{split}
		\|u-\pi_N^L  u\|_{L^2(\Omega)}&\le\bigg\{ \frac{4}{(2\mu+2m+1) \pi}
		\frac{ \Gamma( {N-\mu-m})}{\Gamma( {N+\mu+m+1})}\\
&\qquad + \frac{2^{6\mu+8}\Gamma^2(\mu+1)}{\pi^2(4\mu+2)} \frac{(N+1)^2\,\Gamma(2N-2\mu+1)}{(2N+1)^2\,\Gamma(2N+2\mu+3)}
\bigg\}^{\frac 12}\, U^{(\mu,m)}_-.
		\end{split}
		\end{equation}
	\end{itemize}

\end{thm}
\begin{comment}
	\begin{proof}
(i) From Lemma \ref{BoundIntLeg} and \eqref{seminormF0-1}-\eqref{Unend-0}, we obtain
for $\mu\ge -1/2$ and $n\ge \mu+m+1,$
	\begin{equation}\label{BoundForUnEnd}
	\begin{split}
	|\hat u_n^{L}| &\le  \bigg\{\frac{1}{2^{\mu+m+2} \sqrt{\pi}}
	\frac{ (2n+1)\Gamma( ({n-\mu-m})/ 2)}{\Gamma( ({n+\mu+m+3})/2)}
		+\widehat  C_{n,\mu}\sum_{j=0}^{k}\frac{\widehat  C_{n,\mu+j}}{\widehat  C_{n,\mu}}\bigg\}\, U^{(\mu,m)}_-,
	\end{split}
	\end{equation}
where

\begin{equation}\label{newconstants-2}
\begin{split}
\widehat C_{n,\mu}:&= 2^{\mu}\Gamma(\mu+1)\pi ^{-1} \frac{(2n+1)\Gamma(n-\mu)}{\Gamma(n+\mu+2)}.
\end{split}\end{equation}
 We now prove \eqref{FracEndLinfA}.
By  \eqref{BoundForUnEnd} with $|P_n(x)|\le 1,$ we obtain
	\begin{equation}\begin{split}\label{FracEndLinfA-2}
	\big|(u -\pi_N^L  u)(x)\big|&\le \sum_{n=N+1}^{\infty}|\hat u_n^{L}| \le \sum_{n=N+1}^{\infty}\bigg\{\frac{1}{2^{\mu+m+1} \sqrt{\pi}}
	\frac{ (2n+1)\Gamma( ({n-\mu-m}+1)/ 2)}{\Gamma( ({n+\mu+m})/2+1)}
		\\
&\quad+\sum_{j=0}^{m}\widehat  C_{n,\mu+j}\bigg\}\, U^{(\mu,m)}_-.
	\end{split}\end{equation}
Similar to \eqref{FracLinfA-2}-\eqref{Tsigmabnd}, we have
\begin{equation}\begin{split}\label{FracEndLinfA-3}
\sum_{n=N+1}^{\infty}
	\frac{ (2n+1)\Gamma( ({n-\mu-m})/ 2)}{\Gamma( ({n+\mu+m+3})/2)} \le \frac{2}{ (\mu+m-1/2)}\frac{ \Gamma( ({N-\mu-m+1})/ 2)}{\Gamma( ({N+\mu+m})/2)}.
	\end{split}\end{equation}
Moreover, we can show that
\begin{equation}\begin{split}\label{FracEndLinfA-4}
\sum_{n=N+1}^{\infty}\widehat C_{n,\rho+1}&= \frac{2^{\rho+1}\Gamma(\rho+1)}{\pi} \sum_{n=N+1}^{\infty}\frac{(n+1/2)\Gamma(n-\rho)}{\Gamma(n+\rho+2)}\le\frac{2^{\rho+1}\Gamma(\rho+1)}{\pi} \sum_{n=N+1}^{\infty}\frac{\Gamma(n-\rho)}{\Gamma(n+\rho+1)}\\
&=\frac{2^{\rho}\Gamma(\rho+1)}{\pi(\rho-1)} \sum_{n=N+1}^{\infty}\bigg\{(n+\rho)\frac{\Gamma(n-\rho)}{\Gamma(n+\rho+1)}-(n-\rho)\frac{\Gamma(n-\rho)}{\Gamma(n+\rho+1)}\bigg\}
\\
&=\frac{2^{\rho}\Gamma(\rho+1)}{\pi(\rho-1)} \sum_{n=N+1}^{\infty}\bigg\{\frac{\Gamma(n-\rho)}{\Gamma(n+\rho)}-\frac{\Gamma(n-\rho+1)}{\Gamma(n+\rho+1)}\bigg\}
=\frac{2^{\rho}\Gamma(\rho+1)}{\pi(\rho-1)}\frac{\Gamma(N-\rho+1)}{\Gamma(N+\rho+1)}.
\end{split}\end{equation}
Thus the estimate \eqref{FracEndLinfA} is a consequence of  \eqref{FracEndLinfA-2}-\eqref{FracEndLinfA-3} and  \eqref{FracEndLinfA-4} with $\rho=\mu+j.$

(ii) We now turn to \eqref{FracEndL2B}.  By the orthogonality of Legendre polynomials, we find readily from
 \eqref{BoundForUnEnd} and the Cauchy-Schwarz inequality that
	\begin{equation}\label{BoundForUnEnd-4}
	\begin{split}
	\big\|u-\pi_N^{L} u\big\|_{L^2(\Omega)}^2=& \sum_{n=N+1}^{\infty}\frac{2}{2n+1}\big|\hat u_n^{L}\big|^2
	\le   \bigg\{ \sum_{n=N+1}^{\infty}
	 \frac{  (2n+1)\,\Gamma^2( ({n-\mu-m})/ 2)}{2^{2(\mu+m+1)} \,\Gamma^2( ({n+\mu+m}+3)/2)}\\
 & \quad +\sum_{n=N+1}^{\infty}\frac{4\,\widehat  C^2_{n,\mu}}{2n+1}\bigg(\sum_{j=0}^{m}\frac{\widehat  C_{n,\mu+j}}{\widehat  C_{n,\mu}}\bigg)^2\bigg\}\,( U^{(\mu,m)}_-)^2.
	\end{split}
	\end{equation}
	Following the derivations \eqref{FracL2-4}-\eqref{FracL2-5} with $\mu+m$ in place of $\mu$, we have
\begin{equation}\label{BoundForUnEnd-5}
	\begin{split}
\sum_{n=N+1}^{\infty}
	 \frac{ (2n+1)\,\Gamma^2( ({n-\mu-m})/ 2)}{\Gamma^2( ({n+\mu+m+3})/2)}\le \frac{2^{2\mu+2m+4}}{2\mu+2m+1}
\frac{ \Gamma( {N-\mu-m})}{\Gamma( {N+\mu+m+1})}.
	\end{split}
	\end{equation}
We next show that for $n\ge N+1,$
\begin{equation}\label{BoundForUnEnd-10}
	\begin{split}
	\sum_{j=0}^{m}\frac{\widehat  C_{n,\mu+j}}{\widehat  C_{n,\mu}}<\frac{2(N+1)}{2N+1}.
	\end{split}
	\end{equation}
	Observe that for  $1\le i\le j\le m$ and $ \mu+m+1\le N+1\le n, $
	\begin{equation*}\label{BoundForUnEnd-2}
	\begin{split}
(n-\mu-i)(n+\mu+i+1)&\ge (n-\mu-j)(n+\mu+j+1)\ge (n-\mu-m)(n+\mu+m+1)\\
& \ge (n-N)(n+N+1)\ge 2(N+1).
	\end{split}
	\end{equation*}
 Using this inequality and  \eqref{newconstants-2}, we obtain
	\begin{equation}\label{BoundForUnEnd-1}
	\begin{split}
	\sum_{j=0}^{m}\frac{\widehat  C_{n,\mu+j}}{\widehat  C_{n,\mu}}&=
\sum_{j=0}^{m}  \frac{\Gamma(n-\mu-j)}{\Gamma(n-\mu)}
\frac{\Gamma(n+\mu+2)}{\Gamma(n+\mu+j+2)}\\
&=1+\sum_{j=1}^{m}\bigg( \prod_{i=1}^j\frac{1}{(n-\mu-i)(n+\mu+i+1)}\bigg)\\
&= 1+\sum_{j=1}^{m}\Big(\frac{1}{2(N+1)}\Big)^j < \sum_{j=0}^{\infty}\Big(\frac{1}{2(N+1)}\Big)^j=\frac{2(N+1)}{2N+1}.
	\end{split}
	\end{equation}
This yields \eqref{BoundForUnEnd-10}.

Like  \eqref{FracL2-4}, we can get
\begin{equation*}\label{FracEndL2-3}
	\begin{split}
\frac{(n+1/2)\Gamma^2(n-\mu)}{\Gamma^2(n+\mu+2)}&\le\frac{2^{4\mu+3}}{4\mu+2}\bigg\{\frac{\Gamma(2n-2\mu-1)}{\Gamma(2n+2\mu+1)}-
\frac{\Gamma(2n-2\mu)}{\Gamma(2n+2\mu+2)}\bigg\}\\
&\le \frac{2^{4\mu+3}}{4\mu+2}\bigg\{\frac{\Gamma(2n-2\mu-1)}{\Gamma(2n+2\mu+1)}-
\frac{\Gamma(2n-2\mu+1)}{\Gamma(2n+2\mu+3)}\bigg\}.
	\end{split}
	\end{equation*}
Then similar to \eqref{FracL2-5}, we have
\begin{equation}\label{BoundForUnEnd-6}
	\begin{split}
\sum_{n=N+1}^{\infty} \frac{2}{2n+1}
	\widehat  C^2_{n,\mu}	&=\frac{2^{2\mu+2}\Gamma^2(\mu+1)}{\pi^2}\sum_{n=N+1}^{\infty} \!\!\frac{(n+1/2)\Gamma^2(n-\mu)}{\Gamma^2(n+\mu+2)}\\
&\le \frac{2^{6\mu+5}\,\Gamma^2(\mu+1)}{\pi^2(4\mu+2)}\frac{\Gamma(2N-2\mu+1)}{\Gamma(2N+2\mu+3)}.
	\end{split}
	\end{equation}

 Thus, we obtain \eqref{FracEndL2B} from \eqref{BoundForUnEnd-4}-\eqref{BoundForUnEnd-6}.
This completes the proof.
	\end{proof}
	\end{comment}

In contrast to  the interior singularity with a half-order loss, the $L^\infty$-estimate in this case  is optimal. In fact, for the endpoint singularity, the largest error occurs near the boundary where $|P_n(x)|$ attend its maximum  at $x=\pm 1$ (see
Figure \ref{PnWithWeight} (left)), so  the direct summation in e.g., \eqref{FracLinfA-2} will not overestimate.
As an illustration, we consider $u(x)=(1+x)^\mu$ with $\mu\in (k-1,k), k\in \mathbb N$.  From  Theorem \ref{TruncLegEnd} and  \eqref{NRatioGammaV2-1},
we find $\|u-\pi_N^L u\|_{L^\infty (\Omega)}\le C N^{-2\mu}$ and $\|u-\pi_N^L u\|_{L^2 (\Omega)}\le C N^{-2\mu-1}.$
	We tabulate in Table \ref{LinfSetaV2} the errors and convergence order of  Legendre approximations to  $u(x)=(x+1)^\mu$ with various $\mu$, which indicate  the optimal convergence order as predicted.

\begin{table}[!htbp]
	\centering
	\caption{Convergence  of Legendre expansion of  $(1+x)^\mu$.}
	\small
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
		\hline
		\multicolumn{1}{ |c  }{\multirow{2}{*}{$N$} } &
		\multicolumn{4}{ |c| }{Errors in $L^\infty$-norm} & \multicolumn{4}{ |c| }{Errors in $L^2$-norm}      \\ \cline{2-9}
		\multicolumn{1}{ |c  }{}                        &
		\multicolumn{1}{ |c| }{$\mu=0.1$}  & order
		&$ \mu=1.2$  & order  &  $\mu=0.1$  & order
		&$ \mu=1.2$  & order \\ \hline
$ 2^3 $ &6.15e-01 &--& 2.27e-3 &--&8.82e-3 &-- & 2.32e-04 &--  \\
$ 2^4 $ &5.41e-01 & 0.18 & 4.87e-04 & 2.22  &4.11e-03 & 1.10 & 2.64e-05 & 3.14 \\
$ 2^5 $ &4.74e-01 & 0.19 & 9.87e-05 & 2.30  &1.85e-03 & 1.15 & 2.75e-06 & 3.26 \\
$ 2^6 $ &4.14e-01 & 0.20 & 1.94e-05 & 2.35  &8.22e-04 & 1.17 & 2.74e-07 & 3.33 \\
$ 2^7 $ &3.61e-01 & 0.20 & 3.74e-06 & 2.37  &3.61e-04 & 1.19 & 2.67e-08 & 3.36 \\
$ 2^8 $ &3.15e-01 & 0.20 & 7.15e-07 & 2.39  &1.58e-04 & 1.19 & 2.56e-09 & 3.38 \\
		\hline
	\end{tabular} \label{LinfSetaV2}
\end{table}

\begin{comment}
In the end of this section, we conduct the analysis similar to proof of Theorem \ref{01estimate}, but for functions with endpoint singularities.  We demonstrate that the estimate \eqref{FracEndLinfA}  in Theorem \ref{TruncLegEnd} is optimal
for $u(x)=(1+x)^\alpha,$ which is applicable to  the typical singular function $u(x)=(1+x)^\alpha g(x)$ with smooth $g(x)$ and $g(-1)\not=0.$
	 \begin{prop}\label{unxEndSig} Consider $u(x)=(x+1)^\alpha$  on $\bar \Omega$ with $\alpha>0.$   Then we have the following bounds.
	\begin{itemize}
\item[(i)] For $N>\alpha-1,$
\begin{equation}\label{unxEndSig-1}
\begin{split}
\big|(u-\pi_N^L  u)(- 1)\big|&\le\frac{|\sin(\alpha\pi)|\,2^{\alpha}\,\Gamma(\alpha)\Gamma(\alpha+1)} \pi \frac{\Gamma(N-\alpha+1)}{\Gamma(N+\alpha+1)}.
\end{split}
\end{equation}
\item[(ii)] For $N>\alpha,$
\begin{equation}\label{unxEndSig-1+1}
\begin{split}
\big|(u-\pi_N^L  u)(1)\big|&\le\frac{|\sin(\alpha\pi)|\,2^{\alpha}\,\Gamma^2(\alpha+1)} \pi \frac{\Gamma(N-\alpha)}{\Gamma(N+\alpha+1)}.
\end{split}
\end{equation}
\item[(iii)] For $N>\alpha+3,$
 \begin{equation}\label{unxEndSig-0}
\begin{split}
\big|(u -\pi_N^L  u)(0)\big|\le  \frac{|\sin(\alpha\pi)|\, 2^{\alpha+1/2}\,\Gamma^2(\alpha+1)} {\pi^{3/2}} \frac{\Gamma(N-\alpha-3)}{\Gamma(N+\alpha-3/2)}.
\end{split}
\end{equation}
\end{itemize}
 \end{prop}
\begin{proof} The proof given in Appendix  \ref{AppendixD}.
\end{proof}

\begin{figure}[!ht]
	\begin{center}
		{~}\hspace*{-20pt}
\includegraphics[width=0.36\textwidth]{PoWisePosiEndSig}\qquad\quad
		\includegraphics[width=0.39\textwidth]{ErrorEndSigZerosPm1-1}
		\caption{$|(u -\pi_N^L  u)(x)|$  with $u(x)=(1+x)^{1.2}$ and different $N.$ Here, $\mathcal{B}_1,$ $\mathcal{B}_2,$ and $\mathcal{B}_3$ denote the upper bounds in  \eqref{unxEndSig-1}-\eqref{unxEndSig-0}, respectively.}
		\label{FigErrorEndSig}
	\end{center}
\end{figure}
Indeed,  Figure \ref{FigErrorEndSig} illustrates that the maximum point-wise error is attained at $x=-1,$ and the bounds in Proposition \ref{unxEndSig} and Theorem \ref{TruncLegEnd} agree well with the real decay of the errors.
 In fact, it is seen from  \eqref{unxEndSig-9} that   $\hat u_{n}^{L}\,P_{n}(- 1)$ has no sign change, so
$$|\hat u_{n}^{L}P_{n}(- 1)|=|\hat u_{n}^{L}|\max\limits_{|x|\le1}\{|P_n(x)|\}=|\hat u_{n}^{L}|,$$
which is the key to the optimality of the estimate.

\end{comment}

\subsection{Concluding remarks} We presented a new fractional Taylor formula for singular functions whose integer-order  derivatives
up to $k-1$ are absolutely continuous and Caputo fractional derivative of order $\mu\in (k-1,k]$  is of bounded variation. It could be viewed as an ``interpolation" between the  usual Taylor formulas of two consecutive integer orders.   We derived from this remarkable tool a similar fractional representation of the Legendre expansion of this type of functions, which became the cornerstone of the optimal  error estimates for  the Legendre orthogonal projection. The set of results under the fractional AC-BV framework greatly enriched the approximation theory for spectral and $hp$ methods. It set a good example to show how the fractional calculus could impact this classic field, and seamlessly bridge  between the results  valid only for integer cases.  Here we merely discussed the approximation results, but this will pave the way for the analysis of and applications to  singular problems, which  will be a topic  worthy of future deep investigation.

\begin{appendix}

\section{Useful properties of Gamma function}\label{Gmprop}
\renewcommand{\theequation}{A.\arabic{equation}}
\renewcommand{\thelmm}{A.\arabic{lmm}}

\setcounter{equation}{0}
\setcounter{lmm}{0}

Recall the Euler's reflection formula   (cf.\! \cite[Ch.2]{Andrews1999Book}):
\begin{equation}\label{nonitA}
\Gamma(1-a)\Gamma(a)=\frac{\pi} {\sin (\pi a)},\quad a\not=\pm 1,\pm 2,\cdots,
\end{equation}
and the Legendre duplication formula (cf.\! \cite[(5.5.5)]{Olver2010Book}):
	\begin{equation}\begin{split}\label{DuplicationFormula}
	\Gamma(2z)=\pi^{-1/2}2^{2z-1}\Gamma(z)\Gamma(z+1/2).
	\end{split}\end{equation}

From
	\cite[(1.1) and Thm. 10]{Alzer1997MC}, we have  that for $0\le a\le b$, the ratio
	\begin{equation}\label{gammratio}
	{\mathcal R}_b^a(z):=\frac{\Gamma(z + a)}{\Gamma(z + b)},\quad z\ge 0,
	\end{equation}
	is decreasing with respect to $z.$
On the other hand,   the ratio
	\begin{equation}\label{IneqtyLeg-1}
\widehat {\mathcal R}_c(z):=\frac{1}{\sqrt {z+c}}\frac{\Gamma(z+1)}{\Gamma(z+1/2)},
\end{equation}
	is increasing (resp. decreasing) on $[-1/2,\infty)$ (resp.  $(-c,\infty)$),  if $c\ge 1/2$ (resp. $c\le 1/4$),
based on   \cite[Corollary 2]{Bustoz1986MC}.

In the error bounds, the ratio of two Gamma functions appears very often, so the following inequality is useful.
\begin{lemma}\label{Gammalm} Let $b\in (a+m, a+m+1)$ for some  integer $m\ge 0,$ and
set  $b=a+m+\mu$ with $\mu\in (0,1).$
 Then for $z+a>0$ and $z+b>1, $ we have
\begin{equation}\label{NRatioGammaV2-1}
\frac 1 {(z+a)_{m}}\Big(z+b-\frac{3}{2}+\Big(\frac 5 4-\mu\Big)^{1 / 2}\Big)^{-\mu}< \frac{\Gamma(z+a)}{\Gamma(z+b)}< \frac 1 {(z+a)_{m}} \Big(z+b-\frac{\mu+1}{2}\Big)^{-\mu},
\end{equation}
where the Pochhammer symbol: $(c)_{m}=c(c+1)\cdots (c+m-1).$

\end{lemma}
\begin{proof} In fact,  \eqref{NRatioGammaV2-1} can be derived from the bounds in  \cite[(1.3)]{Kershaw83MC-Some}:

\begin{equation}\label{NRatioGamma-0}
\Big(x-\frac{1}{2}+\Big(\nu+\frac{1}{4}\Big)^{1 / 2}\Big)^{\nu-1}<\frac{\Gamma(x+\nu)}{\Gamma(x+1)}<\Big(x+\frac{\nu}{2}\Big)^{\nu-1},
\quad x> 0, \;\;  \nu \in (0,1).
\end{equation}
Indeed,  using the property  $\Gamma(z+1)=z\Gamma(z),$ we can write
\begin{equation*}
\frac{\Gamma(z+a)}{\Gamma(z+b)}= \frac 1 {(z+a)_{m}}\frac{\Gamma(z+a+m)}{\Gamma(z+b)}
= \frac 1 {(z+a)_{m}}\frac{\Gamma(z+b-\mu)}{\Gamma(z+b)}.
\end{equation*}
Then by  \eqref{NRatioGamma-0} with $x=z+b-1$ and $\nu=1-\mu,$  we obtain \eqref{NRatioGammaV2-1} immediately.
\end{proof}

\section{Proof of Lemma  \ref{FracIntPart-0} }\label{AppendixA}
\renewcommand{\theequation}{B.\arabic{equation}}

For $f\in L^1(\Omega)$ and $g\in {\rm AC}(\bar \Omega) $, changing the order of integration by the Fubini's Theorem, we derive from  \eqref{leftintRL}  that

\begin{equation*}\label{FracIntPart-3}
\begin{split}
& \int^{b}_{a}  f(x)  I_{a+}^{\rho} g'(x)\,{\rm d}x=\frac{1}{\Gamma(\rho)}\int_{a}^b  \bigg\{\int _{a}^x \frac{g'(y)}{(x-y)^{1-\rho}} {\rm d}y\bigg\}  f(x) \, {\rm d}x\\
&=\frac{1}{\Gamma(\rho)}\int_{a}^b  \bigg\{\int _{y}^b \frac{f(x)}{(x-y)^{1-\rho}} {\rm d}x\bigg\}  g'(y) \, {\rm d}y=\frac{1}{\Gamma(\rho)}\int_{a}^b  \bigg\{\int _{x}^b \frac{f(y)}{(y-x)^{1-\rho}} {\rm d}y\bigg\}  g'(x) \, {\rm d}x
\\
& =\int^{b}_{a} g'(x)\, I_{b-}^{\rho} f(x)\,{\rm d}x.

\end{split}
\end{equation*}

If $ I_{b-}^{\rho} f(x)\in {\rm BV}(\bar\Omega),$ we  derive from \eqref{IPPW1122} that
\begin{equation*}\label{FracIntPart-4}
\begin{split}
\int^{b}_{a}   f(x) & \,I_{a+}^{\rho} g'(x)\,{\rm d}x  =\int^{b}_{a} g'(x)\,  I_{b-}^{\rho} f(x)\,{\rm d}x = \big\{g(x)\,  I_{b-}^{\rho} f(x)\big\}\big|_{a^+}^{b^-} - \int^{b}_{a} g(x)
\,{{\rm d}\big\{I_{b-}^{\rho} f(x)\big\}}.
\end{split}
\end{equation*}
This yields \eqref{FracIntPart-1}.

We can derive  \eqref{FracIntPart-2} in a similar fashion.

\section{Proof of Proposition  \ref{boundaryvalue} }\label{AppendixA0}
\renewcommand{\theequation}{C.\arabic{equation}}

 Recall the first mean value theorem for the integral (cf.  \cite[p. 354]{Vladimir2016book}):  Let $f,g$ be Riemman integrable on $[c,d]$, $m=\inf\limits_{x\in [c,d]} f(x),$ and  $M=\sup\limits_{x\in [c,d]} f(x)$. If $g$ is nonnegative (or nonpositive) on $[c,d],$ then
	\begin{equation}\label{fgmeanv}
	\int_c^d f(x)g(x){\rm d}x=\kappa \int_c^d g(x){\rm d}x,\quad \kappa\in [m,M].
	\end{equation}
	Recall that (cf. \cite{Samko1993Book}):   for  $\alpha>-1$ and $\mu\in{\mathbb R}^+,$
	\begin{equation}\label{Lintformu}
	I_{a+}^\mu \, (x-a)^\alpha=  \dfrac{\Gamma(\alpha+1)}{\Gamma(\alpha+\mu+1)} (x-a)^{\alpha+\mu}.
	\end{equation}
	For any $x\in [a,a+\delta], $ we derive from \eqref{leftintRL}, \eqref{fgmeanv} and \eqref{Lintformu} that
	\begin{equation}\label{integbB}
	\begin{split}
	I_{a+}^\mu\, u (x) & =\frac 1 {\Gamma(\mu)}\int_{a}^x \frac{u(y)}{(x-y)^{1-\mu}} {\rm d}y= \frac{\kappa(x)}{\Gamma(\mu)} \int_{a}^x \frac{(y-a)^\alpha}{(x-y)^{1-\mu}} {\rm d}y\\
	&= \kappa(x){}_a I_{x}^\mu \, (x-a)^\alpha =\frac{  \Gamma(\alpha+1)\kappa(x)(x-a)^{\mu+\alpha}}{\Gamma(\alpha+\mu+1)} ,
	\end{split}\end{equation}
	where $\kappa(x)\in [m(x),M(x)]$, $m(x)=\inf\limits_{y\in [a,x]} v(y)$, $M(x)=\sup\limits_{y\in [a,x]} v(x)$. We know that
	\begin{equation}\label{integbC}
	\lim\limits_{x\rightarrow a^+}m(x)=v(a),\quad \lim\limits_{x\rightarrow a^+}M(x)=v(a)\Rightarrow \lim\limits_{x\rightarrow a^+}\kappa(x)=v(a).
	\end{equation}
	From \eqref{integbB} and \eqref{integbC}, we obtain \eqref{integbA}.

	This completes the proof.

\section{Proof of Theorem  \ref{01estimate} }\label{AppendixC}
\renewcommand{\theequation}{D.\arabic{equation}}
We start with the exact formula for the Legendre expansion coefficients  of $u(x)=|x|:$
\begin{equation}\label{un|x|}
\begin{split}
 \hat u_{2j}^L= \frac{(-1)^{j+1} (j+1/4) \Gamma(j-1/2)}{\sqrt{\pi}\, (j+1)!},\quad \hat u_{2j+1}^L=0,\quad j\ge 1,
\end{split}
\end{equation}
which can be derived from \eqref{HatUnCases1} with $k=1,$ i.e.,
\begin{equation*}
\begin{split}
\hat u_n^L&=\frac{2n+1}{2} \int^{1}_{-1}(I_{1-}^{2} P_n)(x)\, \rd \big\{ u^{(1)}(x)\big\}=\frac{2n+1}{2^{2}} (I_{1-}^{2} P_n)(0)\\
&
=\frac{2n+1}{2^{5}}
{}_2F_1\Big(\!\! -n+2, n+3;3;\frac{1} 2\Big),\quad n\ge 2,
\end{split}
\end{equation*}
and the value  at $z=1/2$ (cf. \cite[(15.4.28)]{Olver2010Book}):
\begin{align*}\label{Fatzero}
{}_2F_1\Big(a,b;\frac{a+b+1}2;\frac 1 2\Big) =
\frac{\sqrt{\pi}\,\Gamma((a+b+1)/2)}{\Gamma((a+1)/2)\Gamma((b+1)/2)}.
\end{align*}
Then we obtain from \eqref{un|x|} that
\begin{equation}\label{un|x|-4}
\begin{split}
(u -\pi_N^L  u)(0)&=\sum_{j=\lceil \frac {N+1}2 \rceil}^\infty \hat u_{2j}^{L}P_{2j}(0)=- \frac{1 }{\pi}
\sum_{j=\lceil \frac{N+1} 2 \rceil}^\infty  \frac{(j+1/4)\Gamma(j-1/2)\Gamma(j+1/2)} {\Gamma(j+1)\Gamma(j+2)},
\end{split}
\end{equation}
where $\lceil \frac{N+1} 2 \rceil$ is the  smallest integer   $\ge \frac{N+1} 2,$ and we used the known value (cf.
\cite{Szego1975Book}):

\begin{equation*}\label{P2kvalue}
P_{2j}(0)={}_2F_1\Big(\!\! -2j, 2j+1;1;\frac{1} 2\Big)=(-1)^j\frac{\Gamma( j+1/2)}{\sqrt{\pi}\,j!}.
\end{equation*}
From \eqref{gammratio}, we have
\begin{equation}\label{un|x|-5}
\frac{\Gamma(j+1/2)}{\Gamma(j+1)}\le \frac{\Gamma(j)}{\Gamma(j+1/2)}.
\end{equation}
Thus, using \eqref{un|x|-5} and $\Gamma(z+1)=z\Gamma(z),$ we obtain
\begin{equation*}\label{un|x|-6}
\begin{split}
\frac{(j+1/4)\Gamma(j-1/2)\Gamma(j+1/2)} {\Gamma(j+1)\Gamma(j+2)}&=\frac{j+1/4}{j+1} \frac{\Gamma(j-1/2)}{\Gamma(j+1)}\frac{\Gamma(j+1/2)} {\Gamma(j+1)} \le \frac{\Gamma(j-1/2)\Gamma(j)} {\Gamma(j+1/2)\Gamma(j+1)}\\
&=\frac{1}{(j-1/2)j}\le \frac{1}{(j-1)j}=\frac{1}{j-1}-\frac{1}{j} .
\end{split}
\end{equation*}
Then
\begin{equation*}\label{un|x|-7}
\begin{split}
&\sum_{j=\lceil \frac{N+1} 2 \rceil}^\infty  \frac{(j+1/4)\Gamma(j-1/2)\Gamma(j+1/2)} {\Gamma(j+1)\Gamma(j+2)}
\le \sum_{j=\lceil \frac {N+1}2 \rceil}^\infty \Big(\frac{1}{j-1}-\frac{1}{j}\Big) =\frac{1}{\lceil \frac{N+1} 2 \rceil-1}\le \frac{2}{N-1}.
\end{split}
\end{equation*}
From \eqref{un|x|-4} and the above, we get the first result in \eqref{un|x|-1+1}.

We  now prove the second estimate in \eqref{un|x|-1+1}.
As $P_n(\pm1)=(\pm1)^n,$ we derive from \eqref{un|x|} that

\begin{equation}\label{un|x|-0}
\begin{split}
(u -\pi_N^L  u)(\pm 1)&=\sum_{j=\lceil \frac {N+1} 2 \rceil}^\infty \hat u_{2j}^{L}\,P_{2j}(\pm1)=
\sum_{j=\lceil \frac {N+1} 2 \rceil}^\infty  \frac{(-1)^{j+1}}{\sqrt{\pi}}\frac{(j+1/4)\Gamma(j-1/2)} {\Gamma(j+2)}.
\end{split}
\end{equation}
Denoting
\begin{equation*}
\begin{split}
{\mathcal S}_j:=  \frac{(-1)^{j+1}}{\sqrt{\pi}}\frac{(j+1/4)\Gamma(j-1/2)} {\Gamma(j+2)},\quad {\mathcal T}_j:=
\frac{1 }{2\sqrt{\pi}}\frac{\Gamma(j-3/2)} {\Gamma(j)},
\end{split}
\end{equation*}
we have
\begin{equation}\label{un|x|V2-1}
\begin{split}
{\mathcal S}_j+{\mathcal S}_{j+1}&=(-1)^{j+1}\frac{3 }{2\sqrt{\pi}}\frac{(j+3/4)\Gamma(j-1/2)} {\Gamma(j+3)} \le \frac{3 }{2\sqrt{\pi}}\frac{\Gamma(j-1/2)} {\Gamma(j+2)}\\
&\le  \frac{3 }{4\sqrt{\pi}}\Big(\frac{\Gamma(j-3/2)} {\Gamma(j+1)}+\frac{\Gamma(j-1/2)} {\Gamma(j+2)}\Big)= \big({\mathcal T}_j-{\mathcal T}_{j+1}\big)+
\big({\mathcal T}_{j+1}-{\mathcal T}_{j+2}\big),
\end{split}
\end{equation}
where we noted
\begin{equation*}
\begin{split}
 \frac{\Gamma(j-1/2)} {\Gamma(j+2)}\le  \frac{\Gamma(j-3/2)} {\Gamma(j+1)},
\end{split}
\end{equation*}
and
\begin{equation*}
\begin{split}
 \frac{3 }{4\sqrt{\pi}}\frac{\Gamma(j-3/2)} {\Gamma(j+1)}= \frac{1 }{2\sqrt{\pi}}\bigg(j \frac{\Gamma(j-3/2)} {\Gamma(j+1)}-(j-3/2)\frac{\Gamma(j-3/2)} {\Gamma(j+1)}\bigg)={\mathcal T}_j-{\mathcal T}_{j+1}.
\end{split}
\end{equation*}
Thus from \eqref{gammratio} and \eqref{un|x|-0}-\eqref{un|x|V2-1}, we obtain
\begin{equation*}\label{un|x|-2}
\begin{split}
\big|(u -\pi_N^L  u)(\pm 1)\big|&= \big(|{\mathcal S}_{\lceil\frac{N+1} 2 \rceil}+{\mathcal S}_{\lceil\frac{N+1} 2 \rceil+1}|\big)+\cdots
+\big(|{\mathcal S}_{\lceil\frac{N+1} 2 \rceil+2i}+{\mathcal S}_{\lceil\frac{N+1} 2 \rceil+2i+1}|\big)+\cdots\\
 &\le \big\{\big({\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}}-{\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+1}\big)+
\big({\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+1}-{\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+2}\big)\big\}+\cdots\\
&\quad + \big\{\big({\mathcal T}_{{\lceil\frac{N+1} 2\rceil}+2i }-{\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+2i+1}\big)+
\big({\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+2i+1}-{\mathcal T}_{{\lceil\frac{N+1} 2 \rceil}+2i+2}\big)\big\}+\cdots\\
 &= \sum_{j=\lceil \frac{N+1} 2 \rceil}^\infty\big({\mathcal T}_j-{\mathcal T}_{j+1}\big)=\frac{1 }{2\sqrt{\pi}}\frac{\Gamma(\lceil \frac {N+1}2 \rceil-3/2)} {\Gamma(\lceil \frac{N+1}2 \rceil)}\le
\frac{1 }{2\sqrt{\pi}}\frac{\Gamma( N/2-1)} {\Gamma( N/2+1/2)}.
\end{split}
\end{equation*}
This ends the proof.

\end{appendix}