Multidimensional integration and differentiation

Specific Domains Numerics
1

Hi,

I have a series of parameterized functions containing double integrals and function limits and want to find stationary points of those functions. I have a solution using SymPy.jl, but this is only working well for k=1,2,3. For k = 4, for example, this is extremely slow (I run this for various k, vA, vB, and f(k,x)). MWE is below.

I’m wondering if I can implement this with numerical packages. I have tried with HCubature.jl but can’t get it to work, I think because the upper limit of the inner integral contains the variable of integration of the outer integral.

For differentiation I was thinking to use ForwardDiff.jl and IntervalRootFinding.jl to solve the system of first order conditions.

Any pointers how to use/connect these numerical packages would be much appreciated - or perhaps there are better solutions even. Thank you!

using SymPy

function func(vA,vB,wA::Sym,wB::Sym,a::Sym,b::Sym,f_a::Sym,f_b::Sym,vbar)
   return integrate(integrate(f_a*f_b,(b,wB-wA+a,vbar)),(a,0,wA))*(vA-wA) +
          (integrate(integrate(f_a*f_b,(b,0,wB)),(a,wA,vbar)) +
          integrate(integrate(f_a*f_b,(b,0,wB-wA+a)),(a,0,wA)))*(vB-wB)
end

function findS(vA,vB,a::Sym,b::Sym,f_a::Sym,f_b::Sym,vbar,obj_func::Function)
    SymPy.@vars wA wB
    w = [wA, wB]
    return solve(diff.(obj_func(vA,vB,w...,a,b,f_a,f_b,vbar),w), w)
end

SymPy.@vars x a b wA wB
k = 4
f(k,x) = k*x^(k-1)
println(findS(0.4,0.7,a,b,f(k,a),f(k,b),1.0,func))
2

Perform a change of variables to a rectangular domain, as explained here and here and here, for example.

3 Likes
3

Thanks, I’ll try that!

4

I’ve made progress with the integration and this works fine on its own - thanks! However, when differentiating with ForwardDiff.jl and finding roots with IntervalRootFinding.jl I run into trouble. I’m posting two attempts following the examples I found - unfortunately neither of them works yet. Any help would be much appreciated.

using IntervalArithmetic, IntervalRootFinding, ForwardDiff, HCubature

const X = 0..1

function f(x::Float64, k::Int64)
    return k*x^(k-1)
end

function func(wA, wB, vA, vB, f::Function, k::Int64, vbar::Float64)
    return hcubature(x->f(x[1],k) * f(x[2]*vbar + (1 - x[2])*(wB - wA + x[1]), k)
                    * (vbar - wB + wA - x[1]), [0, 0], [wA, 1])[1] * (vA-wA) + 
           (hcubature(x->f(x[1],k)*f(x[2],k), [wA, 0], [vbar, wB])[1] + 
           hcubature(x->f(x[1],k)*f(x[2]*(wB-wA+x[1]),k) * (wB-wA+x[1]), [0, 0], [wA, 1])[1]) * (vB-wB)
end

function findS(vA::Float64,vB::Float64,d::Function,k::Int64,vbar::Float64,objfn::Function)
	g( (wA,wB) ) = objfn(wA, wB, vA, vB, f, k, vbar)
	return IntervalRootFinding.roots(∇(g), X × X)
end

function findS2(vA::Float64,vB::Float64,d::Function,k::Int64,vbar::Float64,objfn::Function)
	∇g = t -> ForwardDiff.gradient(w -> objfn(w[1], w[2], vA, vB, f, k, vbar), [t[1],t[2]])
    return IntervalRootFinding.roots(∇g, X × X)
end

findS(0.8,0.5,f,1,1.0,func)
findS2(0.8,0.5,f,1,1.0,func)