Alicia Ellis
2012-Dec-04 16:06 UTC
[R] Solve system of equations (nleqslv) only returns origin
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?
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)
}
Xstart=c(10000, 200, 0.5, 0.5, 12)
fstart= firstordercond(Xstart)
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]]
Berend Hasselman
2012-Dec-04 17:50 UTC
[R] Solve system of equations (nleqslv) only returns origin
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.