# Multidimensional integration and differentiation

**URL:** <https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163>\
**Category:** Numerics\
**Created:** [February 27, 2021, 10:40pm UTC](https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163 "2021-02-27T22:40:06Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![simonf](https://avatars.discourse-cdn.com/v4/letter/s/f08c70/32.png) [@simonf](https://discourse.julialang.org/u/simonf)\
**Post date:** [February 27, 2021, 10:40pm UTC](https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163/1 "2021-02-27T22:40:06Z")

</div>

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](https://github.com/JuliaMath/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](https://github.com/JuliaIntervals/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!

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

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 27, 2021, 11:04pm UTC](https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163/2 "2021-02-27T23:04:58Z")

</div>

> [@simonf](#):
>
> I’m wondering if I can implement this with numerical packages. I have tried with [HCubature.jl](https://github.com/JuliaMath/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.

Perform a change of variables to a rectangular domain, as explained [here](https://discourse.julialang.org/t/2d-integration-over-non-rectangular-domain-using-cubature/2991) and [here](https://discourse.julialang.org/t/nested-quadgk-long-compilation-time-gets-short-when-plotting-the-argument-beforehand-but-not-when-calling-it/24519/3) and [here](https://github.com/stevengj/HCubature.jl/issues/16), for example.

---

<div class="post-metadata">

**Author:** ![simonf](https://avatars.discourse-cdn.com/v4/letter/s/f08c70/32.png) [@simonf](https://discourse.julialang.org/u/simonf)\
**Post date:** [February 27, 2021, 11:15pm UTC](https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163/3 "2021-02-27T23:15:49Z")

</div>

Thanks, I’ll try that!

---

<div class="post-metadata">

**Author:** ![simonf](https://avatars.discourse-cdn.com/v4/letter/s/f08c70/32.png) [@simonf](https://discourse.julialang.org/u/simonf)\
**Post date:** [March 10, 2021, 10:42pm UTC](https://discourse.julialang.org/t/multidimensional-integration-and-differentiation/56163/4 "2021-03-10T22:42:11Z")

</div>

I’ve made progress with the integration and this works fine on its own - thanks! However, when differentiating with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) and finding roots with [IntervalRootFinding.jl](https://github.com/JuliaIntervals/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.

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

```