Sequences

A Sequence is a structure representing a sequence in a prescribed VectorSpace. More precisely, a Sequence is comprised of the two fields space::VectorSpace and coefficients::AbstractVector where the dimension of the former matches the length of the latter.

julia> a = Sequence(Taylor(1), [1, 2])Sequence in Taylor(1) with coefficients Vector{Int64}: 1 2

The two fields space and coefficients are accessible via the respective functions of the same name.

julia> space(a)Taylor(1)julia> coefficients(a)2-element Vector{Int64}: 1 2

For convenience, the methods zeros, ones, fill and fill! are available:

julia> s = Taylor(1)Taylor(1)julia> zeros(s)Sequence in Taylor(1) with coefficients Vector{Float64}: 0.0 0.0julia> ones(s)Sequence in Taylor(1) with coefficients Vector{Float64}: 1.0 1.0julia> fill(2, s)Sequence in Taylor(1) with coefficients Vector{Int64}: 2 2julia> fill!(zeros(s), 2)Sequence in Taylor(1) with coefficients Vector{Float64}: 2.0 2.0

The coefficients of a Sequence are indexed according to the indices of the space (as given by indices).

julia> a[0:1] # indices(space(a))2-element Vector{Int64}: 1 2

When the space of a Sequence is a CartesianSpace, its coefficients are given as the concatenation of the coefficients associated with each space. The function component extracts a Sequence composing the cartesian space; its coefficients are a view into the parent sequence, so the two stay glued. The function unpack unglues the sequence into the Vector of its components, which lets generic code written for vectors operate on the components seamlessly.

julia> b = Sequence(ScalarSpace() × Taylor(1)^2, [1, 2, 3, 4, 5])Sequence in 𝕂 × Taylor(1)² with coefficients Vector{Int64}: 1 2 3 4 5julia> b[1:5] # indices(space(b))5-element Vector{Int64}: 1 2 3 4 5julia> component(b, 1) # extract the sequence associated with the space ScalarSpace()Sequence in 𝕂 with coefficients SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}: 1julia> component(b, 2) # extract the sequence associated with the space Taylor(1)^2Sequence in Taylor(1)² with coefficients SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}: 2 3 4 5julia> component(component(b, 2), 1)Sequence in Taylor(1) with coefficients SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}: 2 3julia> component(component(b, 2), 2)Sequence in Taylor(1) with coefficients SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}: 4 5julia> unpack(b)2-element Vector{Sequence{T, SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}} where T<:VectorSpace}: Sequence(ScalarSpace(), [1]) Sequence(Taylor(1)², [2, 3, 4, 5])

Similarly, the function eachcomponent returns a Generator whose iterates yield each Sequence composing the cartesian space.

Arithmetic

The addition and subtraction operations are implemented as the + and - functions respectively.

julia> c = Sequence(Taylor(1), [0, 1])Sequence in Taylor(1) with coefficients Vector{Int64}: 0 1julia> d = Sequence(Taylor(2), [1, 2, 1])Sequence in Taylor(2) with coefficients Vector{Int64}: 1 2 1julia> c + dSequence in Taylor(2) with coefficients Vector{Int64}: 1 3 1julia> c - dSequence in Taylor(2) with coefficients Vector{Int64}: -1 -1 -1

The discrete convolution between sequences whose spaces are a SequenceSpace is implemented as the *, mul! and ^ functions. Their bar counterparts mul_bar (unicode alias *\bar<tab>) and pow_bar (unicode alias ^\bar<tab>) give the result projected in the smallest compatible space between the operands; in general, mul_bar is not associative.

julia> c * dSequence in Taylor(3) with coefficients Vector{Int64}: 0 1 2 1julia> c ^ 3Sequence in Taylor(3) with coefficients Vector{Int64}: 0 0 0 1julia> mul_bar(c, d) # project(c * d, Taylor(1))Sequence in Taylor(1) with coefficients Vector{Int64}: 0 1julia> pow_bar(c, 3) # project(c ^ 3, Taylor(1))Sequence in Taylor(1) with coefficients Vector{Int64}: 0 0

To improve performance, the FFT algorithm may be used to compute discrete convolutions via the Convolution Theorem. However, the performance gain is tempered with the loss of accuracy which may stop the decay of the coefficients. To circumvent machine precision limitations, the coefficients beyond a prescribed order are enclosed rigorously by a Banach-algebra estimate.^[1]

