Functions

r_scale(c::Real,β::AbstractVector{<:Real},α::AbstractVector{<:Real})

Given the recursion coefficients (α,β) for a system of orthogonal polynomials that are orthogonal with respect to some positive weight $m(t)$, this function returns the recursion coefficients (α_,β_) for the scaled measure $c m(t)$ for some positive $c$.

rm_compute(weight::Function,lb::Real,ub::Real,Npoly::Int=4,Nquad::Int=10;quadrature::Function=clenshaw_curtis)

Given a positive weight function with domain (lb,ub), i.e. a function $w: [lb, ub ] \rightarrow \mathbb{R}_{\geq 0}$, this function creates Npoly recursion coefficients (α,β).

The keyword quadrature specifies what quadrature rule is being used.

rm_logistic(N::Int)

Creates N recurrence coefficients for monic polynomials that are orthogonal on $(-\infty,\infty)$ relative to $w(t) = \frac{\mathrm{e}^{-t}}{(1 - \mathrm{e}^{-t})^2}$

rm_hermite(N::Int,mu::Real)
rm_hermite(N::Int)

Creates N recurrence coefficients for monic generalized Hermite polynomials that are orthogonal on $(-\infty,\infty)$ relative to $w(t) = |t|^{2 \mu} \mathrm{e}^{-t^2}$

The call rm_hermite(N) is the same as rm_hermite(N,0).

rm_hermite_prob(N::Int)

Creates N recurrence coefficients for monic probabilists' Hermite polynomials that are orthogonal on $(-\infty,\infty)$ relative to $w(t) = \mathrm{e}^{-0.5t^2}$

rm_laguerre(N::Int,a::Real)
rm_laguerre(N::Int)

Creates N recurrence coefficients for monic generalized Laguerre polynomials that are orthogonal on $(0,\infty)$ relative to $w(t) = t^a \mathrm{e}^{-t}$.

The call rm_laguerre(N) is the same as rm_laguerre(N,0).

rm_legendre(N::Int)

Creates N recurrence coefficients for monic Legendre polynomials that are orthogonal on $(-1,1)$ relative to $w(t) = 1$.

rm_legendre01(N::Int)

Creates N recurrence coefficients for monic Legendre polynomials that are orthogonal on $(0,1)$ relative to $w(t) = 1$.

rm_jacobi(N::Int,a::Real,b::Real)
rm_jacobi(N::Int,a::Real)
rm_jacobi(N::Int)

Creates N recurrence coefficients for monic Jacobi polynomials that are orthogonal on $(-1,1)$ relative to $w(t) = (1-t)^a (1+t)^b$.

The call rm_jacobi(N,a) is the same as rm_jacobi(N,a,a) and rm_jacobi(N) the same as rm_jacobi(N,0,0).

rm_jacobi01(N::Int,a::Real,b::Real)
rm_jacobi01(N::Int,a::Real)
rm_jacobi01(N::Int)

Creates N recurrence coefficients for monic Jacobi polynomials that are orthogonal on $(0,1)$ relative to $w(t) = (1-t)^a t^b$.

The call rm_jacobi01(N,a) is the same as rm_jacobi01(N,a,a) and rm_jacobi01(N) the same as rm_jacobi01(N,0,0).

rm_meixner_pollaczek(N::Int,lambda::Real,phi::Real)
rm_meixner_pollaczek(N::Int,lambda::Real)

Creates N recurrence coefficients for monic Meixner-Pollaczek polynomials with parameters λ and ϕ. These are orthogonal on $[-\infty,\infty]$ relative to the weight function $w(t)=(2 \pi)^{-1} \exp{(2 \phi-\pi)t} |\Gamma(\lambda+ i t)|^2$.

The call rm_meixner_pollaczek(n,lambda) is the same as rm_meixner_pollaczek(n,lambda,pi/2).

stieltjes(N::Int,nodes_::AbstractVector{<:Real},weights_::AbstractVector{<:Real};removezeroweights::Bool=true)

Stieltjes procedure–-Given the nodes and weights, the function generates the firstN recurrence coefficients of the corresponding discrete orthogonal polynomials.

Set the Boolean removezeroweights to true if zero weights should be removed.

lanczos(N::Int,nodes::AbstractVector{<:Real},weights::AbstractVector{<:Real};removezeroweights::Bool=true)

