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

Multivariate integration #48

Open Alexander-Murray opened 3 years ago

Alexander-Murray commented 3 years ago

Hey, I tried following the tutorial examples here but I can't seem to figure out the right syntax for multivariate integration. Here is what I am trying:

deg = [3, 5, 6, 4]
d = minimum(deg)
op1 = GaussOrthoPoly(deg[1], addQuadrature=true);
op2 = Uniform01OrthoPoly(deg[2], addQuadrature=true);
op3 = Beta01OrthoPoly(deg[3], 2, 1.2, addQuadrature=true);
ops = [op1, op2, op3];
mop = MultiOrthoPoly(ops, d);
integrate((a,b,c,d) -> a*b*c*d, mop)

however, I get the following error:

type MultiOrthoPoly has no field quad

Stacktrace:
 [1] getproperty(::MultiOrthoPoly{ProductMeasure,Quad{Float64,Array{Float64,1}},Array{AbstractCanonicalOrthoPoly{Array{Float64,1},M,Quad{Float64,Array{Float64,1}}} where M,1}}, ::Symbol) at .\Base.jl:33
 [2] integrate(::Function, ::MultiOrthoPoly{ProductMeasure,Quad{Float64,Array{Float64,1}},Array{AbstractCanonicalOrthoPoly{Array{Float64,1},M,Quad{Float64,Array{Float64,1}}} where M,1}}) at C:\Users\mrra\.julia\packages\PolyChaos\AlY80\src\auxfuns.jl:96
 [3] top-level scope at In[5]:8
 [4] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

Is there any way of doing this with the PolyChaos package?

Alexander-Murray commented 3 years ago

One solution is to use multiple univariate integrals:

deg = [3, 5, 6, 4]
op1 = GaussOrthoPoly(deg[1], addQuadrature=true);
op2 = Uniform01OrthoPoly(deg[2], addQuadrature=true);
op3 = Beta01OrthoPoly(deg[3], 2, 1.2, addQuadrature=true);
ops = [op1, op2, op3];
integrate(c->integrate(b->integrate((a)->a*b*c, op1),op2),op3)

however, I'd be interested to hear about other approaches.

timueh commented 3 years ago

hi Alex, yes, multiple univariate integrals is the way to proceed.

Put differently, there is not support for generic multivariate integration, only integration of product measures for which you can decompose the $n$-variate integral to $n$ univariate integrals. For each of the $n$ univariate integrals you can then use the quad rules.