Multidimensional integration and differentiation

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))
``````

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

3 Likes

Thanks, I’ll try that!

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)

``````