[R] newton.method

Ravi Varadhan rvaradhan at jhmi.edu
Mon Aug 2 15:17:53 CEST 2010


Try this:

for( x in seq(0,2,length=10)) { z <- nleqslv(x,f); print(c(z$x,z$fvec))}

You'll see that it identifies 0.7347242 correctly as a root.  Your function
goes to zero asymptotically as x gets larger.  From numerical algorithms's
point of view, it is essentially 0 for x > 3.  Thus, you have infinite
number of roots larger than 3 or so, and for most starting values Broyden's
method finds those roots.

Moral of the story:  Do due diligence before complaining about algorithms.
You should at least plot the function over a range of x-values, and then add
an `abline(h=0)' to it so you can see what I am talking about.

Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of sammyny
Sent: Monday, August 02, 2010 4:49 AM
To: r-help at r-project.org
Subject: Re: [R] newton.method


I have this function:
function(x)
-0.3*x*exp(-(((log(x)+(0.03+0.3*0.3/2)*0.5)/(0.3*sqrt(0.5)))^2)/2)/(2*sqrt(2
*pi*0.5))
+ 0.03*exp(-0.03*0.5)*pnorm(-(log(x)+(0.03-0.3*0.3/2)*0.5)/(0.3*sqrt(0.5)))

uniroot is giving the correct results.
> uniroot(f,c(0,10))
$root
[1] 0.7347249
$f.root
[1] -1.955740e-07
$iter
[1] 12
$estim.prec
[1] 6.103516e-05

The newton method is not generating good result.
> for( x in seq(0,8,1)) { z <- nleqslv(x,f); print(c(z$x,z$fvec))} 
[1]           NaN 1.340781e+154
[1]           NaN 1.340781e+154
[1]  3.406704e+00 -5.608658e-09
[1]  3.340717e+00 -9.480243e-09
[1]  4.000000e+00 -5.470677e-11
[1]  5.000000e+00 -3.386274e-14
[1]  6.000000e+00 -3.559448e-17
[1]  7.000000e+00 -6.065006e-20
[1]  8.000000e+00 -1.581856e-22


In this case I am looking for largest value of 'x' such that f(x) >= 0

-- 
View this message in context:
http://r.789695.n4.nabble.com/newton-method-tp2306111p2310041.html
Sent from the R help mailing list archive at Nabble.com.

______________________________________________
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.



More information about the R-help mailing list