Maybe you could use complex numbers, so the domain error isn't thrown, then take the real part at the end (after checking that im==0)? Alternatively, you could file an issue at NLsolve stating that the objective function is evaluated outside the bounds. Maybe it is possible improve the algorithm, maybe not.
(Also, to get better feedback and be kinder to the reviewers, try to make a smaller example which does only uses the minimal number of dependencies. And make sure it works as, I think, your example does not work, unless `ns` is something which is exported by one of the used packages I don't have installed.) On Tue, 2015-06-30 at 07:01, David Zentler-Munro <[email protected]> wrote: > '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 > prod=linspace(xmin,xmax,ns) > l1=0.7 > l2=0.3 > assets = linspace(0,amaximum,na); > > 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]))))/ > (beta*((x[1]/(sigma*(1-x[1])))^(gamma/(1-gamma)))*(1/(2-x[1])))) > > ps2= wbar + ((kappa*(1-beta*(1-sigma*((1-x[2])/x[1]))))/ > (beta*((x[2]/(sigma*(1-x[2])))^(gamma/(1-gamma)))*(1/(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)))^(psi/(psi-1)) > K=((r/eta)^(1/(eta-1)))*M > > pd1=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+ > ((prod2*y2)^((psi- 1)/psi)))^(1/(psi-1)))* > ((prod1*y1)^(-1/psi))*prod1 > > pd2=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+ > ((prod2*y2)^((psi-1)/psi)))^(1/(psi-1)))* > ((prod2*y2)^(-1/psi))*prod2 > > fvec[1]=pd1-ps1 > fvec[2]=pd2-ps2 > end > > mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3, 0.3]) > > > However, I get this error > > > [image: error message] > > As I understand it, this occurs because I am trying to take the root of a > negative number. However, the bounds on x[1] and x[2] should stop this > happening. Any advice very welcome. I appreciate the formulas are pretty > ugly so let me know if any further simplifications helpful (I have tried!) > > David
