feat: implement more powers for starting weights

docs/references.rst

@@ -22,6 +22,11 @@ References
     Computer Methods in Applied Mechanics and Engineering, Vol. 194, pp. 743-773, 2005,
     `DOI <https://doi.org/10.1016/j.cma.2004.06.006>`__.
 
+.. [Diethelm2006] K. Diethelm, J. M. Ford, N. J. Ford, M. Weilbeer,
+    *Pitfalls in Fast Numerical Solvers for Fractional Differential Equations*,
+    Journal of Computational and Applied Mathematics, Vol. 186, pp. 482--503, 2006,
+    `DOI <https://doi.org/10.1016/j.cam.2005.03.023>`__.
+
 .. [Shen2011] J. Shen, T. Tang, L.-L. Wang,
     *Spectral Methods - Algorithms, Analysis and Applications*,
     Springer, 2011.

pycaputo/generating_functions.py

@@ -9,44 +9,154 @@
 
 from pycaputo.utils import Array
 
-# {{{ Lubich1986 weights
+# {{{ LMM: Starting Weights
 
 
-def lubich_bdf_starting_weights_count(order: int, alpha: float) -> int:
-    r"""An estimate for the number of starting weights from [Lubich1986]_.
+def lmm_starting_weights(
+    w: Array, sigma: Array, alpha: float, *, atol: float | None = None
+) -> Iterator[Array]:
+    r"""Constructs starting weights for a given set of weights *w* of a
+    fractional-order linear multistep method (FLMM).
+
+    The starting weights are introduced to handle a lack of smoothness of the
+    function :math:`f(x)` being integrated. They are constructed in such a way
+    that they are exact for a series of monomials, which results in an accurate
+    quadrature for functions of the form :math:`f(x) = f_i x^{\sigma_i} +
+    x^{\beta} g(x)`, where :math:`g` is smooth and :math:`\beta < p`
+    (see Remark 3.6 in [Li2020]_). This ensures that the Taylor expansion near
+    the origin is sufficiently accurate.
+
+    Therefore, they are obtained by solving
+
+    .. math::
+
+        \sum_{k = 0}^s w_{mk} k^{\sigma_q} =
+        \frac{\Gamma(\sigma_q +1)}{\Gamma(\sigma_q + \alpha + 1)}
+        m^{\sigma_q + \alpha} -
+        \sum_{k = 0}^s w_{n - k} j^{\sigma_q}
+
+    for a set of :math:`\sigma_q` powers. In the simplest case, we can just set
+    `sigma \in \{0, 1, \dots, p - 1\}` and obtain integer powers. Other values
+    can be chosen depending on the behaviour of :math:`f(x)` near the origin.
+
+    Note that the starting weights are only computed for :math:`m \ge s`.
+    The initial :math:`s` steps are expected to be computed in some other
+    fashion.
+
+    :arg w: convolution weights of an FMM.
+    :arg sigma: an array of monomial powers for which the starting weights are exact.
+    :arg alpha: order of the fractional derivative to approximate.
+    :arg atol: absolute tolerance used for the GMRES solver. If *None*, an
+        exact LU-based solver is used instead (see Section 3.2 in [Diethelm2006]_
+        for a discussion of these methods).
+
+    :returns: the starting weights for every point :math:`x_m` starting with
+        :math:`m \ge s`.
+    """
+
+    from scipy.linalg import lu_factor, lu_solve
+    from scipy.sparse.linalg import gmres
+    from scipy.special import gamma
+
+    s = sigma.size
+    j = np.arange(1, s + 1).reshape(-1, 1)
+
+    A = j**sigma
+    assert A.shape == (s, s)
+
+    if atol is None:
+        lu_p = lu_factor(A)
+
+    for k in range(s, w.size):
+        b = (
+            gamma(sigma + 1) / gamma(sigma + alpha + 1) * k ** (sigma + alpha)
+            - A @ w[k - s : k][::-1]
+        )
+
+        if atol is None:
+            omega = lu_solve(lu_p, b)
+        else:
+            omega, _ = gmres(A, b, atol=atol)
+
+        yield omega
+
+
+def diethelm_bdf_starting_powers(order: int, alpha: float) -> Array:
+    r"""Generate monomial powers for the starting weights from [Diethelm2006]_.
+
+    The construction of the starting weights is given in Section 3.1 from
+    [Diethelm2006]_ and proposes using the powers
+
+    .. math::
+
+        \sigma \in \{i + \alpha j \mid i, j \ge 0, i + \alpha j \le p - 1 \},
+
+    where :math:`p` is the desired order. For certain values of :math:`\alpha`,
+    these monomials can repeat, e.g. for :math:`\alpha = 0.1` we get the same
+    value for :math:`(i, j) = (1, 0)` and :math:`(i, j) = (0, 10)`. This function
+    returns only unique powers.
 
     :arg order: order of the BDF method.
     :arg alpha: order of the fractional derivative.
     """
     from math import ceil
 
-    return ceil(order * (order - 1 + 2 * alpha) / (2 * abs(alpha)))
+    i, j = np.array(
+        [[i, j] for i in range(order) for j in range(ceil((order - i - 1) / alpha) + 1)]
+    ).T
+
+    return np.array(np.unique(i + alpha * j))
+
 
+def lubich_bdf_starting_powers(order: int, alpha: float, *, beta: float = 1.0) -> Array:
+    r"""Generate monomial powers for the starting weights from [Lubich1986]_.
 
