Nothing
### 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")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.