julia> x = Sequence(Taylor(3), interval.([inv(10_000.0 ^ i) for i ∈ 0:3]))Sequence in Taylor(3) with coefficients Vector{Interval{Float64}}: [1.0, 1.0]_com [0.0001, 0.0001]_com [1.0e-8, 1.0e-8]_com [1.0e-12, 1.0e-12]_comjulia> set_conv_algorithm(:fft):fftjulia> x³ = x ^ 3Sequence in Taylor(9) with coefficients Vector{Interval{Float64}}:  [0.999999, 1.00001]_com  [0.000299999, 0.000300001]_com  [5.99999e-8, 6.00001e-8]_com  [9.99947e-12, 1.00006e-11]_com  [9.05207e-16, 1.5993e-15]_com [-6.40001e-19, 6.40001e-19]_com_NG [-6.40001e-23, 6.40001e-23]_com_NG [-6.4e-27, 6.4e-27]_com_NG [-6.4e-31, 6.4e-31]_com_NG [-6.4e-35, 6.4e-35]_com_NGjulia> set_conv_algorithm(:loop) # default algorithm:loopjulia> x³ = x ^ 3 # only rounding errorSequence in Taylor(9) with coefficients Vector{Interval{Float64}}: [1.0, 1.0]_com [0.0003, 0.0003]_com [6.0e-8, 6.00001e-8]_com [1.0e-11, 1.00001e-11]_com [1.19999e-15, 1.20001e-15]_com [1.19999e-19, 1.20001e-19]_com [9.9999e-24, 1.00001e-23]_com [6.0e-28, 6.00001e-28]_com [2.99999e-32, 3.00001e-32]_com [9.9999e-37, 1.00001e-36]_com

Grids and interpolation

The functions to_grid and to_coef convert between coefficient and grid space, using rigorous FFTs when the coefficients are intervals. grid_size(space) gives the number of nodes, as a tuple with one entry per space. The nodes are:

• for Taylor: the roots of unity $e^{\mathrm{i} 2\pi j / (m-1)}$ for $j = 0, ..., m-1$.
• for Fourier (with frequency $\omega$): the equispaced points of the period $2\pi/\omega j / (m-1)$ for $j = 0, ..., m-1$.
• for Chebyshev: the Chebyshev–Lobatto nodes $\cos(\pi j /(m-1))$ for $j = 0, ..., m-1$ (ordered from $1$ down to $-1$).

It is the smallest grid that determines the space, so to_coef(f, space) returns the interpolant of f at those nodes, and to_coef(to_grid(a), space(a)) recovers a.

julia> m = grid_size(Chebyshev(2)) # one entry per factor(3,)julia> g = [cospi((k-1)/(m[1]-1)) for k ∈ 1:m[1]] # values of f(x) = x at the nodes3-element Vector{Float64}:  1.0  0.0 -1.0julia> to_coef(g, Chebyshev(2)) # x = 2 (0.5 T₁(x)) due to the normalizationSequence in Chebyshev(2) with coefficients Vector{ComplexF64}: 0.0 + 0.0im 0.5 + 0.0im 0.0 + 0.0imjulia> to_grid(ans, m) # back to the nodes3-element Vector{ComplexF64}:  1.0 - 0.0im  0.0 + 0.0im -1.0 - 0.0im

Oversampling a grid keeps the interpolant unchanged. Giving fewer sizes than a TensorSpace has factors discretizes only the leading factors lead_space, so that to_grid(a, m) then returns a grid of Sequences on the remaining factors, and to_coef(x_grid, s) interpolates such a grid back into a Sequence on lead_space ⊗ inner_space.

julia> a = Sequence(Chebyshev(2) ⊗ Fourier(1, 1.0), collect(1:9)) # a family of Fourier sequencesSequence in Chebyshev(2) ⊗ Fourier(1, 1.0) with coefficients Vector{Int64}: 1 2 3 4 5 6 7 8 9julia> x_grid = to_grid(a, grid_size(Chebyshev(2))) # one Fourier sequence per Chebyshev–Lobatto node3-element Vector{Sequence{Fourier{Float64}, Vector{ComplexF64}}}: Sequence(Fourier(1, 1.0), ComplexF64[11.0 - 0.0im, 26.0 - 0.0im, 41.0 - 0.0im]) Sequence(Fourier(1, 1.0), ComplexF64[-5.0 + 0.0im, -8.0 + 0.0im, -11.0 + 0.0im]) Sequence(Fourier(1, 1.0), ComplexF64[3.0 - 0.0im, 6.0 - 0.0im, 9.0 - 0.0im])julia> to_coef(x_grid, Chebyshev(2)) # interpolate backSequence in Chebyshev(2) ⊗ Fourier(1, 1.0) with coefficients Vector{ComplexF64}: 1.0 + 0.0im 2.0 + 0.0im 3.0 + 0.0im 4.0 + 0.0im 5.0 + 0.0im 6.0 + 0.0im 7.0 + 0.0im 8.0 + 0.0im 9.0 + 0.0im

A CartesianSpace is discretized componentwise, every component sharing the nodes, so that to_grid(a, m) returns a grid of Sequences in the cartesian space of the remaining factors, and to_coef(x_grid, s) distributes s over the components. The same holds for a LinearOperator whose domain or codomain is cartesian, giving a grid of operators between the inner spaces.