-def lubich_bdf_starting_powers(order: int, alpha: float) -> Array:
-    r"""Generate monomial powers for the starting weights in [Lubich1986]_.
+    The construction of the starting weights is given in Section 4.2 of
+    [Lubich1986]_ and proposes
+
+    .. math::
 
-    In general, we want all monomials :math:`t^\gamma` such that
-    :math:`\gamma = i + \alpha j \le p - 1`, where :math:`i, j \in \mathbb{N}`
-    are non-negative integers.
+        \sigma \in \{i + \beta - 1 \mid q \le p - 1\},
 
-    For certain values of :math:`\alpha`, these monomials can repeat, e.g. for
-    :math:`\alpha = 0.1` we get the same value for :math:`(i, j) = (1, 0)` and
-    :math:`(i, j) = (0, 10)`. This function returns all the powers, so that the
-    caller can decide to use :func:`numpy.unique` or similar functions to
-    remove duplicates.
+    where :math:`\beta` is chosen such that :math:`q + \beta - 1\ le p - 1`
+    and :math:`q + \alpha + \beta - 1 < p`, according to Theorem 2.4 from
+    [Lubich1986]_. In general multiple such :math:`\beta` exist and choosing
+    more and prove beneficial.
 
     :arg order: order of the BDF method.
     :arg alpha: order of the fractional derivative.
-    :returns: an array of ``(i, j)`` indices of shape ``(2, nterms)``.
     """
-    from math import ceil
+    from math import floor
+
+    # NOTE: trying to satisfy
+    # 0 <= q <= p - \beta and 0 <= q < p + 1 - alpha - beta
+    qmax = floor(max(0, min(order - beta, order - alpha - beta + 1))) + 1
 
-    return np.array([
-        [i, j]
-        for i in range(order)
-        for j in range(ceil((order - i - 1) / alpha) + 1)
-        ]).T
+    return np.array([q + beta - 1 for q in range(qmax)])
+
+
+def garrappa_bdf_starting_powers(order: int, alpha: float) -> Array:
+    r"""Generate monomial powers for the starting weights from
+    `FLMM2 <https://www.mathworks.com/matlabcentral/fileexchange/47081-flmm2>`__.
+
+    :arg order: order of the BDF method.
+    :arg alpha: order of the fractional derivative.
+    """
+    from math import ceil, floor
+
+    if order == 2:
+        k = floor(1 / abs(alpha))
+    else:
+        # NOTE: this should be vaguely the cardinality of `diethelm_bdf_starting_powers`
+        k = ceil(order * (order - 1 + 2 * alpha) / (2 * abs(alpha)))
+
+    return np.arange(k) * alpha
+
+
+# }}}
+
+
+# {{{ Lubich1986: BDF Convolution Weights
 
 
 def lubich_bdf_weights(alpha: float, order: int, n: int) -> Array:
@@ -77,7 +187,8 @@ def lubich_bdf_weights(alpha: float, order: int, n: int) -> Array:
 
     :arg alpha: power of the generating function.
     :arg order: order of the method, only :math:`1 \le p \le 6` is supported.
