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]]
______________________________________________
R-help@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.
