https://github.com/cran/nleqslv
Raw File
Tip revision: 435494ed15bede0e00c105012679e77b4a0215d2 authored by Berend Hasselman on 08 August 2016, 17:19:53 UTC
version 3.0.3
Tip revision: 435494e
chquad.Rout.save

R version 2.15.2 Patched (2013-01-21 r61728) -- "Trick or Treat"
Copyright (C) 2013 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 
> # Chebyquad functions (no solution for n=8)
> 
> library("nleqslv")
> 
> chebyquad <- function(x) {
+ 	n <- length(x)
+ 	y <- numeric(n)
+ 
+ 	for(j in 1:n) {
+ 		t1 <- 1.0
+ 		t2 <- 2.0*x[j] - 1.0
+ 		tmp <- 2.0*t2
+ 
+     	for(i in 1:n) {
+ 			y[i] <- y[i] + t2
+ 			t3 <- tmp * t2 - t1
+ 			t1 <- t2
+ 			t2 <- t3
+ 		}
+ 	}
+ 
+ 	y = y / n
+ 
+ 	for(i in 1:n) {
+     	if ( i%%2 == 0 ) {
+ 			y[i] = y[i] + 1.0 / (i * i - 1)
+ 		}
+ 	}
+ 
+ 	y
+ }
> 
> chebyinit <- function(n) {
+     x <- (1:n) / (n + 1)
+ }
> 
> for (n in c(1:7,9)) {  # exclude n=8 since there is no solution
+ 	xstart <- chebyinit(n)
+ 	fstart <- chebyquad(xstart)
+ 
+ 	zz <- nleqslv(xstart, chebyquad, global="dbldog",
+ 	              control=list(ftol=1e-8,xtol=1e-8,trace=0,btol=.01,delta=-2))
+ 	print(all(abs(zz$fvec)<=1e-8))
+ }
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
> 
> for (n in c(1:7,9)) {  # exclude n=8 since there is no solution
+ 	xstart <- chebyinit(n)
+ 	fstart <- chebyquad(xstart)
+ 
+ 	zz <- nleqslv(xstart, chebyquad, global="dbldog", method="Newton",
+ 	              control=list(ftol=1e-8,xtol=1e-8,trace=0,btol=.01,delta=-2))
+ 	print(all(abs(zz$fvec)<=1e-8))
+ }
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
[1] TRUE
> 
> proc.time()
   user  system elapsed 
  0.374   0.046   0.401 
back to top