-    :arg n: number of weights.
+    :arg n: number of weights, which corresponds to the number of discretization
+        intervals.
     """
     if order <= 0:
         raise ValueError(f"Negative orders are not supported: {order}")
@@ -128,72 +239,4 @@ def lubich_bdf_weights(alpha: float, order: int, n: int) -> Array:
     return w
 
 
-def lubich_bdf_starting_weights(
-    w: Array, sigma: Array, alpha: float, *, atol: float | None = None
-) -> Iterator[Array]:
-    r"""Constructs starting weights for a given set of weights *w* from [Lubich1986]_.
-
-    The starting weights are introduced to handle a lack of smoothness of the
-    function :math:`f(x)` being integrated. They are constructed in such a way
-    that they are exact for a series of monomials, which results in an accurate
-    quadrature for functions of the form :math:`f(x) = f_i x^{\sigma_i} +
-    x^{\beta} g(x)`, where :math:`g` is smooth and :math:`\beta < p`
-    (see Remark 3.6 in [Li2020]_). This ensures that the Taylor expansion near
-    the origin is sufficiently accurate.
-
-    Therefore, they are obtained by solving
-
-    .. math::
-
-        \sum_{k = 0}^s w_{mk} k^{\sigma_q} =
-        \frac{\Gamma(\sigma_q +1)}{\Gamma(\sigma_q + \alpha + 1)}
-        m^{\sigma_q + \alpha} -
-        \sum_{k = 0}^s w_{n - k} j^{\sigma_q}
-
-    for :math:`q \in \{0, 1, \dots, s - 1\}`, where :math:`s` is the size of *sigma*.
-    In the simplest case, we can just set `sigma \in \{0, 1, \dots, p - 1\}`
-    and obtain integer powers. Other values can be chosen depending on the
-    behaviour of :math:`f(x)` near the origin.
-
-    Note that the starting weights are only computed for :math:`m \ge s`.
-    The initial :math:`s` steps are expected to be computed in some other
-    fashion.
-
-    :arg w: convolution weights of order :math:`p` defined by [Lubich1986]_.
-    :arg sigma: an array of monomial powers for which the starting weights are exact.
-    :arg alpha: order of the fractional derivative to approximate.
-    :arg atol: absolute tolerance used for the GMRES solver. If *None*, an
-        exact LU-based solver is used instead.
-
-    :returns: the starting weights for every point :math:`x_m` starting with
-        :math:`m \ge s`.
-    """
-
-    from scipy.linalg import lu_factor, lu_solve
-    from scipy.sparse.linalg import gmres
-    from scipy.special import gamma
-
-    s = sigma.size
-    j = np.arange(1, s + 1).reshape(-1, 1)
-
-    A = j**sigma
-    assert A.shape == (s, s)
-
-    if atol is None:
-        lu_p = lu_factor(A)
-
-    for k in range(s, w.size):
-        b = (
-            gamma(sigma + 1) / gamma(sigma + alpha + 1) * k ** (sigma + alpha)
-            - A @ w[k - s : k][::-1]
-        )
-
-        if atol is None:
-            omega = lu_solve(lu_p, b)
-        else:
-            omega, _ = gmres(A, b, atol=atol)
-
-        yield omega
-
-
 # }}}

pycaputo/quadrature/riemann_liouville.py

@@ -491,8 +491,8 @@ def _quad_rl_conv(
         raise TypeError(f"Only uniforms points are supported: '{type(p).__name__}'")
 
     from pycaputo.generating_functions import (
-        lubich_bdf_starting_weights,
-        lubich_bdf_starting_weights_count,
+        lmm_starting_weights,
+        lubich_bdf_starting_powers,
         lubich_bdf_weights,
     )
 
@@ -506,10 +506,10 @@ def _quad_rl_conv(
     w = lubich_bdf_weights(-alpha, m.quad_order, p.n)
 
     if np.isfinite(m.beta):
-        s = lubich_bdf_starting_weights_count(m.quad_order, alpha)
-        sigma = np.arange(s)
-        omegas = lubich_bdf_starting_weights(w, sigma, alpha)
+        sigma = lubich_bdf_starting_powers(m.quad_order, alpha, beta=m.beta)
+        omegas = lmm_starting_weights(w, sigma, alpha)
 
+        s = sigma.size
         for n, omega in enumerate(omegas):
             qc = np.sum(w[: n + s][::-1] * fx[: n + s])
             qs = np.sum(omega * fx[: s + 1])

tests/test_generating_functions.py

@@ -22,9 +22,7 @@
 def test_lubich_bdf_weights(
     order: int, alpha: float, *, visualize: bool = False
 ) -> None:
-    from pycaputo.generating_functions import (
-        lubich_bdf_weights,
-    )
+    from pycaputo.generating_functions import lubich_bdf_weights
     from pycaputo.utils import EOCRecorder, savefig
 
     eoc = EOCRecorder(order=order)