feat(generating_functions): generalize starting weights

pycaputo/generating_functions.py

@@ -129,72 +129,67 @@ def lubich_bdf_weights(alpha: float, order: int, n: int) -> Array:
 
 
 def lubich_bdf_starting_weights(
-    w: Array, s: int, alpha: float, *, beta: float = 1.0, atol: float | None = None
+    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
-    functions 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) = x^{\beta - 1} g(x)`,
-    where :math:`g` is smooth.
+    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^{q + \beta - 1} =
-        \frac{\Gamma(q + \beta)}{\Gamma(q + \beta + \alpha)}
-        m^{q + \beta + \alpha - 1} -
-        \sum_{k = 0}^s w_{n - k} j^{q + \beta - 1}
+        \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}
 
-    where :math:`q \in \{0, 1, \dots, s - 1\}` and :math:`\beta \ne \{0, -1, \dots\}`.
-    In the simplest case, we can take :math:`\beta = 1` and obtain integer
-    powers. Other values can be chosen depending on the behaviour of :math:`f(x)`
-    near the origin.
+    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 defined by [Lubich1986]_.
-    :arg s: number of starting weights.
+    :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 beta: order of the singularity at the origin.
     :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`.
     """
 
-    if s <= 0:
-        raise ValueError(f"Negative s is not allowed: {s}")
-
-    if beta.is_integer() and beta <= 0:
-        raise ValueError(f"Values of beta in 0, -1, ... are not supported: {beta}")
-
     from scipy.linalg import lu_factor, lu_solve
     from scipy.sparse.linalg import gmres
     from scipy.special import gamma
 
-    q = np.arange(s) + beta - 1
+    s = sigma.size
     j = np.arange(1, s + 1).reshape(-1, 1)
 
-    A = j**q
+    A = j**sigma
     assert A.shape == (s, s)
 
     if atol is None:
-        lu, p = lu_factor(A)
+        lu_p = lu_factor(A)
 
     for k in range(s, w.size):
         b = (
-            gamma(q + 1) / gamma(q + alpha + 1) * k ** (q + alpha)
+            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)
+            omega = lu_solve(lu_p, b)
         else:
             omega, _ = gmres(A, b, atol=atol)
 

pycaputo/quadrature.py

@@ -545,7 +545,8 @@ def _quad_rl_conv(
 
     if np.isfinite(m.beta):
         s = lubich_bdf_starting_weights_count(m.quad_order, alpha)
-        omegas = lubich_bdf_starting_weights(w, s, alpha, beta=m.beta)
+        sigma = np.arange(s)
+        omegas = lubich_bdf_starting_weights(w, sigma, alpha)
 
         for n, omega in enumerate(omegas):
             qc = np.sum(w[: n + s][::-1] * fx[: n + s])