mathematicskit.special_functions#
The gamma and beta functions (scipy.special.gamma/beta); Bessel
functions of the first/second kind (scipy.special.jv/yv);
orthogonal polynomial families (Legendre, Chebyshev, Hermite, Laguerre)
via numpy.polynomial; and the discrete Fourier transform via
numpy.fft, plus a from-scratch radix-2 FFT kept as a pedagogical
comparison against a naive DFT and numpy.fft. Signal processing and
transforms: direct and FFT convolution (numpy.convolve,
scipy.signal.fftconvolve); Butterworth, Chebyshev, and window-method
FIR filters (scipy.signal); the Z-transform with pole-zero analysis
and partial-fraction inversion; the Laplace transform with hand-rolled
Talbot and Gaver-Stehfest inversion; and hand-rolled Daubechies
wavelets, the discrete wavelet transform, and the Morlet continuous
wavelet transform.
mathematicskit.special_functions: special functions and transforms.
The gamma and beta functions (scipy.special.gamma/beta) with
Stirling’s asymptotic series; Bessel, Airy, and Mathieu functions;
elliptic integrals, the arithmetic-geometric mean, and Jacobi elliptic
functions; Gauss’s and Kummer’s hypergeometric functions; the error
function and Fresnel integrals; the Riemann zeta function and Euler’s
prime product; the Lambert W function; orthogonal polynomial families
(Legendre, Chebyshev, Hermite, Laguerre) via numpy.polynomial, with
orthogonality verified numerically in tests; and the discrete Fourier
transform via numpy.fft, plus a from-scratch radix-2 FFT kept
specifically as a pedagogical comparison against a naive \(O(n^2)\)
DFT and numpy.fft – the one sub-module in this domain where
hand-rolling is the point.
Signal processing and transforms: direct and FFT convolution and
correlation (numpy.convolve, scipy.signal.fftconvolve);
Butterworth, Chebyshev, and window-method FIR filter design and
frequency responses (scipy.signal); the Z-transform, pole-zero
analysis, and inversion by partial fractions (scipy.signal.residuez);
the Laplace transform by quadrature with hand-rolled Talbot and
Gaver-Stehfest numerical inversion; and hand-rolled Daubechies wavelets,
the discrete wavelet transform, and the Morlet continuous wavelet
transform (no numpy/scipy equivalent since scipy 1.15 removed its
wavelet routines).
- class mathematicskit.special_functions.AiryResult(x, ai, ai_prime, bi, bi_prime)[source]#
Bases:
objectThe Airy functions
Ai,Biand their derivatives at the sample pointsx.
- class mathematicskit.special_functions.ConvolutionComparisonResult(direct_result, fft_result, direct_time, fft_time, max_error)[source]#
Bases:
objectContainer comparing direct \(O(nm)\) convolution with FFT-based convolution.
- Parameters:
- class mathematicskit.special_functions.FFTComparisonResult(naive_result, radix2_result, numpy_result, naive_time, radix2_time, numpy_time, max_error_naive_vs_numpy, max_error_radix2_vs_numpy)[source]#
Bases:
objectContainer comparing naive DFT, hand-rolled radix-2 FFT, and
numpy.fft.- Parameters:
- class mathematicskit.special_functions.FilterCoefficients(b, a)[source]#
Bases:
objectA digital filter’s transfer function \(H(z) = B(z)/A(z)\).
Both polynomials are in powers of \(z^{-1}\), following
scipy.signal.lfilter(): \(B(z) = b_0 + b_1 z^{-1} + \cdots\). An FIR filter hasa == [1.0].
- class mathematicskit.special_functions.FrequencyResponseResult(frequencies, response, magnitude, phase)[source]#
Bases:
objectA filter’s frequency response \(H(e^{i\omega})\) sampled on \([0, \pi)\) (or \([0, f_s/2)\)).
- class mathematicskit.special_functions.FresnelResult(t, s, c)[source]#
Bases:
objectThe Fresnel integrals \(S(t)\) and \(C(t)\) at the sample points
t.
- class mathematicskit.special_functions.JacobiEllipticResult(sn, cn, dn, amplitude)[source]#
Bases:
objectJacobi elliptic functions \(\operatorname{sn}\), \(\operatorname{cn}\), \(\operatorname{dn}\) at
(u, m).
- class mathematicskit.special_functions.PoleZeroResult(zeros, poles, gain, is_stable)[source]#
Bases:
objectPoles, zeros, and gain of a rational transfer function in \(z\).
- class mathematicskit.special_functions.ScalogramResult(times, scales, frequencies, coefficients)[source]#
Bases:
objectA continuous wavelet transform sampled on a time-scale grid.
- class mathematicskit.special_functions.WaveletDecompositionResult(approximation, details, wavelet)[source]#
Bases:
objectA multilevel discrete wavelet decomposition.
- mathematicskit.special_functions.airy_functions(x)[source]#
The Airy functions \(\operatorname{Ai}, \operatorname{Ai}', \operatorname{Bi}, \operatorname{Bi}'\) at
x.\(\operatorname{Ai}(x) = \frac1\pi\int_0^\infty \cos(t^3/3 + xt)\,dt\) decays as \(x\to+\infty\); \(\operatorname{Bi}\) grows. Both satisfy \(y'' = xy\). Via
scipy.special.airy().Examples
>>> from scipy.special import gamma >>> r = airy_functions(0.0) >>> round(float(r.ai), 12) == round(1 / (3 ** (2 / 3) * float(gamma(2 / 3))), 12) True
- mathematicskit.special_functions.apply_filter(filt, x, zero_phase=False)[source]#
Run a signal through a filter by its difference equation.
Solves \(\sum_k a_k\,y[n-k] = \sum_k b_k\,x[n-k]\) with zero initial conditions (
scipy.signal.lfilter()), or filters forward and then backward for zero phase distortion (scipy.signal.filtfilt(), squaring the magnitude response).- Parameters:
filt (
FilterCoefficients)x (
ndarray)zero_phase (
bool)
- Return type:
- Returns:
ndarray
Examples
>>> import numpy as np >>> moving_average = FilterCoefficients(b=np.ones(3) / 3, a=np.array([1.0])) >>> apply_filter(moving_average, np.array([3.0, 3.0, 3.0, 3.0])) array([1., 2., 3., 3.])
- mathematicskit.special_functions.arithmetic_geometric_mean(a, b)[source]#
Gauss’s arithmetic-geometric mean \(\mathrm{AGM}(a, b)\).
Iterate \(a_{n+1} = (a_n + b_n)/2\), \(b_{n+1} = \sqrt{a_n b_n}\); both sequences converge quadratically to a common limit. Via
scipy.special.agm(). See C. F. Gauss, Werke, vol. 3 (1866), 361-403.- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> # Gauss's constant 1/AGM(1, sqrt 2) = 0.8346268... >>> round(1 / float(arithmetic_geometric_mean(1.0, math.sqrt(2.0))), 7) 0.8346268
- mathematicskit.special_functions.bessel_first_kind(nu, x)[source]#
Bessel function of the first kind, \(J_\nu(x)\).
Solves Bessel’s differential equation \(x^2y'' + xy' + (x^2 - \nu^2)y = 0\), regular (finite) at \(x=0\). Via
scipy.special.jv(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 9.1.Examples
>>> round(float(bessel_first_kind(0.0, 0.0)), 6) # J_0(0) = 1 1.0 >>> round(float(bessel_first_kind(1.0, 0.0)), 6) # J_nu(0) = 0 for nu > 0 0.0
- mathematicskit.special_functions.bessel_second_kind(nu, x)[source]#
Bessel function of the second kind (Weber/Neumann function), \(Y_\nu(x)\).
The second, linearly independent solution of Bessel’s equation, singular (diverging to \(-\infty\)) at \(x=0\). Via
scipy.special.yv(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 9.1.Examples
>>> bool(bessel_second_kind(0.0, 0.1) < bessel_second_kind(0.0, 1.0)) # diverges toward x=0 True
- mathematicskit.special_functions.beta_function(a, b)[source]#
The beta function \(B(a,b) = \dfrac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt\).
Via
scipy.special.beta(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 6.2.- Parameters:
- Returns:
float or ndarray
Examples
>>> from scipy import special >>> round(float(beta_function(2.0, 3.0)), 6) == round(float(special.gamma(2.0) * special.gamma(3.0) / special.gamma(5.0)), 6) True
- mathematicskit.special_functions.butterworth_filter(order, cutoff, btype='lowpass', fs=None)[source]#
Design a digital Butterworth filter via
scipy.signal.butter().- Parameters:
- Return type:
- Returns:
FilterCoefficients
Examples
>>> import numpy as np >>> filt = butterworth_filter(4, 0.25) >>> response = frequency_response(filt, worN=[0.25 * np.pi]) >>> round(float(response.magnitude[0] ** 2), 6) # half power at the cutoff 0.5
- mathematicskit.special_functions.chebyshev1_filter(order, ripple_db, cutoff, btype='lowpass', fs=None)[source]#
Design a digital Chebyshev type I filter via
scipy.signal.cheby1().The passband gain ripples between \(1\) and \(10^{-r/20}\) (
r = ripple_db);cutoffis the frequency where the gain last reaches \(10^{-r/20}\).- Parameters:
- Return type:
- Returns:
FilterCoefficients
Examples
>>> import numpy as np >>> filt = chebyshev1_filter(5, 1.0, 0.3) >>> gain_db = 20 * np.log10(frequency_response(filt, worN=[0.3 * np.pi]).magnitude[0]) >>> round(float(gain_db), 6) -1.0
- mathematicskit.special_functions.chebyshev_polynomial(n, x)[source]#
The Chebyshev polynomial (first kind) \(T_n(x) = \cos(n\arccos x)\), orthogonal on \([-1,1]\) with weight \(1/\sqrt{1-x^2}\).
Via
numpy.polynomial.chebyshev.chebval(). The minimal-sup-norm property of (monic-scaled) Chebyshev polynomials is exactly why Chebyshev nodes (seemathematicskit.numerical_analysis.systems.chebyshev) avoid the Runge phenomenon. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.4.- Parameters:
n (
int)x (float or array-like of float, in
[-1, 1])
- Returns:
float or ndarray
Examples
>>> round(float(chebyshev_polynomial(3, 0.5)), 6) # T_3(cos(pi/3)) = cos(pi) = -1 -1.0
- mathematicskit.special_functions.circular_convolve(x, h)[source]#
Circular (periodic) convolution of two equal-length sequences.
\((x \circledast h)[n] = \sum_{k=0}^{N-1} x[k]\,h[(n-k) \bmod N] = \mathcal{F}^{-1}\{X_k H_k\}\), the convolution theorem for the DFT applied directly with
numpy.fft.fft(). Real inputs give a real output.- Parameters:
- Return type:
- Returns:
ndarray
Examples
>>> import numpy as np >>> np.allclose(circular_convolve([1.0, 2.0, 3.0, 4.0], [0.0, 1.0, 0.0, 0.0]), [4.0, 1.0, 2.0, 3.0]) True
- mathematicskit.special_functions.compare_convolution_methods(n, m, seed=0)[source]#
Time and cross-check direct and FFT convolution of random signals of lengths
nandm.Direct convolution scales as \(O(nm)\); FFT convolution as \(O((n+m)\log(n+m))\), so the FFT wins once the kernel is more than a few dozen samples long.
- Parameters:
- Return type:
- Returns:
ConvolutionComparisonResult
Examples
>>> result = compare_convolution_methods(1000, 200, seed=0) >>> result.max_error < 1e-9 True
- mathematicskit.special_functions.compare_fft_methods(n, seed=0)[source]#
Time and cross-check the naive DFT, radix-2 FFT, and
numpy.ffton the same random signal.Demonstrates the \(O(n^2)\) vs. \(O(n\log n)\) gap directly: the naive DFT’s time grows quadratically with n while the radix-2 FFT’s (and
numpy.fft’s) grows only log-linearly, becoming dramatically faster for even moderately large n.- Parameters:
n (
int) – Signal length; must be a power of 2 (forfft_radix2()).seed (
int)
- Return type:
- Returns:
FFTComparisonResult
Examples
>>> result = compare_fft_methods(256, seed=0) >>> result.max_error_naive_vs_numpy < 1e-8 True >>> result.max_error_radix2_vs_numpy < 1e-8 True
- mathematicskit.special_functions.complementary_error_function(x)[source]#
The complementary error function \(\operatorname{erfc}(x) = 1 - \operatorname{erf}(x)\).
Via
scipy.special.erfc(), which stays accurate in the far tail where computing1 - erf(x)directly would cancel to zero.Examples
>>> float(complementary_error_function(10.0)) > 0.0 # 1 - erf(10) is exactly 0.0 in float64 True
- mathematicskit.special_functions.complete_elliptic_integral_first_kind(m)[source]#
The complete elliptic integral of the first kind, \(K(m) = \int_0^{\pi/2} \frac{d\theta}{\sqrt{1 - m\sin^2\theta}}\).
Via
scipy.special.ellipk()(parameter \(m = k^2\)). See Abramowitz & Stegun, Sec. 17.3.Examples
>>> import math >>> abs(float(complete_elliptic_integral_first_kind(0.0)) - math.pi / 2) < 1e-15 True
- mathematicskit.special_functions.complete_elliptic_integral_second_kind(m)[source]#
The complete elliptic integral of the second kind, \(E(m) = \int_0^{\pi/2} \sqrt{1 - m\sin^2\theta}\,d\theta\).
\(4aE(e^2)\) is the perimeter of an ellipse with semi-major axis \(a\) and eccentricity \(e\), the problem that gave elliptic integrals their name. Via
scipy.special.ellipe(). See Abramowitz & Stegun, Sec. 17.3.Examples
>>> float(complete_elliptic_integral_second_kind(1.0)) # degenerate ellipse: a segment 1.0
- mathematicskit.special_functions.confluent_hypergeometric_1f1(a, b, z)[source]#
Kummer’s confluent hypergeometric function \(M(a, b, z) = {}_1F_1(a; b; z) = \sum_{n=0}^\infty \frac{(a)_n}{(b)_n}\frac{z^n}{n!}\).
It solves Kummer’s equation \(zw'' + (b - z)w' - aw = 0\) and arises as the confluent limit \({}_1F_1(a; b; z) = \lim_{c\to\infty} {}_2F_1(a, c; b; z/c)\). Kummer’s transformation reads \(M(a, b, z) = e^z M(b - a, b, -z)\). Via
scipy.special.hyp1f1(). See E. E. Kummer, “De integralibus quibusdam definitis et seriebus infinitis,” Journal für die reine und angewandte Mathematik 17 (1837), 228-242.Examples
>>> import math >>> round(float(confluent_hypergeometric_1f1(2.0, 2.0, 1.0)), 12) == round(math.e, 12) # 1F1(a; a; z) = e^z True
- mathematicskit.special_functions.convolve_direct(x, h, mode='full')[source]#
Linear convolution by the direct \(O(NM)\) sum, via
numpy.convolve().- Parameters:
- Return type:
- Returns:
ndarray
Examples
>>> convolve_direct([1.0, 2.0, 3.0], [0.0, 1.0, 0.5]) array([0. , 1. , 2.5, 4. , 1.5])
- mathematicskit.special_functions.convolve_fft(x, h, mode='full')[source]#
Linear convolution via the convolution theorem, using
scipy.signal.fftconvolve().Both inputs are zero-padded to length \(N + M - 1\) so the circular convolution the DFT computes equals the linear one.
- Parameters:
- Return type:
- Returns:
ndarray
Examples
>>> import numpy as np >>> np.allclose(convolve_fft([1.0, 2.0, 3.0], [0.0, 1.0, 0.5]), [0.0, 1.0, 2.5, 4.0, 1.5]) True
- mathematicskit.special_functions.cross_correlate(x, y, mode='full')[source]#
Cross-correlation \((x \star y)[k] = \sum_n x[n+k]\,\overline{y[n]}\), via
scipy.signal.correlate().Correlation is convolution with the conjugated, time-reversed second signal, so its peak locates the lag at which
ybest matchesx. Withmode="full", output indexicorresponds to lagi - (len(y) - 1)(seescipy.signal.correlation_lags()).Examples
>>> import numpy as np >>> c = cross_correlate([0.0, 0.0, 1.0, 2.0, 0.0], [1.0, 2.0]) >>> int(np.argmax(c)) - 1 # lag of the best match 2
- mathematicskit.special_functions.daubechies_filter(p)[source]#
The orthonormal Daubechies lowpass (scaling) filter with
pvanishing moments.The filter has length \(2p\) and frequency response \(H(\omega) = \sqrt2\bigl(\tfrac{1+e^{-i\omega}}{2}\bigr)^p Q(e^{-i\omega})\), where \(|Q|^2 = P(\sin^2(\omega/2))\), \(P(y) = \sum_{k=0}^{p-1}\binom{p-1+k}{k}y^k\). Substituting \(y = (2 - z - z^{-1})/4\) turns \(P\) into a polynomial whose roots come in reciprocal pairs \((r, 1/r)\); keeping the roots inside the unit circle gives the minimum-phase \(Q\). Root-finding becomes ill-conditioned beyond \(p \approx 20\).
- Parameters:
p (
int) – Number of vanishing moments ("db<p>"),p >= 1.- Return type:
- Returns:
ndarray – Coefficients \(h[0], \ldots, h[2p-1]\), with \(\sum h = \sqrt2\) and \(\sum h^2 = 1\).
Examples
>>> import numpy as np >>> h = daubechies_filter(2) >>> r3 = np.sqrt(3.0) >>> bool(np.allclose(h, np.array([1 + r3, 3 + r3, 3 - r3, 1 - r3]) / (4 * np.sqrt(2.0)))) True
- mathematicskit.special_functions.dft_naive(x)[source]#
The discrete Fourier transform, by direct evaluation of its defining sum.
\(X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i kn/N}\), computed as a dense matrix-vector product – \(O(N^2)\) time. Kept hand-rolled purely to demonstrate the speedup
fft_radix2()achieves over it. See Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 30.1.- Parameters:
x (
ndarray) – Complex (or real) input signal, any length.- Return type:
- Returns:
ndarray, complex
Examples
>>> import numpy as np >>> x = np.array([1.0, 2.0, 3.0, 4.0]) >>> np.allclose(dft_naive(x), np.fft.fft(x)) True
- mathematicskit.special_functions.discrete_wavelet_transform(x, wavelet='haar', level=None)[source]#
Multilevel orthogonal discrete wavelet transform with periodic boundary handling.
Each level computes \(a[k] = \sum_j h[j]\,x[(2k+j) \bmod N]\) and \(d[k] = \sum_j g[j]\,x[(2k+j) \bmod N]\), then recurses on \(a\) (Mallat 1989). The transform is orthogonal: \(\|x\|^2 = \|a\|^2 + \sum \|d_\ell\|^2\).
- Parameters:
x (
ndarray) – 1-D real signal;len(x)must be divisible by2**level.wavelet (
str) – Seewavelet_filters().level (
int|None) – Number of levels; defaults to the largest for which every level’s input is at least as long as the filter.
- Return type:
- Returns:
WaveletDecompositionResult
Examples
>>> import numpy as np >>> result = discrete_wavelet_transform(np.array([4.0, 4.0, 2.0, 0.0]), "haar", level=1) >>> np.round(result.approximation * np.sqrt(2), 12), np.round(result.details[0] * np.sqrt(2), 12) (array([8., 2.]), array([0., 2.]))
- mathematicskit.special_functions.error_function(x)[source]#
The error function \(\operatorname{erf}(x) = \frac{2}{\sqrt\pi}\int_0^x e^{-t^2}\,dt\).
The probability that a normal variable lies within \(x\sqrt2\) standard deviations of its mean. Named by J. W. L. Glaisher, “On a class of definite integrals,” Philosophical Magazine 42 (1871), 294-302. Via
scipy.special.erf().Examples
>>> import math >>> round(float(error_function(1 / math.sqrt(2))), 4) # the "68%" of the 68-95-99.7 rule 0.6827
- mathematicskit.special_functions.euler_product(s, prime_bound)[source]#
Euler’s product \(\prod_{p \le P} (1 - p^{-s})^{-1}\) over primes up to
prime_bound.Euler (1737) showed that for \(s > 1\) this product converges to \(\zeta(s)\) as \(P\to\infty\), by unique prime factorization. See L. Euler, “Variae observationes circa series infinitas,” Commentarii Academiae Scientiarum Petropolitanae 9 (1744), 160-188 (presented 1737).
- Parameters:
- Returns:
float
Examples
>>> round(euler_product(2.0, 10), 6) # primes 2, 3, 5, 7 1.595052
- mathematicskit.special_functions.fft_numpy(x)[source]#
The discrete Fourier transform via
numpy.fft.fft().The primary API for actually computing a DFT in this domain: a highly optimized mixed-radix FFT (handling any length, not just powers of 2). See NumPy’s FFT documentation.
Examples
>>> import numpy as np >>> np.allclose(fft_numpy(np.array([1.0, 0.0, -1.0, 0.0])), np.fft.fft([1.0, 0.0, -1.0, 0.0])) True
- mathematicskit.special_functions.fft_radix2(x)[source]#
The Cooley-Tukey radix-2 fast Fourier transform.
Recursively splits the DFT of length \(N\) (a power of 2) into two DFTs of length \(N/2\) (even- and odd-indexed samples), combined via the “butterfly” identity \(X_k = E_k + e^{-2\pi ik/N}O_k\), \(X_{k+N/2} = E_k - e^{-2\pi ik/N}O_k\) – giving \(O(N\log N)\) time overall. Requires
len(x)to be a power of 2. See Cooley & Tukey (1965), Math. Comp. 19, and Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 30.2.- Parameters:
x (
ndarray) – Complex (or real) input signal;len(x)must be a power of 2.- Return type:
- Returns:
ndarray, complex
Examples
>>> import numpy as np >>> x = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]) >>> np.allclose(fft_radix2(x), np.fft.fft(x)) True
- mathematicskit.special_functions.fir_window_filter(numtaps, cutoff, window='hamming', pass_zero=True, fs=None)[source]#
Design a linear-phase FIR filter by the window method, via
scipy.signal.firwin().The ideal brick-wall impulse response, a sampled sinc, is truncated to
numtapssamples and multiplied bywindow. The result is symmetric, so its phase is exactly linear (a pure delay of(numtaps - 1) / 2samples).- Parameters:
numtaps (
int) – Filter length (odd for a highpass).window (
str|tuple) – Anyscipy.signal.get_window()spec, e.g."hamming"or("kaiser", 8.0).pass_zero (
bool|str) –Truefor lowpass-style (passes DC),Falsefor highpass-style, or a band type string as inscipy.signal.firwin().
- Return type:
- Returns:
FilterCoefficients – With
a == [1.0].
Examples
>>> import numpy as np >>> filt = fir_window_filter(31, 0.4) >>> bool(np.allclose(filt.b, filt.b[::-1])) # symmetric => linear phase True
- mathematicskit.special_functions.frequency_response(filt, worN=512, fs=None)[source]#
Evaluate \(H(e^{i\omega})\) on the unit circle via
scipy.signal.freqz().- Parameters:
filt (
FilterCoefficients)worN (int or array_like) – Number of equally spaced frequencies on \([0, \pi)\), or the frequencies themselves (rad/sample, or units of
fs).
- Return type:
- Returns:
FrequencyResponseResult
Examples
>>> import numpy as np >>> response = frequency_response(butterworth_filter(2, 0.5), worN=[0.0]) >>> round(float(response.magnitude[0]), 12) # unit DC gain 1.0
- mathematicskit.special_functions.fresnel_integrals(t)[source]#
The Fresnel integrals \(S(t) = \int_0^t \sin(\pi\tau^2/2)\,d\tau\) and \(C(t) = \int_0^t \cos(\pi\tau^2/2)\,d\tau\).
The parametric curve \((C(t), S(t))\) is the Euler (Cornu) spiral, which winds into the points \(\pm(\tfrac12, \tfrac12)\) as \(t\to\pm\infty\). Via
scipy.special.fresnel(). See A. Fresnel, “Mémoire sur la diffraction de la lumière,” Mémoires de l’Académie Royale des Sciences de l’Institut de France 5 (1826), 339-475.Examples
>>> r = fresnel_integrals(1.0) >>> round(float(r.c), 6), round(float(r.s), 6) (0.779893, 0.438259)
- mathematicskit.special_functions.gamma_function(x)[source]#
The gamma function \(\Gamma(x) = \int_0^\infty t^{x-1}e^{-t}\,dt\).
Extends the factorial to real (and complex) arguments: \(\Gamma(n) = (n-1)!\) for positive integers
n. Viascipy.special.gamma(). See Abramowitz & Stegun, Handbook of Mathematical Functions, Sec. 6.1.Examples
>>> round(float(gamma_function(5.0)), 6) # 4! = 24 24.0 >>> round(float(gamma_function(0.5)), 6) # sqrt(pi) 1.772454
- mathematicskit.special_functions.hermite_polynomial(n, x)[source]#
The (physicists’) Hermite polynomial \(H_n(x)\), orthogonal on \((-\infty,\infty)\) with weight \(e^{-x^2}\).
Via
numpy.polynomial.hermite.hermval(). Arise as the eigenfunctions of the quantum harmonic oscillator. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.3.Examples
>>> float(hermite_polynomial(0, 2.0)) 1.0 >>> float(hermite_polynomial(1, 2.0)) # H_1(x) = 2x 4.0 >>> round(float(hermite_polynomial(2, 1.0)), 6) # H_2(x) = 4x^2 - 2 2.0
- mathematicskit.special_functions.hypergeometric_2f1(a, b, c, z)[source]#
Gauss’s hypergeometric function \({}_2F_1(a, b; c; z) = \sum_{n=0}^\infty \frac{(a)_n (b)_n}{(c)_n}\frac{z^n}{n!}\).
\((q)_n = q(q+1)\cdots(q+n-1)\) is the rising factorial. The series converges for \(|z| < 1\);
scipy.special.hyp2f1()uses analytic continuation elsewhere. Gauss’s summation theorem (1812) evaluates it at \(z = 1\) for \(\operatorname{Re}(c - a - b) > 0\):\[{}_2F_1(a, b; c; 1) = \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)}.\]See C. F. Gauss, “Disquisitiones generales circa seriem infinitam,” Commentationes Societatis Regiae Scientiarum Gottingensis Recentiores 2 (1813).
- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> # ln(1 + z) = z * 2F1(1, 1; 2; -z) >>> round(0.5 * float(hypergeometric_2f1(1.0, 1.0, 2.0, -0.5)), 12) == round(math.log(1.5), 12) True
- mathematicskit.special_functions.inner_product(f, g, weight, a, b)[source]#
The weighted inner product \(\langle f, g\rangle = \int_a^b f(x)g(x)w(x)\,dx\), via
scipy.integrate.quad().Used to numerically verify an orthogonal polynomial family’s defining property directly: \(\langle p_m, p_n\rangle_w = 0\) for \(m \neq n\). See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.1.
- Parameters:
- Return type:
- Returns:
float
Examples
>>> from mathematicskit.special_functions.systems.orthogonal_polynomials import legendre_polynomial >>> p2 = lambda x: legendre_polynomial(2, x) >>> p3 = lambda x: legendre_polynomial(3, x) >>> abs(inner_product(p2, p3, lambda x: 1.0, -1.0, 1.0)) < 1e-10 True
- mathematicskit.special_functions.inverse_discrete_wavelet_transform(result)[source]#
Reconstruct a signal from its
discrete_wavelet_transform().Applies the transpose of each (orthogonal) analysis step, \(x[(2k+j) \bmod N] \mathrel{+}= a[k]\,h[j] + d[k]\,g[j]\), from the coarsest level to the finest. Zeroing small detail coefficients before reconstructing is wavelet compression/denoising.
- Parameters:
result (
WaveletDecompositionResult)- Return type:
- Returns:
ndarray
Examples
>>> import numpy as np >>> x = np.arange(8.0) >>> bool(np.allclose(inverse_discrete_wavelet_transform(discrete_wavelet_transform(x, "db2")), x)) True
- mathematicskit.special_functions.inverse_laplace_stehfest(F, t, n=14)[source]#
Invert a Laplace transform numerically by the Gaver-Stehfest method.
\[f(t) \approx \frac{\ln 2}{t}\sum_{k=1}^{n} V_k\,F\!\left(\frac{k\ln 2}{t}\right)\]with
stehfest_coefficients()\(V_k\). Only real \(s\) are needed. In double precision the error is smallest near \(n = 12\)-\(16\) (about 5-7 digits for smooth \(f\)); it cannot resolve oscillating or discontinuous \(f\).- Parameters:
- Return type:
- Returns:
ndarray – Same shape as
t.
Examples
>>> import numpy as np >>> f = inverse_laplace_stehfest(lambda s: 1.0 / (s + 1.0), [0.5, 1.0, 2.0]) >>> bool(np.allclose(f, np.exp(-np.array([0.5, 1.0, 2.0])), atol=1e-5)) True
- mathematicskit.special_functions.inverse_laplace_talbot(F, t, m=32)[source]#
Invert a Laplace transform numerically by the fixed Talbot method.
With \(r = 2m/(5t)\) and \(\theta_k = k\pi/m\), the contour \(s(\theta) = r\theta(\cot\theta + i)\) gives (Abate & Valkó 2004, Eq. 36)
\[f(t) \approx \frac{r}{m}\left[\tfrac12 F(r)e^{rt} + \sum_{k=1}^{m-1}\operatorname{Re}\Bigl(e^{t s(\theta_k)}F\bigl(s(\theta_k)\bigr) \bigl(1 + i\sigma(\theta_k)\bigr)\Bigr)\right], \quad \sigma(\theta) = \theta + (\theta\cot\theta - 1)\cot\theta .\]In double precision, accuracy saturates near \(10^{-10}\) for \(m \approx 20\)-\(40\); larger \(m\) amplifies roundoff. Every singularity of \(F\) must lie inside the contour, which crosses the imaginary axis at \(\pm i m\pi/(5t)\). Poles in the right half-plane are not supported, and a pole \(\alpha \pm i\beta\) near the imaginary axis (an oscillating \(f\)) needs \(t < m\pi/(5\beta)\), in practice with some margin: raise
mto reach later times.- Parameters:
- Return type:
- Returns:
ndarray – Same shape as
t.
Examples
>>> import numpy as np >>> f = inverse_laplace_talbot(lambda s: 1.0 / (s + 1.0), [0.5, 1.0, 2.0]) >>> bool(np.allclose(f, np.exp(-np.array([0.5, 1.0, 2.0])), atol=1e-9)) True
- mathematicskit.special_functions.inverse_z_transform(b, a, n)[source]#
The first
nsamples of the causal sequence whose Z-transform is \(B(z)/A(z)\).Expands \(H(z)\) in partial fractions with
scipy.signal.residuez(),\[H(z) = \sum_i \frac{r_i}{(1 - p_i z^{-1})^{m_i}} + \sum_j k_j z^{-j},\]then inverts term by term using the table pair \(1/(1 - p z^{-1})^m \leftrightarrow \binom{n+m-1}{m-1} p^n\) (Oppenheim & Schafer, Sec. 3.4). Repeated poles (within
residuez’s default tolerance,1e-3) take increasing powers \(m_i = 1, 2, \ldots\) in order.- Parameters:
- Return type:
- Returns:
ndarray – Real when
bandaare real.
Examples
>>> import numpy as np >>> np.round(inverse_z_transform([1.0], [1.0, -0.5], 4), 12) # 0.5**n array([1. , 0.5 , 0.25 , 0.125])
- mathematicskit.special_functions.jacobi_elliptic_functions(u, m)[source]#
The Jacobi elliptic functions \(\operatorname{sn}\), \(\operatorname{cn}\), \(\operatorname{dn}\).
Defined by inverting the incomplete elliptic integral of the first kind: if \(u = \int_0^\varphi d\theta/\sqrt{1 - m\sin^2\theta}\) then \(\operatorname{sn} u = \sin\varphi\), \(\operatorname{cn} u = \cos\varphi\), \(\operatorname{dn} u = \sqrt{1 - m\sin^2\varphi}\). They are doubly periodic; on the real axis, \(\operatorname{sn}\) and \(\operatorname{cn}\) have period \(4K(m)\). Via
scipy.special.ellipj(). See C. G. J. Jacobi, Fundamenta nova theoriae functionum ellipticarum (1829).- Parameters:
- Returns:
JacobiEllipticResult
Examples
>>> r = jacobi_elliptic_functions(0.7, 0.5) >>> round(float(r.sn**2 + r.cn**2), 12) 1.0
- mathematicskit.special_functions.laguerre_polynomial(n, x)[source]#
The (simple) Laguerre polynomial \(L_n(x)\), orthogonal on \([0,\infty)\) with weight \(e^{-x}\).
Via
numpy.polynomial.laguerre.lagval(). Arise as the radial part of the hydrogen atom’s wavefunctions. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.5.- Parameters:
n (
int)x (float or array-like of float,
x >= 0)
- Returns:
float or ndarray
Examples
>>> float(laguerre_polynomial(0, 3.0)) 1.0 >>> round(float(laguerre_polynomial(1, 1.0)), 6) # L_1(x) = 1 - x 0.0
- mathematicskit.special_functions.lambert_w(z, branch=0)[source]#
The Lambert W function: the solution \(w\) of \(we^w = z\).
Via
scipy.special.lambertw(). For real \(z \ge -1/e\) the principal branch (branch=0) is real, and for \(-1/e \le z < 0\) the lower branchbranch=-1is a second real solution; in those cases the imaginary part (identically zero) is dropped and a real result is returned. See Corless et al. (1996).- Parameters:
- Returns:
float, complex, or ndarray
Examples
>>> import math >>> round(float(lambert_w(math.e)), 12) # 1 * e^1 = e 1.0 >>> # both real branches at z = -ln(2)/2: W_0 = -ln 2, W_{-1} = -ln 4 >>> round(float(lambert_w(-math.log(2) / 2)), 12) == round(-math.log(2), 12) True >>> round(float(lambert_w(-math.log(2) / 2, branch=-1)), 12) == round(-math.log(4), 12) True
- mathematicskit.special_functions.laplace_transform(f, s)[source]#
The Laplace transform \(F(s) = \int_0^\infty f(t)\,e^{-st}\,dt\), by adaptive quadrature.
Uses
scipy.integrate.quad()on \([0, \infty)\), integrating the real and imaginary parts separately for complex \(s\). \(\operatorname{Re} s\) must exceed the abscissa of convergence of \(f\).- Parameters:
- Return type:
- Returns:
ndarray – Same shape as
s; real when everysis real.
Examples
>>> import numpy as np >>> F = laplace_transform(np.sin, [1.0, 2.0]) # 1 / (s^2 + 1) >>> np.allclose(F, [0.5, 0.2]) True
- mathematicskit.special_functions.legendre_polynomial(n, x)[source]#
The Legendre polynomial \(P_n(x)\), orthogonal on \([-1,1]\) with weight 1.
Via
numpy.polynomial.legendre.legval()on the basis coefficients for degreen. Solutions of Legendre’s differential equation; arise e.g. in the multipole expansion of the electrostatic/gravitational potential. See Arfken, Weber & Harris, Mathematical Methods for Physicists, 7th ed., Ch. 18.2.- Parameters:
n (
int)x (float or array-like of float, in
[-1, 1])
- Returns:
float or ndarray
Examples
>>> import numpy as np >>> round(float(legendre_polynomial(2, 1.0)), 6) # P_n(1) = 1 for all n 1.0 >>> np.round(legendre_polynomial(2, np.array([-1.0, 0.0, 1.0])), 6) array([ 1. , -0.5, 1. ])
- mathematicskit.special_functions.log_gamma_function(x)[source]#
The (natural) log of the gamma function, \(\ln|\Gamma(x)|\).
Via
scipy.special.gammaln()– numerically stable for largex, where \(\Gamma(x)\) itself overflowsfloat64long before its logarithm does. See NIST Digital Library of Mathematical Functions, Sec. 5.4.Examples
>>> round(float(log_gamma_function(171.0)), 4) 706.5731
- mathematicskit.special_functions.mathieu_characteristic_a(m, q)[source]#
Characteristic value \(a_m(q)\) for the even Mathieu function \(\operatorname{ce}_m\).
At \(q = 0\) the equation is \(y'' + ay = 0\) and \(a_m(0) = m^2\). Via
scipy.special.mathieu_a().Examples
>>> float(mathieu_characteristic_a(3, 0.0)) 9.0
- mathematicskit.special_functions.mathieu_characteristic_b(m, q)[source]#
Characteristic value \(b_m(q)\) for the odd Mathieu function \(\operatorname{se}_m\).
\(b_m(0) = m^2\). Via
scipy.special.mathieu_b().Examples
>>> float(mathieu_characteristic_b(2, 0.0)) 4.0
- mathematicskit.special_functions.mathieu_even(m, q, x)[source]#
The even Mathieu function \(\operatorname{ce}_m(x, q)\), with
xin radians.Normalized so that \(\operatorname{ce}_m(x, 0) = \cos mx\) for \(m \ge 1\) (and \(1/\sqrt2\) for \(m = 0\)). Via
scipy.special.mathieu_cem().- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(mathieu_even(2, 0.0, 0.3)), 12) == round(math.cos(0.6), 12) True
- mathematicskit.special_functions.mathieu_odd(m, q, x)[source]#
The odd Mathieu function \(\operatorname{se}_m(x, q)\), with
xin radians.Normalized so that \(\operatorname{se}_m(x, 0) = \sin mx\). Via
scipy.special.mathieu_sem().- Parameters:
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(mathieu_odd(1, 0.0, 0.3)), 12) == round(math.sin(0.3), 12) True
- mathematicskit.special_functions.morlet_cwt(x, scales, dt=1.0, omega0=6.0)[source]#
Continuous wavelet transform with the Morlet wavelet, computed by FFT.
The Morlet wavelet \(\psi(\eta) = \pi^{-1/4}e^{i\omega_0\eta}e^{-\eta^2/2}\) has Fourier transform \(\hat\psi(\omega) = \pi^{-1/4}e^{-(\omega-\omega_0)^2/2}\) for \(\omega > 0\). At each scale \(s\),
\[W(s, t_n) = \mathcal{F}^{-1}\Bigl\{\hat x(\omega_k)\, \sqrt{2\pi s/\delta t}\;\hat\psi(s\omega_k)\Bigr\}[n]\](Torrence & Compo 1998, Eqs. 4, 6), a periodic convolution. Scale \(s\) corresponds to Fourier period \(4\pi s/(\omega_0 + \sqrt{2 + \omega_0^2})\), so a sinusoid of frequency \(f\) peaks near \(s \approx 1.03/f\) for \(\omega_0 = 6\).
- Parameters:
- Return type:
- Returns:
ScalogramResult
Examples
>>> import numpy as np >>> t = np.arange(1024) * 0.01 >>> scales = np.geomspace(0.05, 1.0, 60) >>> result = morlet_cwt(np.sin(2 * np.pi * 5.0 * t), scales, dt=0.01) >>> peak = np.argmax(np.mean(np.abs(result.coefficients) ** 2, axis=1)) >>> bool(abs(result.frequencies[peak] - 5.0) < 0.3) True
- mathematicskit.special_functions.poles_zeros(b, a)[source]#
Poles, zeros, and gain of \(H(z) = B(z)/A(z)\), and whether the causal system is stable.
bandaare zero-padded to equal length before callingscipy.signal.tf2zpk(), so that zeros or poles at \(z = 0\) implied by the \(z^{-1}\) convention are kept.- Parameters:
- Return type:
- Returns:
PoleZeroResult
Examples
>>> result = poles_zeros([1.0], [1.0, -0.9]) >>> result.poles, result.zeros, result.is_stable (array([0.9]), array([0.]), True)
- mathematicskit.special_functions.riemann_zeta(s)[source]#
The Riemann zeta function \(\zeta(s) = \sum_{n=1}^\infty n^{-s}\) (for \(s > 1\)), analytically continued.
Via
scipy.special.zeta(), real arguments only. Euler’s 1735 solution of the Basel problem gives \(\zeta(2) = \pi^2/6\); the continuation gives values such as \(\zeta(0) = -1/2\) and \(\zeta(-1) = -1/12\), and the “trivial zeros” \(\zeta(-2n) = 0\). See NIST Digital Library of Mathematical Functions, Ch. 25.- Parameters:
s (float or array-like of float) – Real argument, \(s \ne 1\).
- Returns:
float or ndarray
Examples
>>> import math >>> round(float(riemann_zeta(2.0)), 12) == round(math.pi**2 / 6, 12) True >>> round(float(riemann_zeta(-1.0)), 12) -0.083333333333
- mathematicskit.special_functions.stehfest_coefficients(n)[source]#
The Gaver-Stehfest weights \(V_1, \ldots, V_n\) (
neven).\[V_k = (-1)^{k + n/2} \sum_{j=\lfloor (k+1)/2 \rfloor}^{\min(k, n/2)} \frac{j^{n/2}\,(2j)!}{(n/2 - j)!\,j!\,(j-1)!\,(k-j)!\,(2j-k)!}\](Stehfest 1970). Computed in exact integer arithmetic, then converted to float; the weights sum to zero and grow like \(10^{n/2}\), which is what limits the method’s accuracy in double precision.
Examples
>>> stehfest_coefficients(4) array([ -2., 26., -48., 24.])
- mathematicskit.special_functions.stirling_factorial(n)[source]#
Stirling’s approximation \(n! \approx \sqrt{2\pi n}\,(n/e)^n\).
The leading term of Stirling’s series (J. Stirling, Methodus Differentialis, 1730, Prop. 28; A. de Moivre, Miscellanea Analytica, 1730). The relative error is about \(1/(12n)\), so the ratio to \(n!\) tends to 1 even though the absolute error grows without bound.
Examples
>>> import math >>> round(float(stirling_factorial(10.0)) / math.factorial(10), 4) 0.9917
- mathematicskit.special_functions.stirling_log_gamma(x, terms=3)[source]#
Stirling’s asymptotic series for \(\ln\Gamma(x)\), truncated after
termscorrections.\[\ln\Gamma(x) \sim \left(x - \tfrac12\right)\ln x - x + \tfrac12\ln(2\pi) + \sum_{k=1}^{K} \frac{B_{2k}}{2k(2k-1)\,x^{2k-1}},\]where \(B_{2k}\) are the Bernoulli numbers (from
scipy.special.bernoulli()). The series is asymptotic, not convergent: for fixed \(x\), adding terms helps only up to roughly \(K \approx \pi x\), after which the error grows. See NIST Digital Library of Mathematical Functions, Eq. 5.11.1.- Parameters:
- Returns:
float or ndarray
Examples
>>> abs(float(stirling_log_gamma(10.0, terms=3)) - float(log_gamma_function(10.0))) < 1e-9 True
- mathematicskit.special_functions.transfer_function(b, a, z)[source]#
Evaluate the rational transfer function \(H(z) = B(z)/A(z)\) in powers of \(z^{-1}\).
- Parameters:
b (
ndarray) – Numerator and denominator coefficients, as inscipy.signal.lfilter().a (
ndarray) – Numerator and denominator coefficients, as inscipy.signal.lfilter().z (complex or ndarray)
- Return type:
- Returns:
complex or ndarray
Examples
>>> complex(transfer_function([1.0], [1.0, -0.5], 2.0)) # 1 / (1 - 0.5/2) (1.3333333333333333+0j)
- mathematicskit.special_functions.wavelet_filters(wavelet)[source]#
The orthonormal lowpass/highpass analysis filter pair for a named wavelet.
- Parameters:
wavelet (
str) –"haar"(same as"db1") or"db<p>"for a Daubechies wavelet withpvanishing moments.- Return type:
- Returns:
(ndarray, ndarray) – The lowpass \(h\) and the quadrature-mirror highpass \(g[n] = (-1)^n h[L-1-n]\).
Examples
>>> import numpy as np >>> h, g = wavelet_filters("haar") >>> bool(np.allclose(h, [2**-0.5, 2**-0.5]) and np.allclose(g, [2**-0.5, -(2**-0.5)])) True
- mathematicskit.special_functions.z_transform(x, z)[source]#
The Z-transform \(X(z) = \sum_{n=0}^{N-1} x[n]\,z^{-n}\) of a finite sequence.
Evaluated as a polynomial in \(z^{-1}\) by Horner’s rule (
numpy.polyval()).- Parameters:
- Return type:
- Returns:
complex or ndarray
Examples
>>> complex(z_transform([1.0, 2.0, 3.0], 2.0)) # 1 + 2/2 + 3/4 (2.75+0j)