julia> a = Sequence((Chebyshev(2) ⊗ Fourier(1, 1.0)) × Chebyshev(2), collect(1:12)) # two familiesSequence in (Chebyshev(2) ⊗ Fourier(1, 1.0)) × Chebyshev(2) with coefficients Vector{Int64}:  1  2  3  4  5  6  7  8  9 10 11 12julia> x_grid = to_grid(a, grid_size(Chebyshev(2))) # the second component is left with no factor3-element Vector{Sequence{CartesianProduct{Tuple{Fourier{Float64}, ScalarSpace}}, Vector{ComplexF64}}}: Sequence(Fourier(1, 1.0) × ScalarSpace(), ComplexF64[11.0 - 0.0im, 26.0 - 0.0im, 41.0 - 0.0im, 56.0 - 0.0im]) Sequence(Fourier(1, 1.0) × ScalarSpace(), ComplexF64[-5.0 + 0.0im, -8.0 + 0.0im, -11.0 + 0.0im, -14.0 + 0.0im]) Sequence(Fourier(1, 1.0) × ScalarSpace(), ComplexF64[3.0 - 0.0im, 6.0 - 0.0im, 9.0 - 0.0im, 12.0 - 0.0im])julia> to_coef(x_grid, Chebyshev(2)) # interpolate backSequence in (Chebyshev(2) ⊗ Fourier(1, 1.0)) × Chebyshev(2) with coefficients Vector{ComplexF64}:  1.0 + 0.0im  2.0 + 0.0im  3.0 + 0.0im  4.0 + 0.0im  5.0 + 0.0im  6.0 + 0.0im  7.0 + 0.0im  8.0 + 0.0im  9.0 + 0.0im 10.0 + 0.0im 11.0 + 0.0im 12.0 + 0.0im

API

AbstractSequence

Sequence{T<:VectorSpace,S<:AbstractVector} <: AbstractSequence

julia> Sequence(Taylor(2), [1, 2, 1]) # 1 + 2x + x^2Sequence in Taylor(2) with coefficients Vector{Int64}: 1 2 1julia> Sequence(Taylor(1) ⊗ Fourier(1, 1.0), [0.5, 0.5, 0.0, 0.0, 0.5, 0.5]) # (1 + x) cos(y)Sequence in Taylor(1) ⊗ Fourier(1, 1.0) with coefficients Vector{Float64}: 0.5 0.5 0.0 0.0 0.5 0.5julia> Sequence([1, 2, 3])Sequence in 𝕂³ with coefficients Vector{Int64}: 1 2 3

algebraicweight(a::Sequence{<:SequenceSpace})

julia> rate(algebraicweight(Sequence(Taylor(10), [inv((1.0 + i)^2) for i in 0:10]))) ≈ 2truejulia> rate.(algebraicweight(Sequence(Taylor(10) ⊗ Fourier(3, 1.0), vec([inv((1.0 + i)^2 * (1.0 + abs(j))^3) for i in 0:10, j in -3:3])))) .≈ (2, 3)(true, true)

component(a::Sequence{<:CartesianSpace}, i)

eachcomponent(a::Sequence{<:CartesianSpace})

geometricweight(a::Sequence{<:SequenceSpace})

julia> rate(geometricweight(Sequence(Taylor(10), [inv(2.0^i) for i in 0:10]))) ≈ 2truejulia> rate.(geometricweight(Sequence(Taylor(10) ⊗ Fourier(3, 1.0), vec([inv(2.0^i * 3.0^abs(j)) for i in 0:10, j in -3:3])))) .≈ (2, 3)(true, true)

unpack(a::Sequence{<:CartesianSpace})

InfiniteSequence{T<:SequenceSpace,S<:AbstractVector,R<:Real,U<:BanachSpace} <: AbstractSequence

julia> InfiniteSequence(Sequence(Taylor(2), [1.0, 2.0, 1.0]), 0.0, 0.1, 0.1, Ell1())Sequence in Taylor(2) with coefficients Vector{Float64}: 1.0 2.0 1.0Norm of the truncated sequence: 4.0Finite error: 0.0Tail error: 0.1Total error: 0.1Banach space: ℓ¹()

fft_size(s::SequenceSpace)
fft_size(s::CartesianSpace)

julia> fft_size(Chebyshev(2)), grid_size(Chebyshev(2))((4,), (3,))julia> fft_size(Chebyshev(2) × Chebyshev(4)) # the coarser component is oversampled(8,)

grid_size(s::SequenceSpace)
grid_size(s::CartesianSpace)

julia> grid_size(Taylor(2)), grid_size(Fourier(2, 1.0)), grid_size(Chebyshev(2))((3,), (5,), (3,))julia> grid_size(Chebyshev(2)^3) # the three components share the nodes(3,)

to_coef!(c::Sequence, x_grid)

to_coef(f::Function, s::SequenceSpace)
to_coef(a::Sequence, s::SequenceSpace)
to_coef(x_grid::AbstractArray, s::SequenceSpace)

to_grid!(x_grid, a::Sequence)

to_grid(a::Sequence, m = grid_size(space(a)))

*(a::Sequence{<:SequenceSpace}, b::Sequence{<:SequenceSpace})

^(a::Sequence{<:SequenceSpace}, n::Int)