# README
Gosl. fun. Special functions, DFT, FFT, Bessel, elliptical integrals, orthogonal polynomials, interpolators
This package implements special functions such as orthogonal polynomials and elliptical functions of first, second and third kind.
Routines to interpolate and/or assist on spectral methods are also available; e.g. FourierInterp, ChebyInterp.
API
# Functions
Atan2p implements a positive version of atan2, in such a way that: 0 ≤ α ≤ 2π.
Atan2pDeg implements a positive version of atan2, in such a way that: 0 ≤ α ≤ 360.
Beta computes the beta function by calling the Lgamma function.
Binomial comptues the binomial coefficient (n k)^T.
Boxcar implements the boxcar function
Boxcar(x;a,b) = Heav(x-a) - Heav(x-b)
│ 0 if x < a or x > b Boxcar(x;a,b) = ┤ 1/2 if x = a or x = b │ 1 if x > a and x < b
Note: a ≤ x ≤ b; i.e.
CarlsonRc computes Carlson’s degenerate elliptic integral according to [1] Computes Rc(x,y) where x must be nonnegative and y must be nonzero.
CarlsonRd computes Carlson’s elliptic integral of the second kind according to [1] Computes Rf(x,y,z) where x,y must be non-negative and at most one can be zero.
CarlsonRf computes Carlson's elliptic integral of the first kind according to [1].
CarlsonRj computes Carlson’s elliptic integral of the third kind according to [1] Computes Rj(x,y,z,p) where x,y,z must be nonnegative, and at most one can be zero.
ChebyshevT directly computes the Chebyshev polynomial of first kind Tn(x) using the trigonometric functions.
ChebyshevTdiff1 computes the first derivative of the Chebyshev function Tn(x)
dTn ———— dx
.
ChebyshevTdiff2 computes the second derivative of the Chebyshev function Tn(x)
d²Tn ————— dx²
.
ChebyshevXgauss computes Chebyshev-Gauss roots considering symmetry
/ (2i+1)⋅π \ X[i] = -cos | —————————— | i = 0 ..
ChebyshevXlob computes Chebyshev-Gauss-Lobatto points using the sin function and considering symmetry
/ π⋅(N-2i) \ X[i] = -sin | —————————— | i = 0 ..
Dft1d computes the discrete Fourier transform (DFT) in 1D.
Elliptic1 computes Legendre elliptic integral of the first kind F(φ,k), evaluated using Carlson’s function Rf [1].
Elliptic2 computes Legendre elliptic integral of the second kind E(φ,k), evaluated using Carlson's functions Rf and Rd [1].
Elliptic3 computes Legendre elliptic integral of the third kind Π(n,φ,k), evaluated using Carlson's functions Rf and Rj.
ExpMix uses Euler's formula to compute exp(-i⋅x) = cos(x) - i⋅sin(x).
ExpPix uses Euler's formula to compute exp(+i⋅x) = cos(x) + i⋅sin(x).
Factorial100 returns the factorial n! up to 100! using the math/big package.
Factorial22 implements the factorial function; i.e.
Hat implements the hat function
--———-- o (xc,y0+h) | / \ h / \ m = h/l | /m \ y0 ——————o o—————————
|< 2l >|
.
HatD1 returns the first derivative of the hat function NOTE: the discontinuity is ignored ⇒ D1(xc-l)=D1(xc+l)=D1(xc)=0.
Heav computes the Heaviside step function (== derivative of Ramp(x))
│ 0 if x < 0 Heav(x) = ┤ 1/2 if x = 0 │ 1 if x > 0
.
ImagPowN computes iⁿ = (√-1)ⁿ
i¹ = i i² = -1 i³ = -i i⁴ = 1 i⁵ = i i⁶ = -1 i⁷ = -i i⁸ = 1 i⁹ = i i¹⁰ = -1 i¹¹ = -i i¹² = 1
.
ImagXpowN computes (x⋅i)ⁿ
(x⋅i)¹ = x¹⋅i (x⋅i)² = -x² (x⋅i)³ = -x³ ⋅i (x⋅i)⁴ = x⁴ (x⋅i)⁵ = x⁵⋅i (x⋅i)⁶ = -x⁶ (x⋅i)⁷ = -x⁷ ⋅i (x⋅i)⁸ = x⁸ (x⋅i)⁹ = x⁹⋅i (x⋅i)¹⁰ = -x¹⁰ (x⋅i)¹¹ = -x¹¹⋅i (x⋅i)¹² = x¹²
.
Logistic implements the sigmoid/logistic function.
LogisticD1 implements the first derivative of the sigmoid/logistic function.
ModBesselI0 returns the modified Bessel function I0(x) for any real x.
ModBesselI1 returns the modified Bessel function I1(x) for any real x.
ModBesselIn returns the modified Bessel function In(x) for any real x and n ≥ 0.
ModBesselK0 returns the modified Bessel function K0(x) for positive real x.
ModBesselK1 returns the modified Bessel function K1(x) for positive real x.
ModBesselKn returns the modified Bessel function Kn(x) for positive x and n ≥ 0.
NegOnePowN computes (-1)ⁿ.
NewAxis builds a new Axis type from a data slice for an InterpType.
NewBiLinear builds a two dimensional bi-linear interpolant Input: xx -- function sample points abscissas yy -- function sample points ordinates f -- function values f(i,j) is stored at f[len(xx)*j + i]
Ref: f(x,y) = x^2 + 2y^2 2 h i j xx = [0.00,0.50,1.00] | yy = [0.00,1.00,2.00] 1 d e f f = [a:0.00,b:0.25,c:1.00, | d:2.00,e:2,25,f:3.00, a____b____c h:8.00,i:8.25,j:9.00] 0 0.5 1.0.
NewChebyInterp returns a new ChebyInterp structure
gaussChebyshev == true:
/ (2⋅j+1)⋅π \ X_j = cos | ——————————— | j = 0 ..
NewDataInterp creates new interpolator for data point sets xx and yy (with same lengths)
Type -- type of interpolator "lin" : linear "poly" : polynomial
p -- order of interpolator xx -- x-data yy -- y-data.
NewFourierInterp allocates a new FourierInterp object
N -- number of terms.
NewGeneralOrthoPoly creates a new orthogonal polynomial
kind -- is the type of orthognal polynomial: "J" or "jac" : Jacobi "L" or "leg" : Legendre "H" or "her" : Hermite "T" or "cheby1" : Chebyshev first kind "U" or "cheby2" : Chebyshev second kind
N -- is the (max) degree of the polynomial.
NewInterpCubic returns a new object.
NewInterpQuad returns a new object.
NewLagIntSet returns a set of LagrangeInterp.
NewLagrangeInterp allocates a new LagrangeInterp
N -- degree
gridType -- type of grid: "uni" : uniform 1D grid kind "cg" : Chebyshev-Gauss 1D grid kind "cgl" : Chebyshev-Gauss-Lobatto 1D grid kind
NOTE: the grid will be generated in [-1, 1]
.
NewSinusoidBasis creates a new Sinusoid object with the "basis" parameters set T -- period; e.g.
NewSinusoidEssential creates a new Sinusoid object with the "essential" parameters set T -- period; e.g.
Pow2 computes x².
Pow3 computes x³.
PowP computes real raised to positive integer xⁿ.
Ramp function => MacAulay brackets.
Rbinomial computes the binomial coefficient with real (non-negative) arguments by calling the Gamma function.
Rect implements the rectangular function
Rect(x) = Boxcar(x;-0.5,0.5)
│ 0 if |x| > 1/2 Rect(x) = ┤ 1/2 if |x| = 1/2 │ 1 if |x| < 1/2
.
Sabs implements a smooth abs function: Sabs(x) = x*x / (sign(x)*x + eps).
SabsD1 returns the first derivative of Sabs.
SabsD2 returns the first derivative of Sabs.
Sign implements the sign function
│ -1 if x < 0 Sign(x) = ┤ 0 if x = 0 │ 1 if x > 0
.
Sinc computes the sine cardinal (sinc) function
Sinc(x) = | 1 if x = 0 | sin(x)/x otherwise
.
Sramp implements a smooth ramp function.
SrampD1 returns the first derivative of Sramp.
SrampD2 returns the second derivative of Sramp.
SuqCos implements the superquadric auxiliary function that uses cos(x).
SuqSin implements the superquadric auxiliary function that uses sin(x).
UintBinomial implements the Binomial coefficient using uint64.
# Constants
BiCubicType defines the bi-cubic type.
BiLinearType defines the bi-linear type.
π = 3.141592653589...
# Structs
Axis implements a type to hold an arbitrarily spaced discrete data.
BiLinear implements a two dimensional interpolant.
ChebyInterp defines a structure for efficient computations with Chebyshev polynomials such as projection or interpolation
Some equations are based on [1,2,3]
References: [1] Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in Single Domains.
DataInterp implements numeric interpolators to be used with discrete data.
FourierInterp performs interpolation using truncated Fourier series
N/2 - 1 ———— +i k X[j] f(x[j]) = \ A[k] ⋅ e with X[j] = 2 π j / N / ———— k = -N/2 Eq (2.1.27) of [1] x ϵ [0, 2π]
where:
N - 1 1 ———— -i k X[j] A[k] = ——— \ f(x[j]) ⋅ e with X[j] = 2 π j / N N / ———— j = 0 Eq (2.1.25) of [1]
NOTE: (1) f=u in [1] and A[k] is the tilde(u[k]) of [1] (2) FFTW says "so you should initialize your input data after creating the plan." Therefore, the plan can be created and reused several times.
GeneralOrthoPoly implements general orthogonal polynomials.
InterpCubic computes a cubic polynomial to perform interpolation either using 4 points or 3 points and a known derivative.
InterpQuad computes a quadratic polynomial to perform interpolation either using 3 points or 2 points and a known derivative.
LagrangeInterp implements Lagrange interpolators associated with a grid X
An interpolant I^X_N{f} (associated with a grid X; of degree N; with N+1 points) is expressed in the Lagrange form as follows:
N X ———— X I {f}(x) = \ f(x[i]) ⋅ ℓ (x) N / i ———— i = 0
where ℓ^X_i(x) is the i-th Lagrange cardinal polynomial associated with grid X and given by:
N N ━━━━ x - X[j] ℓ (x) = ┃ ┃ ————————————— 0 ≤ i ≤ N i ┃ ┃ X[i] - X[j] j = 0 j ≠ i
or, barycentric form:
N λ[i] ⋅ f[i] Σ ——————————— X i=0 x - x[i] I {f}(x) = —————————————————— N N λ[i] Σ ———————— i=0 x - x[i]
with:
λ[i] ———————— N x - x[i] ℓ (x) = ——————————————— i N λ[k] Σ ———————— k=0 x - x[k]
The barycentric weights λk are normalized and computed from ηk as follows:
ηk = Σ ln(|xk-xl|) (k≠l)
a ⋅ b k+N λk = ————— a = (-1) b = exp(m) m = -ηk lf0
lf0 = 2ⁿ⁻¹/n
or, if N > 700:
/ a ⋅ b \ / b \ / b \ λk = | ————— | ⋅ | ——— | ⋅ | ——— | b = exp(m/3) \ lf0 / \ lf1 / \ lf2 /
lf0⋅lf1⋅lf2 = 2ⁿ⁻¹/n
References: [1] Canuto C, Hussaini MY, Quarteroni A, Zang TA (2006) Spectral Methods: Fundamentals in Single Domains.
Sinusoid implements the sinusoid equation:
y(t) = A0 + C1⋅cos(ω0⋅t + θ) [essential-form]
y(t) = A0 + A1⋅cos(ω0⋅t) + B1⋅sin(ω0⋅t) [basis-form]
A1 = C1⋅cos(θ) B1 = -C1⋅sin(θ) θ = arctan(-B1 / A1) if A1<0, θ += π C1 = sqrt(A1² + B1²)
.
# Type aliases
InterpType specifies the type of interpolant.
LagIntSet is groups interpolators together; e.g.
Mm defines a matrix function f(m) of a matrix argument m (matrix matrix)) Input: m -- input matrix Output: M -- output matrix.
Mss defines a matrix function f(a,b) of two scalar arguments (matrix scalar scalar) Input: a -- first input scalar b -- second input scalar Output: m -- output matrix.
Mv defines a matrix function f(v) of a vector argument v (matrix vector)) Input: v -- input vector Output: f -- output matrix.
Ss defines a scalar function f(s) of a scalar argument s (scalar scalar) Input: s -- input scalar Returns: scalar.
Sss defines a scalar function f(r,s) of two scalar arguments (scalar scalar scalar) Input: s1, s2 -- input scalar Returns: scalar.
Sv defines a scalar function f(v) of a vector argument v (scalar vector) Input: v -- input vector Returns: scalar.
Svs defines a scalar function f(v,s) of a vector and a scalar Input: s -- the scalar v -- the vector Returns: scalar.
Tt defines a triplet (matrix) function f(t) of a triplet (matrix) argument t (triplet triplet) Input: t -- input triplet Output: f -- output triplet.
Tv defines a triplet (matrix) function f(v) of a vector argument v (triplet vector) Input: v -- input vector Output: f -- output triplet.
Vs defines a vector function f(s) of a scalar argument s (vector scalar) Input: s -- input scalar Output: f -- output vector.
Vss defines a vector function f(a,b) of two scalar arguments (vector scalar scalar) Input: a -- first input scalar b -- second input scalar Output: f -- output vector.
Vv defines a vector function f(v) of a vector argument v (vector vector) Input: v -- input vector Output: f -- output vector.
Vvss defines two vector functions u(a,b) and v(a,b) of two scalar arguments (vector vector scalar scalar) Input: a -- first input scalar b -- second input scalar Output: u -- first output vector v -- second output vector.
Vvvss defines three vector functions u(a,b), v(a,b) and w(a,b) of two scalar arguments (vector vector vector scalar scalar) Input: a -- first input scalar b -- second input scalar Output: u -- first output vector v -- second output vector w -- second output vector.