roots-vig.R
In rootoned: Roots of one-dimensional functions in R-only code

### R code from vignette source 'roots-vig.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: chunk1
###################################################
rm(list=ls())
railss<-function(xx){ # SPMA rail problem to provide sumsquares of residuals
    r<-xx[1]
    p<-xx[2] # get the two parameters
    e1<-r*sin(p)-2640
    e2<-r*p-2640.5
    ss<-e1*e1+e2*e2
}
x<-c(10000, 0.5) # start
cat("Function at start = ",railss(x),"\n")
ans<-optim(x, railss)
print(ans)
rr<-ans$par[1]
pp<-ans$par[2]
GG<-4 + 8.5/12
AA<-  rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")


###################################################
### code chunk number 2: chunk2
###################################################
xx<-ans$par # new start
ans2<-optim(xx, railss)
print(ans2)
rr<-ans2$par[1]
pp<-ans2$par[2]
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")


###################################################
### code chunk number 3: chunk3
###################################################
xx<-ans$par # new start
ans2<-optim(xx, railss, control=list(maxit=20000))
print(ans2)
rr<-ans2$par[1]
pp<-ans2$par[2]
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")


###################################################
### code chunk number 4: chunk4
###################################################
xx<-ans$par # new start
ans3<-optim(xx, railss, method='BFGS', control=list(maxit=20000))
print(ans3)
rr<-ans3$par[1]
pp<-ans3$par[2]
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")


###################################################
### code chunk number 5: chunk5
###################################################
xy<-c(70000, .04) # new start
ans4<-optim(xy, railss, control=list(maxit=20000))
print(ans4)
rr<-ans4$par[1]
pp<-ans4$par[2]
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")

cat("\n Try from last solution")
ans5<-optim(c(rr, pp), railss, control=list(maxit=20000))
print(ans5)
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")


###################################################
### code chunk number 6: chunk6
###################################################
tint<-c(0.01, 0.1) # interval for root
g<-function(p) {
 fun<-2640.5 * sin(p)/p - 2640
}
tryp<- uniroot(g, tint)
print(tryp)
pp<-tryp$root
rr <- 2640.5/pp
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")
e1<-rr*sin(pp)-2640
e2<-rr*pp-2640.5
ss<-e1*e1+e2*e2
cat("ss loss function = ", ss, "\n")



###################################################
### code chunk number 7: chunk7
###################################################
tinr<-c(50000, 90000) # interval for root
f<-function(r) {
 fun<-r * sin(2640.5/r) - 2640
}
tryr<- uniroot(f, tinr)
print(tryr)
rr<-tryr$root
pp <- 2640.5/rr
AA<-rr*(1-cos(pp)) 
cat("rail distance = ",AA," + G =",AA+GG,"\n")
e1<-rr*sin(pp)-2640
e2<-rr*pp-2640.5
ss<-e1*e1+e2*e2
cat("ss loss function = ", ss, "\n")



###################################################
### code chunk number 8: chunk8
###################################################
cat("Tangent function false root\n")
mytan<-function(xdeg){ # tangent in degrees
    xrad<-xdeg*pi/180.0 # conversion to radians
    tt<-tan(xrad)
}

cat("find root between 80 and 100 degrees\n")
tint<-c(80,100)
rtan1<-uniroot(mytan, tint)
print(rtan1)
cat("\n\nTighten the tolerance\n")
rtan2<-uniroot(mytan, tint, tol=.Machine$double.eps)
print(rtan2)


###################################################
### code chunk number 9: chunk9
###################################################
cat("Gaussian derivative\n")
der<-function(x,mu,sigma){
    dd<-(-1)*(1/sqrt(2 * pi * sigma^2)) * (exp(-(x - mu)^2/(2 * sigma^2)) * 
    ((x - mu)/(sigma^2)))
}
r1<-uniroot(der, lower=2, upper=6, mu=4, sigma=1)
r.1<-uniroot(der, lower=2, upper=6, mu=4, sigma=.1)
r.01<-uniroot(der, lower=2, upper=6, mu=4, sigma=.01)
r.001<-uniroot(der, lower=2, upper=6, mu=4, sigma=.001)
sig<-c(1, .1, .01, .001)
roo<-c(r1$root, r.1$root, r.01$root, r.001$root)
tabl<-data.frame(sig, roo)
print(tabl)