Lanczos procedure–-given the nodes and weights, the function generates the first N recurrence coefficients of the corresponding discrete orthogonal polynomials.

Set the Boolean removezeroweights to true if zero weights should be removed.

The script is adapted from the routine RKPW in W.B. Gragg and W.J. Harrod, The numerically stable reconstruction of Jacobi matrices from spectral data, Numer. Math. 44 (1984), 317-335.

mcdiscretization(N::Int,quads::Vector{},discretemeasure::AbstractMatrix{<:Real}=zeros(0,2);discretization::Function=stieltjes,Nmax::Integer=300,ε::Float64=1e-8,gaussquad::Bool=false)

This routine returns $N$ recurrence coefficients of the polynomials that are orthogonal relative to a weight function $w$ that is decomposed as a sum of $m$ weights $w_i$ with domains $[a_i,b_i]$ for $i=1,\dots,m$,

\[w(t) = \sum_{i}^{m} w_i(t) \quad \text{with } \operatorname{dom}(w_i) = [a_i, b_i].\]

For each weight $w_i$ and its domain $[a_i, b_i]$ the function mcdiscretization() expects a quadrature rule of the form nodes::AbstractVector{<:Real}, weights::AbstractVector{<:Real} = myquadi(N::Int) all of which are stacked in the parameter quad quad = [ myquad1, ..., myquadm ] If the weight function has a discrete part (specified by discretemeasure) it is added on to the discretized continuous weight function.

The function mcdiscretization() performs a sequence of discretizations of the given weight $w(t)$, each discretization being followed by an application of the Stieltjes or Lanczos procedure (keyword discretization in [stieltjes, lanczos]) to produce approximations to the desired recurrence coefficients. The function applies to each subinterval $i$ an N-point quadrature rule (the $i$th entry of quad) to discretize the weight function $w_i$ on that subinterval. If the procedure converges to within a prescribed accuracy ε before N reaches its maximum allowed value Nmax. If the function does not converge, the function prompts an error message.

The keyword gaussquad should be set to true if Gauss quadrature rules are available for all$m$ weights $w_i(t)$ with $i = 1, \dots, m$.

For further information, please see W. Gautschi "Orthogonal Polynomials: Approximation and Computation", Section 2.2.4.

showpoly(coeffs::Vector{<:Real};sym::String,digits::Integer)

Show the monic polynomial with coefficients coeffs in a human-readable way. The keyword sym sets the name of the variable, and digits controls the number of shown digits.

julia> using PolyChaos

julia> showpoly([1.2, 2.3, 3.4456])
x^3 + 3.45x^2 + 2.3x + 1.2

julia> showpoly([1.2, 2.3, 3.4456], sym = "t", digits = 2)
t^3 + 3.45t^2 + 2.3t + 1.2

showpoly(d::Integer,α::Vector{<:Real},β::Vector{<:Real}; sym::String,digits::Integer)
showpoly(d::Range,α::Vector{<:Real},β::Vector{<:Real};sym::String,digits::Integer) where Range <: OrdinalRange

Show the monic polynomial of degree/range d that has the recurrence coefficients α, β.

julia> using PolyChaos

julia> α, β = rm_hermite(10)
([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.77245, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5])

julia> showpoly(3, α, β)
x^3 - 1.5x

julia> showpoly(0:2:10, α, β)
1
x^2 - 0.5
x^4 - 3.0x^2 + 0.75
x^6 - 7.5x^4 + 11.25x^2 - 1.88
x^8 - 14.0x^6 + 52.5x^4 - 52.5x^2 + 6.56
x^10 - 22.5x^8 + 157.5x^6 - 393.75x^4 + 295.31x^2 - 29.53

Tailored to types from PolyChaos.jl

showpoly(d::Union{Integer,Range},op::AbstractOrthoPoly;sym::String,digits::Integer) where Range <: OrdinalRange

Show the monic polynomial of degree/range d of an AbstractOrthoPoly.

julia> using PolyChaos

julia> op = GaussOrthoPoly(5);

julia> showpoly(3, op)
x^3 - 3.0x

julia> showpoly(0:(op.deg), op; sym = "t")
1
t
t^2 - 1.0
t^3 - 3.0t
t^4 - 6.0t^2 + 3.0
t^5 - 10.0t^3 + 15.0t

