Uniform-grid and graded-grid L1 discretization

Input $0<\alpha <1$, $T>0$, an integer $N\ge 1$, $\lambda \ge 0$, an initial value $U^0=u_0$, and forcing values $F^1,\dots ,F^N$. Choose either the uniform grid $t_k=kh$ with $h=T\mathrm{/}N$, or a strictly increasing grid $0=t_0<t_1<\cdots <t_N=T$. A power-graded choice is $t_k=T(k\mathrm{/}N{)}^r$ with a specified $r\ge 1$. Write $h_k=t_k-t_{k-1}$ and, for $1\le k\le n\le N$, define the exact kernel weight

$b_{n,k}:=\frac{(t_n-t_{k-1}{)}^{1-\alpha }-(t_n-t_k{)}^{1-\alpha }}{\mathrm{\Gamma }(2-\alpha )h_k}\mathrm{.}$

For a grid sequence $V^0,\dots ,V^N$, the L1 approximation to the Caputo derivative at $t_n$ is

${\delta }_{\mathcal{T}}^{\alpha }{V}^n:={\sum }_{k=1}^n{b}_{n,k}(V^k-V^{k-1})\mathrm{.}$

For $n=1,\dots ,N$ in increasing order, compute the value given by the scalar recurrence

$U^n:=\frac{F^n+b_{n,1}{U}^0+{\sum }_{k=1}^{n-1}(b_{n,k+1}-b_{n,k})U^k}{b_{n,n}+\lambda }\mathrm{.}$

Return $U^0,\dots ,U^N$. On a uniform grid, writing $a_j=(j+1{)}^{1-\alpha }-j^{1-\alpha }$ gives

${\delta }_h^{\alpha }{U}^n=\frac{h^{-\alpha }}{\mathrm{\Gamma }(2-\alpha )}{\sum }_{k=1}^n{a}_{n-k}(U^k-U^{k-1})\mathrm{.}$

An input outside the stated parameter or grid domain returns invalid-input. A finite-precision implementation that cannot certify a scalar solve returns solve-failed without reporting a value. The graded-grid procedure is an implementation definition only; no optimal graded-mesh rate is part of this item.