Non-linear system of equations Julia

Question

I'm trying to solve a large number (50) of non-linear simultaneous equations in Julia. For the moment I'm just trying to make this work with 2 equations to get the syntax right etc. However, I've tried a variety of packages/tools - NLsolve, nsolve in SymPy and NLOpt in JuMP (where I ignore the objective function and just enter equality constraints)- without much luck. I guess I should probably focus on making it work in one. I'd appreciate any advice on choice of packages and if possible code.

Here's how I tried to do it in NLsolve (using it in mcpsolve mode so I can impose constraints on the variables I am solving for - x[1] and x[2] - which are unemployment rates and so bounded between zero and 1) :

using Distributions
using Devectorize
using Distances
using StatsBase
using NumericExtensions
using NLsolve

beta = 0.95                                                                
xmin= 0.73                                                                 
xmax = xmin+1                                                             
sigma = 0.023                                                            
eta = 0.3                                        
delta = 0.01                                                                                                
gamma=0.5                                                                   
kappa = 1                                                                  
psi=0.5
ns=50
prod=linspace(xmin,xmax,ns)
l1=0.7
l2=0.3                                            
wbar=1
r=((1/beta)-1-1e-6 +delta)


## Test code

function f!(x, fvec)

    ps1= wbar + (kappa*(1-beta*(1-sigma*((1-x[1])/x[1]))))
    ps2= wbar + (kappa*(1-beta*(1-sigma*((1-x[2])/x[2]))))

    prod1=prod[1]
    prod2=prod[50]
    y1=(1-x[1])*l1
    y2=(1-x[2])*l2
    M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))
    K=((r/eta)^(1/(eta-1)))*M

    pd1=(1-eta)*(K^eta)*(M^(-eta))*prod1
    pd2=(1-eta)*(K^eta)*(M^(-eta))*prod2

    fvec[1]=pd1-ps1
    fvec[2]=pd2-ps2
end

mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3, 0.3])

I get this error message:

Any suggestions are very welcome! I appreciate the formulas are pretty ugly so let me know if any further simplifications helpful (I have tried!).

Answer 1

I've thought you were giving initial conditions outside the range of bounds, for I tried mcpsolve(f!,[0.0,0.0],[0.0,0.0],[0.3, 0.3]) and it worked.

However, I've also tried other combinations:

mcpsolve(f!,[0.4,0.4], [0.0,0.0], [0.3, 0.3]) did worked

mcpsolve(f!,[0.4,0.4], [0.3,0.3], [1.0,1.0]) did not

mcpsolve(f!,[0.6,0.6], [1.0,1.0], [0.3,0.3]) did not

Have you checked these values on your test?