Thanks @pfitzseb for providing this functionality.

showbasis(α::Vector{<:Real},β::Vector{<:Real};sym::String,digits::Integer)

Show all basis polynomials given the recurrence coefficients α, β. The keyword sym sets the name of the variable, and digits controls the number of shown digits.

julia> using PolyChaos

julia> α, β = rm_hermite(5);

julia> showbasis(α, β)
1
x
x^2 - 0.5
x^3 - 1.5x
x^4 - 3.0x^2 + 0.75
x^5 - 5.0x^3 + 3.75x

Tailored to types from PolyChaos.jl

showbasis(op::AbstractOrthoPoly;sym::String,digits::Integer)

Show all basis polynomials of an AbstractOrthoPoly.

julia> using PolyChaos

julia> op = LegendreOrthoPoly(4);

julia> showbasis(op)
1
x
x^2 - 0.33
x^3 - 0.6x
x^4 - 0.86x^2 + 0.09

rec2coeff(deg::Int,a::Vector{<:Real},b::Vector{<:Real})
rec2coeff(a,b) = rec2coeff(length(a),a,b)

Get the coefficients of the orthogonal polynomial of degree up to deg specified by its recurrence coefficients (a,b). The function returns the values $c_i^{(k)}$ from

\[p_k (t) = t^d + \sum_{i=0}^{k-1} c_i t^i,\]

where $k$ runs from 1 to deg.

The call rec2coeff(a,b) outputs all possible recurrence coefficients given (a,b).

Univariate

evaluate(n::Int,x::Array{<:Real},a::AbstractVector{<:Real},b::AbstractVector{<:Real})
evaluate(n::Int,x::Real,a::AbstractVector{<:Real},b::AbstractVector{<:Real})
evaluate(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly)
evaluate(n::Int,x::Real,op::AbstractOrthoPoly)

Evaluate the n-th univariate basis polynomial at point(s) x The function is multiple-dispatched to facilitate its use with the composite type AbstractOrthoPoly

If several basis polynomials (stored in ns) are to be evaluated at points x, then call

evaluate(ns::AbstractVector{<:Int},x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(ns,x,op.α,op.β)
evaluate(ns::AbstractVector{<:Int},x::Real,op::AbstractOrthoPoly) = evaluate(ns,[x],op)

If all basis polynomials are to be evaluated at points x, then call

evaluate(x::AbstractVector{<:Real},op::AbstractOrthoPoly) = evaluate(collect(0:op.deg),x,op)
evaluate(x::Real,op::AbstractOrthoPoly) = evaluate([x],op)

which returns an Array of dimensions (length(x),op.deg+1).

Multivariate

evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(n::AbstractVector{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly)
evaluate(n::AbstractVector{<:Int},x::AbstractVector{<:Real},op::MultiOrthoPoly)

Evaluate the n-th p-variate basis polynomial at point(s) x The function is multiply dispatched to facilitate its use with the composite type MultiOrthoPoly

If several basis polynomials are to be evaluated at points x, then call

evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},a::Vector{<:AbstractVector{<:Real}},b::Vector{<:AbstractVector{<:Real}})
evaluate(ind::AbstractMatrix{<:Int},x::AbstractMatrix{<:Real},op::MultiOrthoPoly)

where ind is a matrix of multi-indices.

If all basis polynomials are to be evaluated at points x, then call

evaluate(x::AbstractMatrix{<:Real},mop::MultiOrthoPoly) = evaluate(mop.ind,x,mop)

which returns an array of dimensions (mop.dim,size(x,1)).

computeSP2(n::Integer,β::AbstractVector{<:Real})
computeSP2(n::Integer,op::AbstractOrthoPoly) = computeSP2(n,op.β)
computeSP2(op::AbstractOrthoPoly) = computeSP2(op.deg,op.β)

Computes the nregular scalar products aka 2-norms of the orthogonal polynomials, namely

\[\|ϕ_i\|^2 = \langle \phi_i,\phi_i\rangle \quad \forall i \in \{ 0,\dots,n \}.\]

Notice that only the values of β of the recurrence coefficients (α,β) are required. The computation is based on equation (1.3.7) from Gautschi, W. "Orthogonal Polynomials: Computation and Approximation". Whenever there exists an analytical expression for β, this function should be used.

