Class: Minimization::BrentRootFinder

Inherits:
Object
  • Object
show all
Defined in:
lib/multidim/brent_root_finder.rb

Overview

Brent's root finding method.

Usage

f = lambda {|x| x**2} root_finder = Minimization::BrentRootFinder.new(100000) min_x = root_finder.find_root(-1000,1000, f)

Constant Summary collapse

MAX_ITERATIONS_DEFAULT = 
10e6
EPSILON = 
10e-10

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(max_iterations = nil) ⇒ BrentRootFinder 



39
40
41
42
43
44
 
# File 'lib/multidim/brent_root_finder.rb', line 39

def initialize(max_iterations=nil)
  @iterations = 0
  if (@max_iterations.nil?)
    @max_iterations = MAX_ITERATIONS_DEFAULT
  end
end

Instance Attribute Details

#max_iterations ⇒ Object

Returns the value of attribute max_iterations.



35
36
37
 
# File 'lib/multidim/brent_root_finder.rb', line 35

def max_iterations
  @max_iterations
end

Instance Method Details

#f(x) ⇒ Object 



46
47
48
 
# File 'lib/multidim/brent_root_finder.rb', line 46

def f(x)
  return @f.call(x)
end

#find_root(lower_bound, upper_bound, f) ⇒ Object

Find root in interval (lower_bound, upper_bound) of function f

Parameters:

  • lower_bound: Lower bound of the minimization search
  • upper_bound: Upper bound of the minimization search
  • f: Function to find roots


56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
 
# File 'lib/multidim/brent_root_finder.rb', line 56

def find_root(lower_bound, upper_bound, f)
  @f = f
  lower  = lower_bound
  f_upper = f(lower_bound)
  upper  = upper_bound
  f_lower = f(upper_bound)
  c  = lower
  fc = f_upper
  d  = upper - lower
  e  = d

  absolute_accuracy = EPSILON
  relative_accuracy = EPSILON

  loop do
    @iterations += 1
    if (fc.abs < f_lower.abs)
      lower  = upper
      upper  = c
      c  = lower
      f_upper = f_lower
      f_lower = fc
      fc = f_upper
    end

    tolerance = 2 * relative_accuracy * upper.abs + absolute_accuracy
    m   = 0.5 * (c - upper)

    if (m.abs <= tolerance or f_lower.abs < EPSILON or @iterations > @max_iterations)
      return upper
    end
    if (e.abs < tolerance or f_upper.abs <= f_lower.abs)
      # use bisection
      d = m
      e = d
    else 
      # use inverse cubic interpolation
      s = f_lower / f_upper
      if (lower == c)
        p = 2 * m * s
        q = 1 - s
      else
        q = f_upper / fc
        r = f_lower / fc
        p = s * (2 * m * q * (q - r) - (upper - lower) * (r - 1))
        q = (q - 1) * (r - 1) * (s - 1)
      end
      if (p > 0)
        q = -q
      else 
        p = -p
      end
      s = e
      e = d
      if (p >= 1.5 * m * q - (tolerance * q).abs or p >= (0.5 * s * q).abs)
        # interpolation failed, fall back to bisection
        d = m
        e = d
      else
        d = p / q
      end
    end
    # Update the best estimate of the root and bounds on each iteration
    lower  = upper
    f_upper = f_lower

    if (d.abs > tolerance)
      upper += d
    elsif (m > 0)
      upper += tolerance
    else
      upper -= tolerance
    end
    f_lower = f(upper)
    if ((f_lower > 0 and fc > 0) or (f_lower <= 0 and fc <= 0))
      c  = lower
      fc = f_upper
      d  = upper - lower
      e  = d
    end
  end
end