rm(list=ls()) # e.g. Find the minimum of f(x) = x^2/10 -2*sin(x) in [0, 4]. f = function(x) { return(x^2/10 - 2*sin(x)) } # Golden section search minimizes a function f() whose opposite -f() # is unimodal over an interval [lower, upper]. That is, there exists # some minimizing value m such that x <= m implies f() is decreasing # and x >= m implies f() is increasing. This search finds m. # # Here is the code without any graphing. golden.section.search = function(f, lower, upper, tol=0.001, ...) { phi = (1 + sqrt(5))/2 # golden.ratio (allows cutting 1/2 of calls to f()) repeat { range = upper - lower lower.middle = upper - (phi - 1)*range upper.middle = lower + (phi - 1)*range if (range < tol) { return(lower) } if (f(lower.middle, ...) < f(upper.middle, ...)) { upper = upper.middle } else { lower = lower.middle } } } # Here is the code again, this time code added to make a graph and # show the progress of the search. golden.section.search = function(f, lower, upper, tol=0.001, ...) { phi = (1 + sqrt(5))/2 # golden.ratio (allows cutting 1/2 of calls to f()) curve(f, lower, upper) # for graph i = 0 # for graph repeat { range = upper - lower lower.middle = upper - (phi - 1)*range upper.middle = lower + (phi - 1)*range if (i < 4) { # for graph points(x=c(lower, lower.middle, upper.middle, upper), y=c(i,i,i,i), col=i+1) } if (range < tol) { abline(v=lower, col="orange") # for graph return(lower) } cat(sep="", "lower=", lower, ", lower.middle=", lower.middle, ", upper.middle=", upper.middle, ", upper=", upper, "\n") if (f(lower.middle, ...) < f(upper.middle, ...)) { upper = upper.middle } else { lower = lower.middle } i = i + 1 # for graph scan(what=character(), n=1, quiet=TRUE) # Require "Enter" to move. } } g = golden.section.search(f, 0, 4) print(g)