The function is multiple-dispatched to facilitate its use with AbstractOrthoPoly.

Univariate

computeSP(a::AbstractVector{<:Integer},α::AbstractVector{<:Real},β::AbstractVector{<:Real},nodes::AbstractVector{<:Real},weights::AbstractVector{<:Real};issymmetric::Bool=false)
computeSP(a::AbstractVector{<:Integer},op::AbstractOrthoPoly;issymmetric=issymmetric(op))

Multivariate

computeSP( a::AbstractVector{<:Integer},
           α::AbstractVector{<:AbstractVector{<:Real}},β::AbstractVector{<:AbstractVector{<:Real}},
           nodes::AbstractVector{<:AbstractVector{<:Real}},weights::AbstractVector{<:AbstractVector{<:Real}},
           ind::AbstractMatrix{<:Integer};
           issymmetric::BitArray=falses(length(α)))
computeSP(a::AbstractVector{<:Integer},op::AbstractVector,ind::AbstractMatrix{<:Integer})
computeSP(a::AbstractVector{<:Integer},mOP::MultiOrthoPoly)

Computes the scalar product

\[\langle \phi_{a_1},\phi_{a_2},\cdots,\phi_{a_n} \rangle,\]

where n = length(a). This requires to provide the recurrence coefficients (α,β) and the quadrature rule (nodes,weights), as well as the multi-index ind. If provided via the keyword issymmetric, symmetry of the weight function is exploited. All computations of the multivariate scalar products resort back to efficient computations of the univariate scalar products. Mathematically, this follows from Fubini's theorem.

The function is dispatched to facilitate its use with AbstractOrthoPoly and its quadrature rule Quad.

fejer(N::Int)

Fejer's first quadrature rule.

fejer2(n::Int)

Fejer's second quadrature rule according to Waldvogel, J. Bit Numer Math (2006) 46: 195.

clenshaw_curtis(n::Int)

Clenshaw-Curtis quadrature according to Waldvogel, J. Bit Numer Math (2006) 46: 195.

quadgp(weight::Function,lb::Real,ub::Real,N::Int=10;quadrature::Function=clenshaw_curtis,bnd::Float64=Inf)

general purpose quadrature based on Gautschi, "Orthogonal Polynomials: Computation and Approximation", Section 2.2.2, pp. 93-95

Compute the N-point quadrature rule for weight with support (lb, ub). The quadrature rule can be specified by the keyword quadrature. The keyword bnd sets the numerical value for infinity.

gauss(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real})
gauss(α::AbstractVector{<:Real},β::AbstractVector{<:Real})
gauss(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})
gauss(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly})

Gauss quadrature rule, also known as Golub-Welsch algorithm

gauss() generates the N Gauss quadrature nodes and weights for a given weight function. The weight function is represented by the N recurrence coefficients for the monic polynomials orthogonal with respect to the weight function.

radau(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
radau(α::AbstractVector{<:Real},β::AbstractVector{<:Real},end0::Real)
radau(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)
radau(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},end0::Real)

Gauss-Radau quadrature rule. Given a weight function encoded by the recurrence coefficients (α,β)for the associated orthogonal polynomials, the function generates the nodes and weights (N+1)-point Gauss-Radau quadrature rule for the weight function having a prescribed node end0 (typically at one of the end points of the support interval of w, or outside thereof).

