I have two bivariate functions (more generally, n n-variable functions) f_1(x_1,x_2) and f_2(x_1,x_2). I need to solve the following system of equations: \frac{\partial f_1}{\partial x_1} = 0 and \frac{\partial f_2}{\partial x_2} = 0 (the context here is a Cournot-type duopoly model). The functions are too complex to differentiate by hand, so I would like to use AD tools to solve this. However, I cannot come up with a good way to grab the partial derivatives and use them as expressions inside tools such as NLsolve or JuMP. What would be the correct way to do this?
My current approach is to calculate the Jacobian matrix, and point to the corresponding entries when solving the system. This gives me the correct solution, but feels more like a clumsy workaround than an actual solution, in particular for larger problems.
using ForwardDiff
using NLsolve
f1(x1,x2) = x1*(50-2*(x1+x2))-10-2x1
f2(x1,x2) = x2*(50-2*(x1+x2))-10-2x2
function jac(f1,f2,a)
F(x1,x2) = [f1(x1,x2), f2(x1,x2)]
return ForwardDiff.jacobian(x -> F(x[1],x[2]),a)
end
function f!(F,x)
F[1] = jac(f1,f2,x)[1,1] # df1/dx1 == 0
F[2] = jac(f1,f2,x)[2,2] # df2/dx2 == 0
end
x0 = [5.,5.] # Initial guess
res = nlsolve(f!, x0)
println(res.zero) # Correct solution: [8.0,8.0]
The example above is from: http://faculty.econ.ucdavis.edu/faculty/bonanno/teaching/200C/Cournot.pdf