Class: Finrb::Numerical::Brent

Inherits:
Object
  • Object
show all
Defined in:
lib/finrb/numerical/brent.rb,
sig/finrb.rbs

Overview

Brent-Dekker root solver for a continuous function on a sign-changing interval. The interpolation steps are safeguarded by bisection.

Algorithm: R. P. Brent, Algorithms for Minimization Without Derivatives, Chapter 4 (1973). See also the GNU GSL root-finding documentation: https://www.gnu.org/software/gsl/doc/html/roots.html

Instance Method Summary collapse

Constructor Details

#initialize(tolerance:, relative_tolerance: tolerance, max_iterations: DEFAULT_MAX_ITERATIONS) ⇒ Brent

Returns a new instance of Brent.

Parameters:

  • tolerance: (number)
  • relative_tolerance: (number) (defaults to: tolerance)
  • max_iterations: (Integer) (defaults to: DEFAULT_MAX_ITERATIONS)

Raises:

  • (ArgumentError) 


18
19
20
21
22
23
24
25
 
# File 'lib/finrb/numerical/brent.rb', line 18

def initialize(tolerance:, relative_tolerance: tolerance, max_iterations: DEFAULT_MAX_ITERATIONS)
  @absolute_tolerance = decimal(tolerance)
  @relative_tolerance = decimal(relative_tolerance)
  @max_iterations = Integer(max_iterations)

  raise(ArgumentError, 'Tolerance must be positive.') unless @absolute_tolerance.positive? && @relative_tolerance.positive?
  raise(ArgumentError, 'Maximum iterations must be positive.') unless @max_iterations.positive?
end

Instance Method Details

#solve(function, lower:, upper:) ⇒ decimal

Parameters:

  • lower: (number)
  • upper: (number)

Returns:

  • (decimal)

Raises:

  • (ArgumentError) 


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
 
# File 'lib/finrb/numerical/brent.rb', line 27

def solve(function, lower:, upper:)
  left = decimal(lower)
  right = decimal(upper)
  raise(ArgumentError, 'Lower bound must be less than upper bound.') if left >= right

  left_value = evaluate(function, left)
  right_value = evaluate(function, right)

  return left if left_value.zero?
  return right if right_value.zero?
  raise(ConvergenceError, 'Root is not bracketed.') unless opposite_signs?(left_value, right_value)

  left, right, left_value, right_value = best_approximation_last(left, right, left_value, right_value)

  previous = left
  previous_value = left_value
  penultimate = previous
  bisected = true

  @max_iterations.times do
    tolerance = @absolute_tolerance + (@relative_tolerance * right.abs)
    return right if right_value.zero? || (right - left).abs <= tolerance

    candidate =
      if distinct_values?(left_value, right_value, previous_value)
        inverse_quadratic(left, right, previous, left_value, right_value, previous_value)
      else
        right - (right_value * (right - left) / (right_value - left_value))
      end

    bound = ((left * 3) + right) / 4
    outside_safe_interval = candidate <= [bound, right].min || candidate >= [bound, right].max
    insufficient_progress =
      (candidate - right).abs >= if bisected
                                   ((right - previous).abs / 2)
                                 else
                                   ((previous - penultimate).abs / 2)
                                 end
    bracket_too_small =
      if bisected
        (right - previous).abs < tolerance
      else
        (previous - penultimate).abs < tolerance
      end

    if outside_safe_interval || insufficient_progress || bracket_too_small
      candidate = (left + right) / 2
      bisected = true
    else
      bisected = false
    end

    candidate_value = evaluate(function, candidate)
    penultimate = previous
    previous = right
    previous_value = right_value

    if opposite_signs?(left_value, candidate_value)
      right = candidate
      right_value = candidate_value
    else
      left = candidate
      left_value = candidate_value
    end

    left, right, left_value, right_value = best_approximation_last(left, right, left_value, right_value)
  end

  raise(ConvergenceError, "Calculation did not converge after #{@max_iterations} iterations.")
end