lobatto(N::Int,α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
lobatto(α::AbstractVector{<:Real},β::AbstractVector{<:Real},endl::Real,endr::Real)
lobatto(N::Int,op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)
lobatto(op::Union{OrthoPoly,AbstractCanonicalOrthoPoly},endl::Real,endr::Real)

Gauss-Lobatto quadrature rule. Given a weight function encoded by the recurrence coefficients for the associated orthogonal polynomials, the function generates the nodes and weights of the (N+2)-point Gauss-Lobatto quadrature rule for the weight function, having two prescribed nodes endl, endr (typically the left and right end points of the support interval, or points to the left resp. to the right thereof).

Univariate

mean(x::AbstractVector,op::AbstractOrthoPoly)

Multivariate

mean(x::AbstractVector,mop::MultiOrthoPoly)

compute mean of random variable with PCE x

Univariate

var(x::AbstractVector,op::AbstractOrthoPoly)
var(x::AbstractVector,t2::Tensor)

Multivariate

var(x::AbstractVector,mop::MultiOrthoPoly)
var(x::AbstractVector,t2::Tensor)

compute variance of random variable with PCE x

Univariate

std(x::AbstractVector,op::AbstractOrthoPoly)

Multivariate

std(x::AbstractVector,mop::MultiOrthoPoly)

compute standard deviation of random variable with PCE x

Univariate

sampleMeasure(n::Int,name::String,w::Function,dom::Tuple{<:Real,<:Real},symm::Bool,d::Dict;method::String="adaptiverejection")
sampleMeasure(n::Int,m::Measure;method::String="adaptiverejection")
sampleMeasure(n::Int,op::AbstractOrthoPoly;method::String="adaptiverejection")

Draw n samples from the measure m described by its

• name
• weight function w,
• domain dom,
• symmetry property symm,
• and, if applicable, parameters stored in the dictionary d. By default, an adaptive rejection sampling method is used (from AdaptiveRejectionSampling.jl), unless it is a common random variable for which Distributions.jl is used.

The function is dispatched to accept AbstractOrthoPoly.

Multivariate

sampleMeasure(n::Int,m::ProductMeasure;method::Vector{String}=["adaptiverejection" for i=1:length(m.name)])
sampleMeasure(n::Int,mop::MultiOrthoPoly;method::Vector{String}=["adaptiverejection" for i=1:length(mop.meas.name)])

Multivariate extension, which provides an array of samples with n rows and as many columns as the multimeasure has univariate measures.

evaluatePCE(x::AbstractVector{<:Real},ξ::AbstractVector{<:Real},α::AbstractVector{<:Real},β::AbstractVector{<:Real})

Evaluation of polynomial chaos expansion

\[\mathsf{x} = \sum_{i=0}^{L} x_i \phi_i{\xi_j},\]

where L+1 = length(x) and $x_j$ is the $j$th sample where $j=1,\dots,m$ with m = length(ξ).

Univariate

samplePCE(n::Int,x::AbstractVector{<:Real},op::AbstractOrthoPoly;method::String="adaptiverejection")

Combines sampleMeasure and evaluatePCE, i.e. it first draws n samples from the measure, then evaluates the PCE for those samples.

Multivariate

samplePCE(n::Int,x::AbstractVector{<:Real},mop::MultiOrthoPoly;method::Vector{String}=["adaptiverejection" for i=1:length(mop.meas.name)])

calculateAffinePCE(α::AbstractVector{<:Real})

Computes the affine PCE coefficients $x_0$ and $x_1$ from recurrence coefficients $�lpha$.

convert2affinePCE(mu::Real, sigma::Real, op::AbstractCanonicalOrthoPoly; kind::String)

Computes the affine PCE coefficients $x_0$ and $x_1$ from

\[X = a_1 + a_2 \Xi = x_0 + x_1 \phi_1(\Xi),\]

where $\phi_1(t) = t-\alpha_0$ is the first-order monic basis polynomial.

Works for subtypes of AbstractCanonicalOrthoPoly. The keyword kind in ["lbub", "μσ"] specifies whether p1 and p2 have the meaning of lower/upper bounds or mean/standard deviation.

nw(q::EmptyQuad)
nw(q::AbstractQuad)
nw(opq::AbstractOrthoPoly)
nw(opq::AbstractVector)
nw(mop::MultiOrthoPoly)

returns nodes and weights in matrix form

coeffs(op::AbstractOrthoPoly)
coeffs(op::AbstractVector)
coeffs(mop::MultiOrthoPoly)

returns recurrence coefficients of in matrix form

integrate(f::Function,nodes::AbstractVector{<:Real},weights::AbstractVector{<:Real})
integrate(f::Function,q::AbstractQuad)
integrate(f::Function,opq::AbstractOrthoPoly)

integrate function f using quadrature rule specified via nodes, weights. For example $\int_0^1 6x^5 = 1$ can be solved as follows:

julia> opq = Uniform01OrthoPoly(3) # a quadrature rule is added by default

julia> integrate(x -> 6x^5, opq)
0.9999999999999993

• function $f$ is assumed to return a scalar.
• interval of integration is "hidden" in nodes.

issymmetric(m::AbstractMeasure)
issymmetric(op::AbstractOrthoPoly)

Is the measure symmetric (around any point in the domain)?