22# author: Miguel Munoz Zuniga
33# ref: https://en.wikipedia.org/wiki/Brent%27s_method
44# tags: Inversion
5- # options: ytarget='0.0';ytol='3.e-8';xtol='1.e-8';max_iterations='100'
5+ # options: ytarget='0.0';ytol='0.1';xtol='0.01';max_iterations='100'
6+ # options.help: ytarget='Output target value';ytol='Convergence precision on output value';xtol='Convergence precision on input value';max_iterations='Maximum iterations number'
67# input: x=list(min=0,max=1)
78# output: y=0.01
89
@@ -25,6 +26,13 @@ getInitialDesign <- function(brent, input, output) {
2526 if (length(input )!= 1 ) stop(" Cannot find root of more than 1D function" )
2627 brent $ i <- 0
2728 brent $ input <- input
29+
30+ # Rescale xtol in [0,1]
31+ Xname = names(brent $ input )[1 ]
32+ xminptol = matrix (brent $ input [[ Xname ]]$ min + brent $ xtol ,ncol = 1 )
33+ names(xminptol ) <- Xname
34+ brent $ xtol01 = to01(xminptol ,brent $ input ) # Rescale xtol
35+
2836 brent $ exit <- - 1 # Reason end of algo
2937 x = matrix (c(0 , 1 , 1 ),ncol = 1 )
3038 names(x ) <- names(input )
@@ -77,7 +85,8 @@ getNextDesign <- function(brent, X, Y) {
7785 fc = fa
7886 }
7987
80- tol1 = 2 . * brent $ ytol * abs(b ) + 0.5 * brent $ xtol # Convergence check tolerance.
88+ # tol1 = 2. * brent$ytol * abs(b) + 0.5 * brent$xtol01 # Convergence check tolerance.
89+ tol1 = 0.5 * brent $ xtol01 # Convergence check tolerance.
8190 xm = .5 * (c - b )
8291 if ((abs(xm ) < = tol1 ) | (fb == 0 )) {
8392 # stop if fb = 0 return root b or tolerance reached
@@ -186,15 +195,15 @@ to01 = function(X, inp) {
186195# #############################################################################################
187196# @test
188197# f <- function(X) matrix(Vectorize(function(x) {((x+5)/15)^3})(X),ncol=1)
189- #
198+ #
190199# options = list(ytarget=0.3,ytol=3.e-8,xtol=1.e-8,max_iterations=100)
191200# b = Brent(options)
192- #
201+ #
193202# X0 = getInitialDesign(b, input=list(x=list(min=-5,max=10)), NULL)
194203# Y0 = f(X0)
195204# Xi = X0
196205# Yi = Y0
197- #
206+ #
198207# finished = FALSE
199208# while (!finished) {
200209# Xj = getNextDesign(b,Xi,Yi)
@@ -206,5 +215,5 @@ to01 = function(X, inp) {
206215# Yi = rbind(Yi,Yj)
207216# }
208217# }
209- #
210- # print(displayResults(b,Xi,Yi))
218+ #
219+ # print(displayResults(b,Xi,Yi))
0 commit comments