from __future__ import annotations
from fractions import Fraction
from functools import cache
import numpy as np
import scipy.linalg as sla
import sympy
[docs]
@cache
def gregory_coef(k):
r"""Gregory coefficients.
See https://en.wikipedia.org/wiki/Gregory_coefficients.
Starting from :math:`G_0`, the first few coefficients are:
:math:`-\frac{1}{2}, \frac{1}{12}, -\frac{1}{24}, \frac{19}{720}, \ldots`
.. math::
G_k = \frac{-1}{(k+1)!} \left[ \frac{d^{k+1}}{dz^{k+1}}\frac{z}{\ln(1+z)}\right]_{z=0}
Parameters
----------
k : int
(starts from zero)
Returns
-------
:class:`sympy.core.numbers.Rational`
kth Gregory coefficient
"""
assert k >= 0
if k == 0:
sympy.sympify(Fraction(-1, 2))
grn = sympy.sympify(Fraction(1, k + 2))
for r in range(k):
grn -= abs(gregory_coef(r)) * Fraction(1, k + 1 - r)
return grn * (-1) ** (k + 1)
[docs]
@cache
def gregory_weights_exact(k):
r"""Gregory quadrature weights (exact).
See Table 9 in [NESSI]_.
.. math::
\mathbf{w} =
\begin{bmatrix}
1 & 0 & 0 & 0 & \cdots \\
-1 & 1 & 0 & 0 & \cdots \\
1 & -2 & 1 & 0 & \cdots \\
-1 & 3 & -3 & 1 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots \\
\end{bmatrix}
\begin{bmatrix}
G_0 \\
G_1 \\
G_2 \\
G_3 \\
\vdots
\end{bmatrix}
+ \begin{bmatrix}
1 \\
1 \\
1 \\
1 \\
\vdots
\end{bmatrix}
Parameters
----------
k : int
Starts from 0; gives degree of polynomial approx.
Returns
-------
w : (k+1,) ndarray
Order k Gregory quadrature weights, exact.
"""
gr_coef = np.array([gregory_coef(r) for r in range(k + 1)])
wts_exact = sla.invpascal(k + 1, kind="upper", exact=True) @ gr_coef + 1
wts_exact.setflags(write=False)
return wts_exact
[docs]
@cache
def gregory_weights(k):
"""Gregory quadrature weights (double precision).
See Table 9 in [NESSI]_.
Parameters
----------
k : int
Starts from 0; gives degree of polynomial approx.
Returns
-------
w : (k+1,) ndarray
Order n Gregory quadrature weights (double precision)
"""
float_weights = np.array(gregory_weights_exact(k), dtype=np.float64)
float_weights.setflags(write=False)
return float_weights
[docs]
@cache
def sigma_weights_exact(k):
"""Sigma quadrature weights (exact).
Parameters
----------
k : int
Starts from 0; gives degree of polynomial approx.
Returns
-------
sigma : (k+1, 2k+2) ndarray
Order k sigma quadrature weights, exact.
"""
omega = gregory_weights_exact(k)
sigma = np.zeros((k + 1, 2 * k + 2), dtype=object)
for i in range(1, k + 2):
sigma[i - 1, : k + 1] += omega
sigma[i - 1, i : k + 1] += np.flip(omega[i:]) - 1
sigma[i - 1, k + 1 : k + 1 + i] = np.flip(omega[:i])
sigma.setflags(write=False)
return sigma
[docs]
@cache
def sigma_weights(k):
"""Sigma quadrature weights.
Parameters
----------
k : int
Starts from 0; gives degree of polynomial approx.
Returns
-------
sigma : (k+1, 2k+2) ndarray
Order k sigma quadrature weights, double precision.
"""
sigma = np.array(sigma_weights_exact(k), dtype=np.float64)
sigma.setflags(write=False)
return sigma
[docs]
@cache
def sympy_vandermonde(k):
r"""Vandermonde matrix.
.. math::
V^{(k)} = \begin{bmatrix}
1 & 0 & 0 & \cdots & 0 \\
1 & 1 & 1 & \cdots & 1 \\
1 & 2^1 & 2^2 & \cdots & 2^n \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & k & k^2 & \cdots & k^k \\
\end{bmatrix}
Parameters
----------
k : int
Polynomial degree.
Returns
-------
v : (k+1, k+1) :class:`sympy.matrices.dense.DenseMatrix`
Vandermonde matrix constructed from points 0, 1, ..., k
"""
return (
sympy.sympify(np.vander(range(k + 1), increasing=True))
.tomatrix()
.as_immutable()
)
[docs]
@cache
def sympy_inv_vandermonde(k):
"""Inverse of Vandermonde matrix.
Parameters
----------
k : int
Polynomial degree.
Returns
-------
vinv : (k+1, k+1) :class:`sympy.matrices.dense.DenseMatrix`
Inverse of Vandermonde matrix constructed from points 0, 1, ..., n
"""
return sympy_vandermonde(k).inv().as_immutable()
[docs]
@cache
def poly_interpolation_weights_exact(k):
"""Polynomial interpolation weights (exact).
.. math::
P^{(k)} = (V^{(k)})^{-1}
Parameters
----------
k : int
Polynomial degree.
Returns
-------
P : (k+1,k+1) ndarray
"""
vinv = sympy.matrix2numpy(sympy_inv_vandermonde(k))
vinv.setflags(write=False)
return vinv
[docs]
@cache
def poly_interpolation_weights(k):
"""Double precision polynomial interpolation weights.
Same as `poly_interpolation_weights_exact` but in double precision.
Parameters
----------
k : int
Polynomial degree.
Returns
-------
P : (k+1,k+1) ndarray
"""
wts = np.array(poly_interpolation_weights_exact(k), dtype=np.float64)
wts.setflags(write=False)
return wts
[docs]
@cache
def poly_integration_weights_exact(k):
r"""Polynomial integration weights.
.. math::
I^{(k)} = \begin{bmatrix}
0 & 0 & 0 & \cdots & 0 \\
\frac{1}{1} & \frac{1}{2} & \frac{1}{3} & \cdots & \frac{1}{k+1} \\
\frac{2^1}{1} & \frac{2^2}{2} & \frac{2^3}{3} & \cdots & \frac{2^{k+1}}{k+1} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\frac{k^1}{1} & \frac{k^2}{2} & \frac{k^3}{3} & \cdots & \frac{k^{k+1}}{k+1} \\
\end{bmatrix} (V^{(k)})^{-1}
Parameters
----------
k : int
Polynomial degree.
Returns
-------
I : (k+1,k+1) ndarray
Polynomial integration weights.
"""
aux = sympy.Matrix(
[[Fraction(r ** (j + 1), j + 1) for j in range(k + 1)] for r in range(k + 1)]
)
return np.asarray(aux * sympy_inv_vandermonde(k))
[docs]
@cache
def poly_integration_weights(k):
"""Double precision polynomial integration weights.
Same as `poly_integration_weights_exact` but in double precision.
Parameters
----------
k : int
Polynomial degree.
Returns
-------
I : (k+1,k+1) ndarray
Polynomial integration weights.
"""
wts = np.array(poly_integration_weights_exact(k), dtype=np.float64)
wts.setflags(write=False)
return wts
[docs]
@cache
def poly_differentiation_weights_exact(k):
r"""Polynomial differentiation weights.
See equation 85 in [NESSI]_.
.. math::
D^{(k)}_{ml} = \sum_{a=1}^k P^{(k)}_{al}a m^{a-1}.
Parameters
----------
k : int
Polynomial degree.
Returns
-------
D : (k+1,k+1) ndarray
Polynomial differentiation weights.
"""
aux = sympy.Matrix(
[[a * r ** (a - 1) for a in range(1, k + 1)] for r in range(k + 1)]
)
wts = sympy.matrix2numpy(aux * sympy_inv_vandermonde(k)[1:, :])
wts.setflags(write=False)
return wts
[docs]
@cache
def poly_differentiation_weights(k):
"""Double precision polynomial differentiation weights.
Same as `poly_differentiation_weights_exact` but in double precision.
Parameters
----------
k : int
Polynomial degree.
Returns
-------
D : (k+1,k+1) ndarray
Polynomial differentiation weights.
"""
wts = np.array(poly_differentiation_weights_exact(k), dtype=np.float64)
wts.setflags(write=False)
return wts
[docs]
@cache
def poly_backdifferentiation_weights_exact(n):
return -poly_differentiation_weights_exact(n)[0, :]
[docs]
@cache
def poly_backdifferentiation_weights(n, rev=True):
wts = np.array(poly_backdifferentiation_weights_exact(n), dtype=np.float64)
if rev:
wts = np.ascontiguousarray(np.flip(wts))
wts.setflags(write=False)
return wts
[docs]
def master_gregory_weights(k, n):
"""Quadrature weights.
Gives row n of the matrix in Eq. 98 from [NESSI]_.
Parameters
----------
k : int
Polynomial degree.
n : int
Row (or upper limit of integration).
Returns
-------
w : (n+1,) ndarray
Quadrature weights.
"""
if n <= k:
return poly_integration_weights(k)[n]
if n <= 2 * k:
return sigma_weights(k)[n - k - 1]
gwt = gregory_weights(k)
wt = np.ones(n + 1, dtype=np.float64)
wt[: k + 1] = gwt
wt[-k - 1 :] = gwt[::-1]
return wt
[docs]
@cache
def conv_integral_r_exact(m, k):
r"""Equation 104 [NESSI]_.
.. math::
R^{(k)}_{m;r,s} = \sum_{a,b=0}^k P^{(k)}_{ar} P^{(k)}_{bs} \int_0^m (m-x)^a x^b \, dx
Parameters
----------
m : int
Nonnegative
k : int
Nonnegative
Returns
-------
R : (k+1, k+1) ndarray
Equation 104 weights for fixed m, k
"""
vinv = sympy.matrix2numpy(sympy_inv_vandermonde(k))
arr = integral_intermediate2(m, k)
ret = vinv.T @ arr @ vinv
ret.setflags(write=False)
return ret
[docs]
@cache
def conv_integral_r(m, k):
r"""Double-precision version of `conv_integral_r_exact`.
Parameters
----------
m : int
Nonnegative
k : int
Nonnegative
Returns
-------
R : (k+1, k+1) ndarray
Equation 104 weights for fixed m, k
See Also
--------
conv_integral_r_exact : Exact version
"""
ret = np.array(conv_integral_r_exact(m, k), dtype=np.float64)
ret.setflags(write=False)
return ret