From 2e6efce0d929d0a18132d64f2d97757e9c5a10ae Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 21 Sep 2023 12:00:43 +0300 Subject: [PATCH] feat: implement more powers for starting weights --- docs/references.rst | 5 + pycaputo/generating_functions.py | 223 ++++++++++++++--------- pycaputo/quadrature/riemann_liouville.py | 10 +- tests/test_generating_functions.py | 4 +- 4 files changed, 144 insertions(+), 98 deletions(-) diff --git a/docs/references.rst b/docs/references.rst index 33eb39c..9588faf 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -22,6 +22,11 @@ References Computer Methods in Applied Mechanics and Engineering, Vol. 194, pp. 743-773, 2005, `DOI `__. +.. [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 `__. + .. [Shen2011] J. Shen, T. Tang, L.-L. Wang, *Spectral Methods - Algorithms, Analysis and Applications*, Springer, 2011. diff --git a/pycaputo/generating_functions.py b/pycaputo/generating_functions.py index 87aedf9..9579270 100644 --- a/pycaputo/generating_functions.py +++ b/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 `__. + + :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 - - # }}} diff --git a/pycaputo/quadrature/riemann_liouville.py b/pycaputo/quadrature/riemann_liouville.py index 991154e..f1ea1ed 100644 --- a/pycaputo/quadrature/riemann_liouville.py +++ b/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]) diff --git a/tests/test_generating_functions.py b/tests/test_generating_functions.py index dd0af9a..ef898ab 100644 --- a/tests/test_generating_functions.py +++ b/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)