Computation of Scalar Products
By now, we are able to construct orthogonal polynomials, and to construct quadrature rules for a given nonnegative weight function, respectively. Now we combine both ideas to solve integrals involving the orthogonal polynomials
\[\langle \phi_{i_1} \phi_{i_2} \cdots \phi_{i_{m-1}}, \phi_{i_m} \rangle = \int \phi_{i_1}(t) \phi_{i_2}(t) \cdots \phi_{i_{m-1}}(t) \phi_{i_m}(t) w(t) \mathrm{d} t,\]
both for the univariate and multivariate case. The integrand is a polynomial (possibly multivariate) that can be solved exactly with the appropriate Gauss quadrature rules.
To simplify notation we drop the integration interval. It is clear from the context.
Univariate Polynomials
Classical Polynomials
Let's begin with a univariate basis for some classical orthogonal polynomial
using PolyChaos
deg, n = 4, 20
s_α, s_β = 2.1, 3.2
opq = Beta01OrthoPoly(deg, s_α, s_β; Nrec=n, addQuadrature=true)
Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}(4, [0.39622641509433965, 0.45308865339881105, 0.4732655766681396, 0.482729089351984, 0.4879233481934926, 0.49108064278342917, 0.49314292190864784, 0.4945640779897571, 0.4955849084142159, 0.4963428640512603, 0.49692106736331404, 0.4973721930243669, 0.49773093800397555, 0.4980209139779084, 0.4982586420286101, 0.4984559630451982, 0.4986215434156863, 0.4987618443402984, 0.4988817625639983, 0.4989850639437675], [1.0, 0.03797318144060756, 0.049528193464418745, 0.054496397032519545, 0.05707494398498378, 0.05858224318366976, 0.05953877522858778, 0.06018346810103027, 0.0606384719605267, 0.06097150420824787, 0.061222558056958205, 0.06141648535528466, 0.061569388295544744, 0.061692071403646374, 0.06179200321800959, 0.061874480034347076, 0.06194334217794196, 0.06200142915302492, 0.06205087705273863, 0.062093317746215744], Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.00878890525374891, 0.028450106116514506, 0.05852125536671723, 0.09833813997192814, 0.14702350691009355, 0.20350504892083368, 0.266538831062297, 0.3347366435399401, 0.40659656306145764, 0.48053602819819125, 0.5549266957500689, 0.6281303108338084, 0.6985348047519305, 0.7645898365922597, 0.8248410256009169, 0.8779622119533655, 0.9227853738738273, 0.9583291958712525, 0.9838388819882817], [0.00115167784765143, 0.007007259624729887, 0.020314892689531894, 0.04140255942833847, 0.06777171117046146, 0.09471990538153915, 0.11673118518989155, 0.12913645448496422, 0.12947566141190417, 0.11812235962361355, 0.09800772591714006, 0.07359961541684255, 0.04953526816082873, 0.029389573634464723, 0.014963964676256315, 0.006252710931950127, 0.0019796669915308576, 0.00040351473743713024, 3.429268092282879e-5]))
By setting addQuadrature = true
(which is default), an $n$-point Gauss quadrature rule is create relative to the underlying measure opq.measure
, where $n$ is the number of recurrence coefficients stored in opq.α
and opq.β
.
To compute the squared norms
\[\| \phi_k \|^2 = \langle \phi_k, \phi_k \rangle = \int \phi_k(t) \phi_k(t) w(t) \mathrm{d} t\]
of the basis we call computeSP2()
normsq = computeSP2(opq)
5-element Vector{Float64}:
1.0
0.03797318144060756
0.0018807430768498865
0.00010249372143217383
5.849823409553853e-6
For the general case
\[\langle \phi_{i_1} \phi_{i_2} \cdots \phi_{i_{m-1}}, \phi_{i_m} \rangle = \int \phi_{i_1}(t) \phi_{i_2}(t) \cdots \phi_{i_{m-1}}(t) \phi_{i_m}(t) w(t) \mathrm{d} t,\]
there exists a type Tensor
that requires only two arguments: the dimension $m \geq 1$, and an AbstractOrthoPoly
m = 3
t = Tensor(3,opq)
Tensor{Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}}(3, [1 ] = 1.0
[6 ] = 0.0379732
[11] = 0.00188074
[16] = 0.000102494
[21] = 5.84982e-6
[22] = 0.00215924
[23] = 0.00188074
⋮
[45] = 5.84982e-6
[48] = 7.80614e-6
[49] = 7.09802e-7
[53] = 4.73123e-7
[64] = 9.40423e-7
[65] = 5.08904e-7
[69] = 6.53232e-8
[85] = 3.55404e-8, PolyChaos.var"#getfun#45"{Int64, Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}, SparseArrays.SparseVector{Float64, Int64}}(3, Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}(4, [0.39622641509433965, 0.45308865339881105, 0.4732655766681396, 0.482729089351984, 0.4879233481934926, 0.49108064278342917, 0.49314292190864784, 0.4945640779897571, 0.4955849084142159, 0.4963428640512603, 0.49692106736331404, 0.4973721930243669, 0.49773093800397555, 0.4980209139779084, 0.4982586420286101, 0.4984559630451982, 0.4986215434156863, 0.4987618443402984, 0.4988817625639983, 0.4989850639437675], [1.0, 0.03797318144060756, 0.049528193464418745, 0.054496397032519545, 0.05707494398498378, 0.05858224318366976, 0.05953877522858778, 0.06018346810103027, 0.0606384719605267, 0.06097150420824787, 0.061222558056958205, 0.06141648535528466, 0.061569388295544744, 0.061692071403646374, 0.06179200321800959, 0.061874480034347076, 0.06194334217794196, 0.06200142915302492, 0.06205087705273863, 0.062093317746215744], Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.00878890525374891, 0.028450106116514506, 0.05852125536671723, 0.09833813997192814, 0.14702350691009355, 0.20350504892083368, 0.266538831062297, 0.3347366435399401, 0.40659656306145764, 0.48053602819819125, 0.5549266957500689, 0.6281303108338084, 0.6985348047519305, 0.7645898365922597, 0.8248410256009169, 0.8779622119533655, 0.9227853738738273, 0.9583291958712525, 0.9838388819882817], [0.00115167784765143, 0.007007259624729887, 0.020314892689531894, 0.04140255942833847, 0.06777171117046146, 0.09471990538153915, 0.11673118518989155, 0.12913645448496422, 0.12947566141190417, 0.11812235962361355, 0.09800772591714006, 0.07359961541684255, 0.04953526816082873, 0.029389573634464723, 0.014963964676256315, 0.006252710931950127, 0.0019796669915308576, 0.00040351473743713024, 3.429268092282879e-5])), [1 ] = 1.0
[6 ] = 0.0379732
[11] = 0.00188074
[16] = 0.000102494
[21] = 5.84982e-6
[22] = 0.00215924
[23] = 0.00188074
⋮
[45] = 5.84982e-6
[48] = 7.80614e-6
[49] = 7.09802e-7
[53] = 4.73123e-7
[64] = 9.40423e-7
[65] = 5.08904e-7
[69] = 6.53232e-8
[85] = 3.55404e-8), Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}(4, [0.39622641509433965, 0.45308865339881105, 0.4732655766681396, 0.482729089351984, 0.4879233481934926, 0.49108064278342917, 0.49314292190864784, 0.4945640779897571, 0.4955849084142159, 0.4963428640512603, 0.49692106736331404, 0.4973721930243669, 0.49773093800397555, 0.4980209139779084, 0.4982586420286101, 0.4984559630451982, 0.4986215434156863, 0.4987618443402984, 0.4988817625639983, 0.4989850639437675], [1.0, 0.03797318144060756, 0.049528193464418745, 0.054496397032519545, 0.05707494398498378, 0.05858224318366976, 0.05953877522858778, 0.06018346810103027, 0.0606384719605267, 0.06097150420824787, 0.061222558056958205, 0.06141648535528466, 0.061569388295544744, 0.061692071403646374, 0.06179200321800959, 0.061874480034347076, 0.06194334217794196, 0.06200142915302492, 0.06205087705273863, 0.062093317746215744], Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.00878890525374891, 0.028450106116514506, 0.05852125536671723, 0.09833813997192814, 0.14702350691009355, 0.20350504892083368, 0.266538831062297, 0.3347366435399401, 0.40659656306145764, 0.48053602819819125, 0.5549266957500689, 0.6281303108338084, 0.6985348047519305, 0.7645898365922597, 0.8248410256009169, 0.8779622119533655, 0.9227853738738273, 0.9583291958712525, 0.9838388819882817], [0.00115167784765143, 0.007007259624729887, 0.020314892689531894, 0.04140255942833847, 0.06777171117046146, 0.09471990538153915, 0.11673118518989155, 0.12913645448496422, 0.12947566141190417, 0.11812235962361355, 0.09800772591714006, 0.07359961541684255, 0.04953526816082873, 0.029389573634464723, 0.014963964676256315, 0.006252710931950127, 0.0019796669915308576, 0.00040351473743713024, 3.429268092282879e-5])))
To get the desired entries, Tensor
comes with a get()
function that is called for some index $a \in \mathbb{N}_0^m$ that has the entries $a = [i_1, i_2, \dots, i_m]$. For example
t.get([1,2,3])
0.00010249372143217402
Or using comprehension
T = [ t.get([i1,i2,i3]) for i1=0:dim(opq)-1, i2=0:dim(opq)-1, i3=0:dim(opq)-1]
5×5×5 Array{Float64, 3}:
[:, :, 1] =
1.0 0.0 0.0 0.0 0.0
0.0 0.0379732 0.0 0.0 0.0
0.0 0.0 0.00188074 0.0 0.0
0.0 0.0 0.0 0.000102494 0.0
0.0 0.0 0.0 0.0 5.84982e-6
[:, :, 2] =
0.0 0.0379732 0.0 0.0 0.0
0.0379732 0.00215924 0.00188074 0.0 0.0
0.0 0.00188074 0.000144891 0.000102494 0.0
0.0 0.0 0.000102494 8.86598e-6 5.84982e-6
0.0 0.0 0.0 5.84982e-6 5.36411e-7
[:, :, 3] =
0.0 0.0 0.00188074 0.0 0.0
0.0 0.00188074 0.000144891 0.000102494 0.0
0.00188074 0.000144891 0.000127149 1.0934e-5 5.84982e-6
0.0 0.000102494 1.0934e-5 7.80614e-6 7.09802e-7
0.0 0.0 5.84982e-6 7.09802e-7 4.73123e-7
[:, :, 4] =
0.0 0.0 0.0 0.000102494 0.0
0.0 0.0 0.000102494 8.86598e-6 5.84982e-6
0.0 0.000102494 1.0934e-5 7.80614e-6 7.09802e-7
0.000102494 8.86598e-6 7.80614e-6 9.40423e-7 5.08904e-7
0.0 5.84982e-6 7.09802e-7 5.08904e-7 6.53232e-8
[:, :, 5] =
0.0 0.0 0.0 0.0 5.84982e-6
0.0 0.0 0.0 5.84982e-6 5.36411e-7
0.0 0.0 5.84982e-6 7.09802e-7 4.73123e-7
0.0 5.84982e-6 7.09802e-7 5.08904e-7 6.53232e-8
5.84982e-6 5.36411e-7 4.73123e-7 6.53232e-8 3.55404e-8
Notice that we can cross-check the results.
using LinearAlgebra
normsq == diag(T[:,:,1]) == diag(T[:,1,:]) == diag(T[1,:,:])
true
Also, normsq
can be computed analogously in Tensor
format
t2 = Tensor(2, opq)
normsq == [ t2.get([i, i]) for i in 0:dim(opq)-1]
true
Arbitrary Weights
Of course, the type OrthoPoly
can be constructed for arbitrary weights $w(t)$. In this case we have to compute the orthogonal basis and the respective quadrature rule. Let's re-work the above example by hand.
using SpecialFunctions
supp = (0, 1)
w(t) = (t^(s_α-1)*(1-t)^(s_β-1)/SpecialFunctions.beta(s_α,s_β))
my_meas = Measure("my_meas", w, supp, false)
my_opq = OrthoPoly("my_op", deg, my_meas; Nrec=n, addQuadrature = true)
OrthoPoly{Vector{Float64}, Measure, Quad{Float64, Vector{Float64}}}("my_op", 4, [0.39622641509448137, 0.45308865339922233, 0.47326557666900376, 0.4827290893535496, 0.48792334819610145, 0.4910806427875327, 0.4931429219148361, 0.4945640779987818, 0.495584908427019, 0.4963428640690031, 0.49692106738740605, 0.4973721930564992, 0.4977309380461525, 0.4980209140324845, 0.49825864209832377, 0.49845596313321483, 0.4986215435256319, 0.4987618444763114, 0.4988817627307632, 0.4989850641465735], [0.9999999999996537, 0.037973181440564815, 0.049528193464296905, 0.054496397032274234, 0.05707494398455889, 0.05858224318299423, 0.059538775227572845, 0.06018346809956606, 0.060638471958479156, 0.060971504205455586, 0.0612225580532287, 0.061416485350391216, 0.06156938828922262, 0.06169207139558928, 0.06179200320786561, 0.06187448002171521, 0.06194334216236786, 0.062001429133996595, 0.06205087702968264, 0.06209331771849251], Measure("my_meas", Main.w, (0.0, 1.0), false, Dict{Any, Any}()), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.008788905307716177, 0.02845010617907043, 0.0585212554299101, 0.09833814003385198, 0.14702350696982944, 0.20350504897777927, 0.2665388311160098, 0.33473664359009403, 0.40659656310782666, 0.4805360282406439, 0.5549266957885687, 0.6281303108684106, 0.6985348047827785, 0.7645898366195809, 0.824841025625018, 0.877962211974623, 0.9227853738926778, 0.9583291958881769, 0.9838388820037522], [0.001151677857036659, 0.007007259641384841, 0.0203148927103142, 0.0414025594490785, 0.06777171118707008, 0.09471990539100715, 0.11673118519090996, 0.1291364544780523, 0.12947566139904687, 0.11812235960758868, 0.09800772590075595, 0.07359961540231305, 0.04953526814944095, 0.02938957362654602, 0.0149639646714071, 0.00625271092938539, 0.001979666990405483, 0.00040351473706229986, 3.429268084906047e-5]))
Now we can compute the squared norms $\| \phi_k \|^2$
my_normsq = computeSP2(my_opq)
5-element Vector{Float64}:
0.9999999999996537
0.037973181440551666
0.0018807430768424916
0.00010249372143130948
5.849823409460971e-6
And the tensor
my_t = Tensor(m, my_opq)
my_T = [ my_t.get([i1,i2,i3]) for i1=0:dim(opq)-1,i2=0:dim(opq)-1,i3=0:dim(opq)-1]
5×5×5 Array{Float64, 3}:
[:, :, 1] =
1.0 0.0 0.0 0.0 0.0
0.0 0.0379732 0.0 0.0 0.0
0.0 0.0 0.00188074 0.0 0.0
0.0 0.0 0.0 0.000102494 0.0
0.0 0.0 0.0 0.0 5.84982e-6
[:, :, 2] =
0.0 0.0379732 0.0 0.0 0.0
0.0379732 0.00215924 0.00188074 0.0 0.0
0.0 0.00188074 0.000144891 0.000102494 0.0
0.0 0.0 0.000102494 8.86598e-6 5.84982e-6
0.0 0.0 0.0 5.84982e-6 5.36411e-7
[:, :, 3] =
0.0 0.0 0.00188074 0.0 0.0
0.0 0.00188074 0.000144891 0.000102494 0.0
0.00188074 0.000144891 0.000127149 1.0934e-5 5.84982e-6
0.0 0.000102494 1.0934e-5 7.80614e-6 7.09802e-7
0.0 0.0 5.84982e-6 7.09802e-7 4.73123e-7
[:, :, 4] =
0.0 0.0 0.0 0.000102494 0.0
0.0 0.0 0.000102494 8.86598e-6 5.84982e-6
0.0 0.000102494 1.0934e-5 7.80614e-6 7.09802e-7
0.000102494 8.86598e-6 7.80614e-6 9.40423e-7 5.08904e-7
0.0 5.84982e-6 7.09802e-7 5.08904e-7 6.53232e-8
[:, :, 5] =
0.0 0.0 0.0 0.0 5.84982e-6
0.0 0.0 0.0 5.84982e-6 5.36411e-7
0.0 0.0 5.84982e-6 7.09802e-7 4.73123e-7
0.0 5.84982e-6 7.09802e-7 5.08904e-7 6.53232e-8
5.84982e-6 5.36411e-7 4.73123e-7 6.53232e-8 3.55404e-8
Let's compare the results:
abs.(normsq-my_normsq)
5-element Vector{Float64}:
3.462785613805863e-13
5.5892790395972725e-14
7.394909337654632e-15
8.643530769597563e-16
9.28813978311703e-17
norm(T-my_T)
3.600949810975133e-13
The possibility to create quadrature rules for arbitrary weights should be reserved to cases different from classical ones.
Multivariate Polynomials
For multivariate polynomials the syntax for Tensor
is very much alike, except that we are dealing with the type MultiOrthoPoly
now.
mop = MultiOrthoPoly([opq, my_opq], deg)
MultiOrthoPoly{ProductMeasure, Quad{Float64, Vector{Float64}}, Vector{AbstractOrthoPoly{M, Quad{Float64, Vector{Float64}}} where M<:AbstractMeasure}}(["Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}", "my_op"], 4, 15, [0 0; 1 0; … ; 1 3; 0 4], ProductMeasure(PolyChaos.var"#w#39"{Vector{AbstractOrthoPoly{M, Quad{Float64, Vector{Float64}}} where M<:AbstractMeasure}}(AbstractOrthoPoly{M, Quad{Float64, Vector{Float64}}} where M<:AbstractMeasure[Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}(4, [0.39622641509433965, 0.45308865339881105, 0.4732655766681396, 0.482729089351984, 0.4879233481934926, 0.49108064278342917, 0.49314292190864784, 0.4945640779897571, 0.4955849084142159, 0.4963428640512603, 0.49692106736331404, 0.4973721930243669, 0.49773093800397555, 0.4980209139779084, 0.4982586420286101, 0.4984559630451982, 0.4986215434156863, 0.4987618443402984, 0.4988817625639983, 0.4989850639437675], [1.0, 0.03797318144060756, 0.049528193464418745, 0.054496397032519545, 0.05707494398498378, 0.05858224318366976, 0.05953877522858778, 0.06018346810103027, 0.0606384719605267, 0.06097150420824787, 0.061222558056958205, 0.06141648535528466, 0.061569388295544744, 0.061692071403646374, 0.06179200321800959, 0.061874480034347076, 0.06194334217794196, 0.06200142915302492, 0.06205087705273863, 0.062093317746215744], Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.00878890525374891, 0.028450106116514506, 0.05852125536671723, 0.09833813997192814, 0.14702350691009355, 0.20350504892083368, 0.266538831062297, 0.3347366435399401, 0.40659656306145764, 0.48053602819819125, 0.5549266957500689, 0.6281303108338084, 0.6985348047519305, 0.7645898365922597, 0.8248410256009169, 0.8779622119533655, 0.9227853738738273, 0.9583291958712525, 0.9838388819882817], [0.00115167784765143, 0.007007259624729887, 0.020314892689531894, 0.04140255942833847, 0.06777171117046146, 0.09471990538153915, 0.11673118518989155, 0.12913645448496422, 0.12947566141190417, 0.11812235962361355, 0.09800772591714006, 0.07359961541684255, 0.04953526816082873, 0.029389573634464723, 0.014963964676256315, 0.006252710931950127, 0.0019796669915308576, 0.00040351473743713024, 3.429268092282879e-5])), OrthoPoly{Vector{Float64}, Measure, Quad{Float64, Vector{Float64}}}("my_op", 4, [0.39622641509448137, 0.45308865339922233, 0.47326557666900376, 0.4827290893535496, 0.48792334819610145, 0.4910806427875327, 0.4931429219148361, 0.4945640779987818, 0.495584908427019, 0.4963428640690031, 0.49692106738740605, 0.4973721930564992, 0.4977309380461525, 0.4980209140324845, 0.49825864209832377, 0.49845596313321483, 0.4986215435256319, 0.4987618444763114, 0.4988817627307632, 0.4989850641465735], [0.9999999999996537, 0.037973181440564815, 0.049528193464296905, 0.054496397032274234, 0.05707494398455889, 0.05858224318299423, 0.059538775227572845, 0.06018346809956606, 0.060638471958479156, 0.060971504205455586, 0.0612225580532287, 0.061416485350391216, 0.06156938828922262, 0.06169207139558928, 0.06179200320786561, 0.06187448002171521, 0.06194334216236786, 0.062001429133996595, 0.06205087702968264, 0.06209331771849251], Measure("my_meas", Main.w, (0.0, 1.0), false, Dict{Any, Any}()), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.008788905307716177, 0.02845010617907043, 0.0585212554299101, 0.09833814003385198, 0.14702350696982944, 0.20350504897777927, 0.2665388311160098, 0.33473664359009403, 0.40659656310782666, 0.4805360282406439, 0.5549266957885687, 0.6281303108684106, 0.6985348047827785, 0.7645898366195809, 0.824841025625018, 0.877962211974623, 0.9227853738926778, 0.9583291958881769, 0.9838388820037522], [0.001151677857036659, 0.007007259641384841, 0.0203148927103142, 0.0414025594490785, 0.06777171118707008, 0.09471990539100715, 0.11673118519090996, 0.1291364544780523, 0.12947566139904687, 0.11812235960758868, 0.09800772590075595, 0.07359961540231305, 0.04953526814944095, 0.02938957362654602, 0.0149639646714071, 0.00625271092938539, 0.001979666990405483, 0.00040351473706229986, 3.429268084906047e-5]))]), AbstractMeasure[Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Measure("my_meas", Main.w, (0.0, 1.0), false, Dict{Any, Any}())]), AbstractOrthoPoly{M, Quad{Float64, Vector{Float64}}} where M<:AbstractMeasure[Beta01OrthoPoly{Vector{Float64}, Beta01Measure, Quad{Float64, Vector{Float64}}}(4, [0.39622641509433965, 0.45308865339881105, 0.4732655766681396, 0.482729089351984, 0.4879233481934926, 0.49108064278342917, 0.49314292190864784, 0.4945640779897571, 0.4955849084142159, 0.4963428640512603, 0.49692106736331404, 0.4973721930243669, 0.49773093800397555, 0.4980209139779084, 0.4982586420286101, 0.4984559630451982, 0.4986215434156863, 0.4987618443402984, 0.4988817625639983, 0.4989850639437675], [1.0, 0.03797318144060756, 0.049528193464418745, 0.054496397032519545, 0.05707494398498378, 0.05858224318366976, 0.05953877522858778, 0.06018346810103027, 0.0606384719605267, 0.06097150420824787, 0.061222558056958205, 0.06141648535528466, 0.061569388295544744, 0.061692071403646374, 0.06179200321800959, 0.061874480034347076, 0.06194334217794196, 0.06200142915302492, 0.06205087705273863, 0.062093317746215744], Beta01Measure(PolyChaos.var"#102#103"{Float64, Float64}(2.1, 3.2), (0.0, 1.0), false, 2.1, 3.2), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.00878890525374891, 0.028450106116514506, 0.05852125536671723, 0.09833813997192814, 0.14702350691009355, 0.20350504892083368, 0.266538831062297, 0.3347366435399401, 0.40659656306145764, 0.48053602819819125, 0.5549266957500689, 0.6281303108338084, 0.6985348047519305, 0.7645898365922597, 0.8248410256009169, 0.8779622119533655, 0.9227853738738273, 0.9583291958712525, 0.9838388819882817], [0.00115167784765143, 0.007007259624729887, 0.020314892689531894, 0.04140255942833847, 0.06777171117046146, 0.09471990538153915, 0.11673118518989155, 0.12913645448496422, 0.12947566141190417, 0.11812235962361355, 0.09800772591714006, 0.07359961541684255, 0.04953526816082873, 0.029389573634464723, 0.014963964676256315, 0.006252710931950127, 0.0019796669915308576, 0.00040351473743713024, 3.429268092282879e-5])), OrthoPoly{Vector{Float64}, Measure, Quad{Float64, Vector{Float64}}}("my_op", 4, [0.39622641509448137, 0.45308865339922233, 0.47326557666900376, 0.4827290893535496, 0.48792334819610145, 0.4910806427875327, 0.4931429219148361, 0.4945640779987818, 0.495584908427019, 0.4963428640690031, 0.49692106738740605, 0.4973721930564992, 0.4977309380461525, 0.4980209140324845, 0.49825864209832377, 0.49845596313321483, 0.4986215435256319, 0.4987618444763114, 0.4988817627307632, 0.4989850641465735], [0.9999999999996537, 0.037973181440564815, 0.049528193464296905, 0.054496397032274234, 0.05707494398455889, 0.05858224318299423, 0.059538775227572845, 0.06018346809956606, 0.060638471958479156, 0.060971504205455586, 0.0612225580532287, 0.061416485350391216, 0.06156938828922262, 0.06169207139558928, 0.06179200320786561, 0.06187448002171521, 0.06194334216236786, 0.062001429133996595, 0.06205087702968264, 0.06209331771849251], Measure("my_meas", Main.w, (0.0, 1.0), false, Dict{Any, Any}()), Quad{Float64, Vector{Float64}}("golubwelsch", 19, [0.008788905307716177, 0.02845010617907043, 0.0585212554299101, 0.09833814003385198, 0.14702350696982944, 0.20350504897777927, 0.2665388311160098, 0.33473664359009403, 0.40659656310782666, 0.4805360282406439, 0.5549266957885687, 0.6281303108684106, 0.6985348047827785, 0.7645898366195809, 0.824841025625018, 0.877962211974623, 0.9227853738926778, 0.9583291958881769, 0.9838388820037522], [0.001151677857036659, 0.007007259641384841, 0.0203148927103142, 0.0414025594490785, 0.06777171118707008, 0.09471990539100715, 0.11673118519090996, 0.1291364544780523, 0.12947566139904687, 0.11812235960758868, 0.09800772590075595, 0.07359961540231305, 0.04953526814944095, 0.02938957362654602, 0.0149639646714071, 0.00625271092938539, 0.001979666990405483, 0.00040351473706229986, 3.429268084906047e-5]))])
mt2 = Tensor(2,mop)
mt3 = Tensor(3,mop)
mT2 = [ mt2.get([i,i]) for i=0:dim(mop)-1 ]
15-element Vector{Float64}:
0.9999999999996537
0.03797318144059441
0.037973181440551666
0.0018807430768492354
0.00144196250871918
0.0018807430768424916
0.00010249372143213834
7.141779810028215e-5
7.141779810010645e-5
0.00010249372143130948
5.8498234095518276e-6
3.892012680461295e-6
3.5371945211048704e-6
3.892012680434202e-6
5.849823409460971e-6
Notice that mT2
carries the elements of the 2-dimensional tensors for the univariate bases opq
and my_opq
. The encoding is given by the multi-index mop.ind
mop.ind
15×2 Matrix{Int64}:
0 0
1 0
0 1
2 0
1 1
0 2
3 0
2 1
1 2
0 3
4 0
3 1
2 2
1 3
0 4
To cross-check the results we can distribute the multi-index back to its univariate indices with the help of findUnivariateIndices
.
ind_opq = findUnivariateIndices(1,mop.ind)
ind_my_opq = findUnivariateIndices(2,mop.ind)
5-element Vector{Int64}:
1
3
6
10
15
mT2[ind_opq] - normsq
5-element Vector{Float64}:
-3.462785613805863e-13
-1.3149203947904198e-14
-6.51171824794794e-16
-3.54940686217442e-17
-2.025255776885032e-18
mT2[ind_my_opq] - my_normsq
5-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0