1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
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
138
139
|
# = brent_root_finder.rb -
# Minimization- Minimization algorithms on pure Ruby
# Copyright (C) 2010 Claudio Bustos
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
#
# This algorith was adopted and ported into Ruby from GNU GSL library's
# brent.c [https://github.com/ampl/gsl/blob/master/roots/brent.c]. Therefore this
# file is under under the terms of the GNU General Public License.
# Research Paper - http://maths-people.anu.edu.au/~brent/pd/rpb005.pdf
#
# Brent method for finding roots
module Minimization
# 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)
#
class BrentRootFinder
attr_accessor :max_iterations
MAX_ITERATIONS_DEFAULT = 10e6
EPSILON = 10e-10
def initialize(max_iterations=nil)
@iterations = 0
if (@max_iterations.nil?)
@max_iterations = MAX_ITERATIONS_DEFAULT
end
end
def f(x)
return @f.call(x)
end
# Find root in interval (lower_bound, upper_bound) of function f
# == Parameters:
# * <tt>lower_bound</tt>: Lower bound of the minimization search
# * <tt>upper_bound</tt>: Upper bound of the minimization search
# * <tt>f</tt>: Function to find roots
#
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
end
end
|
