SciML / PolyChaos.jl

A Julia package to construct orthogonal polynomials, their quadrature rules, and use it with polynomial chaos expansions.
https://docs.sciml.ai/PolyChaos/stable/
MIT License
115 stars 26 forks source link

computeTensorizedSP #113

Open vavrines opened 1 year ago

vavrines commented 1 year ago

The computeTensorizedSP seems to drop entries less than 1E-10, This is not very reasonable when using Float64 and can cause weird behaviors when making use of tensor products.

MWE

using PolyChaos

op = Uniform_11OrthoPoly(20, Nrec = 40)
t2 = Tensor(2, op)

l1 = [t2.get([i, i]) for i=0:op.deg]
l2 = computeSP2(op)

l1 is not equal to l2

julia> l1
21-element Vector{Float64}:
 ⋮
 3.6023076661943104e-10
 0.0
 0.0
 0.0
 0.0
julia> l2
21-element Vector{Float64}:
 ⋮
 3.6023076661943104e-10
 9.013566368226456e-11
 2.2551316627840712e-11
 5.64173617647297e-12
 1.4113161166911746e-12

computeTensorizedSP is the default method used to construct Tensor, see https://github.com/SciML/PolyChaos.jl/blob/e20be072e31754f1f7627628c65b43b4111b6455/src/typesTensor.jl#L15 In the following example, Tensor.get() would not work at high orders as the denominators are zero.

function ODEgalerkin(du,u,p,t)
   du[:] = [ sum( p[j+1]*u[k+1]*t3.get([j,k,m])/t2.get([m,m]) for j=0:L for k=0:L) for m=0:L ] 
end