pytrunc.utils module#

Numerical utilities.

This module provides the numerical routines used by the phase and truncation modules: the Legendre polynomials and their derivatives, the Bessel functions of the first kind, and the Lobatto quadrature abscissas, weights and integration.

pytrunc.utils.bessel_j1(x: ArrayLike, acc: float = 1e-08, max_iter: int = 50) → ndarray[tuple[Any, ...], dtype[float64]][source]#

The Bessel first kind function J1(x) of order 1

Parameters:

xndarray: The variable x of J1(x), it must be 1-D
accfloat, optional: The tolerance for numerical errors. Default is 1e-8
max_iterint, optional: The maximum number of iterations trying to improve the error accuracy. Default is 50

Returns:

ndarray: The values of the Bessel function J1(x), a 1-D ndarray with the same shape as x

Notes

The scipy equivalent -> scipy.special.j1(x)

Examples

>>> import numpy as np
>>> from pytrunc.utils import bessel_j1
>>> x = np.array([0.0, 1.0, 2.0])
>>> bessel_j1(x)
array([0.        , 0.44005059, 0.57672481])

pytrunc.utils.bessel_j1_derivative(x: ArrayLike, acc: float = 1e-08, max_iter: int = 50) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute the Bessel first kind derivative of order 1 d(J1(x))

Parameters:

xndarray: The variable x of d(J1(x)), it must be 1-D
accfloat, optional: The tolerance for numerical errors. Default is 1e-8
max_iterint, optional: The maximum number of iterations trying to improve the error accuracy. Default is 50

Returns:

ndarray: The values of the Bessel function derivative d(J1(x)), a 1-D ndarray with the same shape as x

Notes

The scipy equivalent -> scipy.special.jvp(1,x)

Examples

>>> import numpy as np
>>> from pytrunc.utils import bessel_j1_derivative
>>> x = np.array([0.0, 1.0, 2.0])
>>> bessel_j1_derivative(x)
array([ 0.5       ,  0.3251471 , -0.06447162])

pytrunc.utils.bessel_j1_roots(nb_roots: int, acc: float = 1e-08, max_iter: int = 50) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Find roots of Bessel first kind function j1(x) using Newton-Raphson iteration

First k approximations equation j1_roots ~ pi * (k + 1/4). See Baricz et al. (2025)

Parameters:

nb_rootsint: The number of j1(x)=0 to find
accfloat, optional: The tolerance for numerical errors. Default is 1e-8
max_iterint, optional: The maximum number of iterations trying to improve the error accuracy. Default is 50

Returns:

ndarray: The roots of j1(x), a 1-D ndarray of size nb_roots

Notes

The scipy equivalent -> scipy.special.jn_zeros(1,x)

References

Baricz, Á., Kumar, P., & Ponnusamy, S. (2025). Asymptotic behavior of zeros of Bessel function derivatives. arXiv preprint arXiv:2510.12353.

Examples

>>> from pytrunc.utils import bessel_j1_roots
>>> bessel_j1_roots(3)
array([ 3.83170597,  7.01558667, 10.17346814])

pytrunc.utils.bessel_jn(x: ArrayLike, n: int, acc: float = 1e-08, max_iter: int = 50) → ndarray[tuple[Any, ...], dtype[float64]][source]#

The Bessel first kind function Jn(x) of order n

Parameters:

xndarray: The variable x of Jn(x), it must be 1-D
nint: The Bessel first kind function order
accfloat, optional: The tolerance for numerical errors. Default is 1e-8
max_iterint, optional: The maximum number of iterations trying to improve the error accuracy. Default is 50

Returns:

ndarray: The values of the Bessel function Jn(x), a 1-D ndarray with the same shape as x

Notes

The scipy equivalent -> scipy.special.jn(n,x)

Examples

>>> import numpy as np
>>> from pytrunc.utils import bessel_jn
>>> x = np.array([0.0, 1.0, 2.0])
>>> bessel_jn(x, n=2)
array([0.        , 0.11490348, 0.35283403])

pytrunc.utils.integrate_lobatto(f: ndarray[tuple[Any, ...], dtype[float64]], x: ndarray[tuple[Any, ...], dtype[float64]], lp: int | None = None, xk: ndarray[tuple[Any, ...], dtype[float64]] | None = None, wk: ndarray[tuple[Any, ...], dtype[float64]] | None = None, assume_sorted: bool = False) → float[source]#

Integrate using lobatto quadrature

Parameters:

fndarray: The ordinates of the function (array to be integrated), it must be 1-D
xndarray: The abscissas, it must be 1-D
lpint or None, optional: The number of lobatto points for the integration. Default is None, meaning lp is equal to len(x)
xkndarray or None, optional: Force the Lobatto quadrature abscissas, it must be 1-D and of size lp. Considered only if wk is also provided
wkndarray or None, optional: Force the Lobatto weights, it must be 1-D and of size lp. Considered only if xk is also provided
assume_sortedbool, optional: If True, the x array is assumed to be sorted in ascending order. Default is False

Returns:

float: The estimated integral calculated using the Lobatto quadrature

