[R] Solve system of equations (nleqslv) only returns origin

[Berend Hasselman]

On 04-12-2012, at 17:06, Alicia Ellis wrote:

> I'm solving 4 complex equations simultaneously.  Code is below.  The code
> returns only zero's for the solution though there should also be a non-zero
> result.  I'm pretty confident that the equations are correct because they
> are straight from a published paper and I checked them pretty thoroughly.
> The parameter values I used are from the published paper as well.  Any
> suggestions for how I can find the maximum non-zero solution instead of
> going straight to the minimum?
> 

Are you trying to minimize a function by solving the first order condition?
If yes then you should try functions such optim, nlminb and there are many more.

See below for more comments.

> Thanks!
> Alicia
> 
> install.packages("nleqslv",
> lib="Users/Alicia/Documents/R.Codes/R.Packages/")
> 
> require(nleqslv)
> 
> ###### Global Parameters ############
> beeta=0.8
> pq=10000
> L=12600
> theta=0.6
> psale=0.6
> mu=psale*(1-theta)
> alphah=0.15
> Cg=6240
> Cs=2820
> A= 100
> D=0.0001
> greekp=0.43
> K=100000
> 
> ##### Species Parameters ##########
> 
> b1=0.38
> p1=16654
> v1 = 0.28
> N1=6000
> g1=1
> delta1=1
> 
> b2=0.4
> p2=2797
> v2 = 0.31
> N2=10000
> g2=1
> delta2=1
> 
> ### Define functions with vector x = c(Lg, Ls, gamma1, gamma2, lamda)
> 
> firstordercond <- function (x) {
> 
>    y=numeric(4)
> 
> 
> y[1]=(alphah/x[3])-(x[5]*((p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])
> +
>      delta1*theta*(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2])))
> 
> 
> y[2]=(alphah/x[4])-(x[5]*((p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])
> +
>      delta2*theta*(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])))
> 
> 
> y[3]=((alphah*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
> +
> ((alphah*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2]))
> +
>      x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
> (b1*(1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*((N1/A)*g1^(greekp))*x[1]^(b1-1))
> +
> 
> (b2*(1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*((N2/A)*g2^(greekp))*x[1]^(b2-1))
> -
>    (delta1*theta*x[3]*((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)) -
> (delta2*theta*x[4]*((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))-Cg*(((N1/A)*g1^(greekp))*b1*x[1]^(b1-1)+((N2/A)*g2^(greekp))*b2*x[1]^(b2-1)))
> 
> 
> y[4]=((alphah*(2*v1*N1*D))/(((N1/A)*g1^(greekp))*x[1]^b1+(2*v1*N1*D)*x[2]))
> + ((alphah*(2*v2*N2*D))/(((N2/A)*g2^(greekp))*x[1]^b2+(2*v2*N2*D)*x[2])) +
>      x[5]*(((-1*beeta*pq*(L-x[1]-x[2])^(beeta-1)) +
> ((1-x[3])*(p1-(((theta+mu)*(((N1/A)*g1^greekp*x[1]^b1)+K))+((theta+mu)*(((1-exp(-2*D*v1*N1))*x[2])+K))))*(2*v1*N1*D))
> +
> 
> ((1-x[4])*(p2-(((theta+mu)*(((N2/A)*g2^greekp*x[1]^b2)+K))+((theta+mu)*(((1-exp(-2*D*v2*N2))*x[2])+K))))*(2*v2*N2*D))
> - (delta1*theta*x[3]*(2*v1*N1*D)) -
>      (delta2*theta*x[4]*(2*v2*N2*D)))-Cs*((2*v1*N1*D)+(2*v2*N2*D)))
> 
>      result=(x)
> 

what is this statement supposed to do?
You are actually returning the input.
And why like that?
Just x or return(x) is quite sufficient.

You should return the vector y i.e. function values.
But y has length 4 and x has length 4.

So where is the fifth value for y?

>    }
> 
> Xstart=c(10000, 200, 0.5, 0.5, 12)
> fstart= firstordercond(Xstart)
> 

If you print fstart you will see that it is identical to Xstart.

You need to rethink you firstordercond() function.
It's not correct.

Berend

> nleqslv(Xstart, firstordercond)
> 
> 
> 
> -- 
> Alicia Ellis
> Postdoc
> Gund Institute for Ecological Economics
> University of Vermont
> 617 Main Street
> Burlington, VT  05405
> (802) 656-1046
> http://www.uvm.edu/~aellis5/
> <http://entomology.ucdavis.edu/faculty/scott/aellis/>
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.