###################################################
### code chunk number 10: chunk10
###################################################
z <- c(10, -3, 0, 1)
simpol<-function(x){ # calculate polynomial z at x 
   zz <- c(10, -3, 0, 1)
   ndeg<-length(zz)-1 # degree of polynomial
   val<-zz[ndeg+1]
   for (i in 1:ndeg){
       val<-val*x+zz[ndeg+1-i]
   }
   val
}
tint<-c(-5,5)
cat("roots of polynomial specified by ")
print(z)
require(polynom)
allroots<-polyroot(z)
print(allroots)
cat("single root from uniroot on interval ",tint[1],",",tint[2],"\n")
rt1<-uniroot(simpol, tint)
print(rt1)


###################################################
### code chunk number 11: chunk11
###################################################
z <- c(16, -8, 1)
simpol2<-function(x){ # calculate polynomial z at x 
   val <- (x-4)^2
}
tint<-c(-5,5)
cat("roots of polynomial specified by ")
print(z)
library(polynom)
allroots<-polyroot(z)
print(allroots)
cat("single root from uniroot on interval ",tint[1],",",tint[2],"\n")
rt1<-try(uniroot(simpol2, tint))
print(rt1)
cat("\n Try a cubic\n")
cub<-function(z) {val<-(z-4)^3}
cc<-c(-64, 48, -12, 1)
croot<-polyroot(cc)
croot
ans<-uniroot(cub, lower=2, upper=6)
ans


###################################################
### code chunk number 12: chunk12
###################################################
mrate<-function(R){
   val<-0
   den<-1
   fact<-1/6
   term<-1
   rr<-R/200
   repeat { # main loop
      term<-term*fact*rr/den
      vallast<-val
      val<-val+term
#      cat("term =",term,"  val now ",val,"\n")
      if (val == vallast) break
      fact<-(fact - 1)
      den<-den+1
      if (den > 1000) stop("Too many terms in mrate")
   }
   val*100
}
A<-function(I,Rval){
   A<-(1+I/100)^6-(1+R/200)
}

for (r2 in 0:24){
   R<-r2/2 # rate per year
   i.formula<-100*((1+R/200)^(1/6)-1)
   i.root<-uniroot(A,c(0,20),tol=.Machine$double.eps,Rval=R)$root
   i.mrate<-mrate(R)
   cat(R,"   ",i.mrate,"   Diffs:",
         i.formula-i.mrate,"  ",i.root-i.mrate,"\n")
}


###################################################
### code chunk number 13: chunk13
###################################################
cat("exponential example\n")
require("rootoned")
alpha<-1.0
efn<-function(x) { exp(-alpha*x) - 0.2 }
zfn<-function(x) { x*0 }
tint<-c(0,100)
curve(efn, from=tint[1], to=tint[2])
curve(zfn, add=TRUE, col='red')
rform<- -log(0.2)/alpha
resr<-root1d(efn,tint,tol=1e-10)
cat("alpha = ",alpha,"\n")
cat("root(formula)=",rform,"  root1d:",resr$root," tol=",resr$rtol,"  fval=",
         resr$froot,"  in   ",resr$fcount,"\n\n") 
alpha=0.02
curve(efn, add=TRUE, col='blue')
resr<-root1d(efn,tint,tol=1e-10, trace=TRUE)
rform<- -log(0.2)/alpha
cat("alpha = ",alpha,"\n")
cat("root(formula)=",rform,"  root1d:",resr$root," tol=",resr$rtol,"  fval=",
         resr$froot,"  in   ",resr$fcount,"\n\n")