Examples

>>> import numpy as np
>>> from pytrunc.utils import integrate_lobatto
>>> theta = np.linspace(0.0, np.pi, 1801)
>>> integrate_lobatto(np.sin(theta), theta, assume_sorted=True)
1.999999497717348

pytrunc.utils.legendre_polynomials(n: int, x: ArrayLike) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Use the recursion formulas to compute the Legendre polynomials Pn(x)

see Eq. 9 in Michels (1963)

Parameters:

nint: The Legendre polynomial order
xndarray: The x values of Pn(x), it must be 1-D

Returns:

ndarray: The Legendre series, a 1-D ndarray with the same shape as x

Notes

The numpy equivalent -> numpy.polynomial.legendre.Legendre.basis(n)(x)

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Wiscombe, W. J. (1977). The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. Journal of Atmospheric Sciences, 34(9), 1408-1422.

Examples

>>> import numpy as np
>>> from pytrunc.utils import legendre_polynomials
>>> x = np.array([-1.0, 0.0, 1.0])
>>> legendre_polynomials(2, x)
array([ 1. , -0.5,  1. ])

pytrunc.utils.legendre_polynomials_derivative(n: int, x: ArrayLike) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Use the recursion formulas to compute the derivative Legendre polynomials d(Pn(x))

see Eq. 10 in Michels (1963)

Parameters:

nint: The Legendre polynomial order
xndarray: The x values of d(Pn(x)), it must be 1-D

Returns:

ndarray: The derivative Legendre series, a 1-D ndarray with the same shape as x

Notes

The numpy equivalent -> numpy.polynomial.legendre.Legendre.basis(n).deriv(1)(x)

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Wiscombe, W. J. (1977). The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. Journal of Atmospheric Sciences, 34(9), 1408-1422.

Examples

>>> import numpy as np
>>> from pytrunc.utils import legendre_polynomials_derivative
>>> x = np.array([-1.0, 0.0, 1.0])
>>> legendre_polynomials_derivative(2, x)
array([-3.,  0.,  3.])

pytrunc.utils.legendre_polynomials_derivative_roots(n: int, acc: float = 1e-08, max_iter: int = 50) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Find roots of legendre polynomial derivative d(Pn(x)), for x > -1 and x < 1

Use of the Newton-Raphson iteration as in Michels (1963), see Eq. 7 and 8

Parameters:

nint: The Legendre polynomial order
accfloat, optional: The tolerance for numerical errors. Default is 1e-8
max_iterint, optional: The maximum number of iterations trying to improve the error accuracy. Default is 50

Returns:

ndarray: The roots of the legendre polynomial derivative d(Pn(x)), a 1-D ndarray of size n-1

Notes

The numpy equivalent -> Legendre.basis(n).deriv().roots()
Faster than the numpy equivalent!

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Examples

>>> from pytrunc.utils import (
...     legendre_polynomials_derivative_roots,
... )
>>> legendre_polynomials_derivative_roots(4)
array([-0.65465367,  0.        ,  0.65465367])

pytrunc.utils.legendre_polynomials_second_derivative(n: int, x: ArrayLike) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Use the recursion formulas to compute the second derivative Legendre polynomials d²(Pn(x))

see Eq. 11 in Michels (1963)

Parameters:

nint: The Legendre polynomial order
xndarray: The x values of d²(Pn(x)), it must be 1-D

Returns:

ndarray: The second derivative Legendre series, a 1-D ndarray with the same shape as x

Notes

The numpy equivalent -> numpy.polynomial.legendre.Legendre.basis(n).deriv(2)(x)

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Wiscombe, W. J. (1977). The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. Journal of Atmospheric Sciences, 34(9), 1408-1422.

Examples

>>> import numpy as np
>>> from pytrunc.utils import (
...     legendre_polynomials_second_derivative,
... )
>>> x = np.array([-1.0, 0.0, 1.0])
>>> legendre_polynomials_second_derivative(3, x)
array([-15.,   0.,  15.])

pytrunc.utils.quadrature_lobatto(abscissa_min: float = -1, abscissa_max: float = 1, n: int = 100) → tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]][source]#

Compute the abscissas (sample points) and weights for Lobatto quadrature

Parameters:

abscissa_minfloat, optional: The min fixed abscissa. Default is -1
abscissa_maxfloat, optional: The max fixed abscissa. Default is 1
nint, optional: The number of abscissas / weights. Default is 100

Returns:

abscissasndarray: The Lobatto quadrature abscissas, a 1-D ndarray of size n
weightsndarray: The Lobatto quadrature weights, a 1-D ndarray of size n

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Wiscombe, W. J. (1977). The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. Journal of Atmospheric Sciences, 34(9), 1408-1422.

Examples

>>> from pytrunc.utils import quadrature_lobatto
>>> xk, wk = quadrature_lobatto(n=4)
>>> xk
array([-1.       , -0.4472136,  0.4472136,  1.       ])
>>> wk
array([0.16666667, 0.83333333, 0.83333333, 0.16666667])

pytrunc.utils module#

This Page