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