[R] newton method for single nonlinear equation

[Berend Hasselman]

Roslina Zakaria wrote:
> 
> newton.inputsingle <- function(pars,n)
> {  runi    <- runif(974, min=0, max=1)
>    lendt   <- length(runi)
>    ## Parameter to estimate
>    z <- vector(length=lendt, mode= "numeric")
>    z  <- pars[1]
>    
>    ## Constant value  
>    
>    alp  <- 2.0165 ; rho <- 0.868; 
>    c    <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)
>    
>    for (i in 1:n)
>    {  t1   <- exp(-pars[1]/(1-rho))                       
>       t2   <- (pars[1]*(1-rho)/(2*sqrt(rho)))^(alp-0.5)   
>       bes1 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-0.5)  
>       bes2 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp-1.5)
>       bes3 <- besselI(pars[1]*sqrt(rho)/(1-rho),alp+0.5)
>  
>    ## Equation
>       f   <- c*t1*t2*bes1 - runi
>    
>    ## derivative
>       fprime   <- c*t1*t2*( -bes1/(1-rho) + (alp-0.5)*bes1/pars[1] +
> sqrt(rho)*(bes2-bes3)/(2*(1-rho)))
>       z[i+1] <- z[i] - f/fprime 
>       }
>       z
> }
>  
> pars <- 0.5          
> newton.inputsingle(pars,5)
>  
> The output :
>  
>> pars <- 0.5          
>> newton.inputsingle(pars,5)
> [1]  0.5000000 -0.4826946 -1.4653892 -2.4480838 -3.4307784 -4.4134730
> Warning messages:
> 1: In z[i + 1] <- z[i] - f/fprime :
>   number of items to replace is not a multiple of replacement length
> 

The warning message is the result of assigning a vector (f/fprime) to a
scalar.
You are also redefining an R built-in function: c

Furthermore it is unclear what you are trying to do.

I would advise: Keep It Simple and Stupid (KISS).
To help you along I have lifted the relevant function out of your
newton.inputsingle
(or at least what I think your function is)

zztest <- function(x) {
   ## Constant value
    alp  <- 2.0165 ; rho <- 0.868;
    czz  <- sqrt(pi)/(gamma(alp)*(1-rho)^alp)

    t1   <- exp(-x/(1-rho))
    t2   <- (x*(1-rho)/(2*sqrt(rho)))^(alp-0.5)
    bes1 <- besselI(x*sqrt(rho)/(1-rho),alp-0.5)
    bes2 <- besselI(x*sqrt(rho)/(1-rho),alp-1.5)
    bes3 <- besselI(x*sqrt(rho)/(1-rho),alp+0.5)

   ## Equation
    f   <- czz*t1*t2*bes1 - .1
}

pars <- 0.5
zztest(pars)
plot(zztest,0,10)
abline(0,0,lty=2)
uniroot(zztest, c(0,2))
uniroot(zztest, c(4,8))

The plot indicates how your function behaves.
The two uniroot calls solve for the two roots that appear in the plot.
You should be able to proceed on your own from here.
BTW: I do hope that this not homework.

good luck

Berend

-- 
View this message in context: http://n4.nabble.com/newton-method-for-single-nonlinear-equation-tp1289991p1310861.html
Sent from the R help mailing list archive at Nabble.com.