Roslina Zakaria wrote:
Hi I would like to solve a system of nonlinear equations below using
dfsane function
mn - 142.36; vr - 9335.69 ; sk - 0.81; kur - 0.25
test_fn - function(p)
{ f - rep(NA, length(p))
f[1] - p[1]*(p [2]+p[3])- mn
f[2] - - vr + 2*p[1]*p[2]*p[3]*(p[4]-1)+p[1]*(p[2]+p[3])^2
f[3] - - sk + (p[1]*(p
[2]+p[3])^3*(p[1]+1)*(p[1]+2)-6*p[1]*p[2]*p[3]*(p[2]+ p[3])*(1- p[4])*( p
[1]+1))/(p[1]*(p[2]+p[3])^2*(p[2]+p[3])-2*p[1]*p[2]*p[3]*(1-p[4]))^1.5
f[4] - - kur + (p[1]*(p[2]+p[3])^4*(p[1]+3)*(p[1]+2)*(p[1]+1)
-12*p[1]*p[2]*p[3]*(1-p[4])*(p[2]+p[3])^2*(p[1]+2)*(p[1]+1)+12*p[1]*p[2]^2*p[3]^2*(1-p[4])^2*(p[1]+1))/(p[1]*(p[2]+p[3])^2*(p[2]+p[3])-2*p[1]*p[2]*p[3]*(1-p[4]))^2
f
}
p0 - c(1.3, 50,60,0.8)
dfsane(par=p0, fn = test_fn, control=list(trace=FALSE))
I tried solving your system of equations with package nleqslv.
It reports an ill conditioned/singular jacobian and cannot solve the system.
I did a rough calculation of the jacobian in the initial point you provided.
n - length(p0)
J - matrix(data=rep(0,n*n),nrow=n)
eps - 1e-8
f0 - test_fn(p0)
J[,1] - (test_fn(p0+c(eps,0 ,0 ,0 )) - f0)/eps
J[,2] - (test_fn(p0+c(0 ,eps,0 ,0 )) - f0)/eps
J[,3] - (test_fn(p0+c(0 ,0 ,eps,0 )) - f0)/eps
J[,4] - (test_fn(p0+c(0 ,0 ,0 ,eps)) - f0)/eps
J
[,1] [,2] [,3] [,4]
[1,] 1.10e+02 1.30e+00 1.30e+00 0.00e+00
[2,] 1.09e+04 2.548002e+02 2.63e+02 7.80e+03
[3,] 2.014600e-03 -7.268630e-05 -7.097656e-05 2.569023e-03
[4,] 4.256762e-04 -3.317069e-05 -3.225198e-05 1.378761e-03
Columns 2 and 3 are almost identical.
p[2] and p[3] seem to be interchangeable.
If you do svd(J) then you will see the singular values
$d
[1] 1.340860e+04 6.401512e+01 2.129760e-04 5.305005e-10
smallest singular value/largest singular value is very close to machine
precision.
Your system of equations is almost singular.
Berend
--
View this message in context:
http://r.789695.n4.nabble.com/BB-package-tp2074937p2076045.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
