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
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
op = OrthoPoly("beta01",deg,Dict(:shape_a=>s_α,:shape_b=>s_β);Nrec=n)
To add the corresponding quadrature rule there is the composite struct OrthoPolyQ
whose simplest constructor reads
opq = OrthoPolyQ(op,n)
By default, an $n$-point Gauss quadrature rule is create relative to the underlying measure op.meas
, where $n$ is the number of recurrence coefficients stored in op.α
and op.β
. The type OrthoPolyQ
has just two fields: an OrthoPoly
, and a Quad
.
To compute the squared norms
of the basis we call computeSP2()
normsq = computeSP2(opq)
For the general case
there exists a type Tensor
that requires only two arguments: the dimension$m \geq 1$, and an OrthoPolyQ
m = 3
t = Tensor(3,opq)
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])
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]
Notice that we can cross-check the results.
using LinearAlgebra
@show normsq == LinearAlgebra.diag(T[:,:,1])
@show normsq == LinearAlgebra.diag(T[:,1,:])
@show normsq == LinearAlgebra.diag(T[1,:,:])
Also, normsq
can be computed analogously in Tensor
format
t2 = Tensor(2,opq)
@show normsq == [ t2.get([i,i]) for i=0:dim(opq)-1]
Arbitrary Weights
Of course, the type OrthoPolyQ
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)
function w(t)
supp[1]<=t<=supp[2] ? (t^(s_α-1)*(1-t)^(s_β-1)/SpecialFunctions.beta(s_α,s_β)) : error("$t not in support")
end
my_meas = Measure("my_meas",w,supp,false,Dict())
my_op = OrthoPoly("my_op",deg,my_meas;Nrec=n)
my_quad = Quad(n,my_op)
my_opq = OrthoPolyQ(my_op,my_quad)
Now we can compute the squared norms $\| \phi_k \|^2$
my_normsq = computeSP2(my_opq)
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]
Let's compare the results:
@show abs.(normsq-my_normsq)
@show norm(T-my_T)
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)
mt2 = Tensor(2,mop)
mt3 = Tensor(3,mop)
mT2 = [ mt2.get([i,i]) for i=0:dim(mop)-1 ]
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
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)
@show mT2[ind_opq] - normsq
@show mT2[ind_my_opq] - my_normsq;