You could have your function return NaN when one or more arguments are out-of- bounds.
If you want to find a solution very near such a boundary, your performance will presumably stink. I'm not sure what the state of the art here is (for minimization, I'd suggest an interior-point method). Best, --Tim On Tuesday, June 30, 2015 06:18:09 AM David Zentler-Munro wrote: > Thanks Mauro. Fair point about complex code, and yes "ns" was missing. I've > tried to simplify code as follows below (still get same error). As far as I > can see it does seem like NLsolve is evaluating outside the bounds so, > unless anyone has other suggestions, I will file an issue. Here's the > updated code > > 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]) > > On Tuesday, June 30, 2015 at 1:32:03 PM UTC+1, Mauro wrote: > > 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] <javascript:>> 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
