The Newton-Raphson Method in Ruby

README.md

Newton-Raphson Method in Ruby

As part of an exercise in figuring out how to derive the continuous equation that is sampled by a discrete function, I went ahead and opened my Numerical Methods textbook and started going through the Newton-Raphson Method. This, plus my first ever tweet storm are the output.

newton.rb

#!/usr/bin/env ruby

# Represents a set in the Newton-Raphson Method.
Step = Struct.new(:iteration, :x_n, :fx, :fp) do
  def to_s
    [
      iteration.to_s.rjust(3),
      x_n.round(5).to_s.rjust(12),
      fx.round(5).to_s.rjust(12),
      fp.round(5).to_s.rjust(12)
    ].join
  end
end

# Represents a final result in the Newton-Raphson Method.
Result = Struct.new(:message) do
  def to_s
    length = message.length
    total_length = 39
    padding = ' ' * ((total_length - length) / 2)

    [
      '-' * total_length,
      "#{padding}#{message.upcase}#{padding}"
    ].join("\n")
  end
end

# Run the Newton-Raphson Method algorithm for a given set of parameters.
#
# @param x_0 [Numeric] - the initial estimation of the root
# @param func [#call] - the function for which we are estimating a root
# @param derivative [#call] - the derivative of `func`
# @param max_iterations [Integer] - the maximum number of iterations to try to find the root
# @param epsilon [Numeric] - the threshold under which we presume convergence
# @param delta [Numeric] - the threshold under which we presume we cannot get closer
#
# @return [Array<Step|Result>] the lines to output
def newton(x_0, func, derivative, max_iterations, epsilon = 0.01, delta = 0.01)
  output_lines = []
  max_iterations = Integer(max_iterations)
  x_n = x_0

  (0..max_iterations).each do |iteration|
    f_x = func.(x_n)
    f_prime = derivative.(x_n)

    if f_prime.abs < delta
      output_lines << Result.new('small derivative')
      break
    end

    output_lines << Step.new(iteration, x_n, f_x, f_prime)

    difference = f_x / f_prime
    x_n -= difference

    if difference.abs < epsilon
      output_lines << Result.new('convergence')
      break
    end
  end

  output_lines
end

# -----------------------------------------------------------------------
# Here, we set up the functions we want to test. We don't have a symbolic
# interpreter so you have to derive the derivative yourself. Huzzah,
# calculus!

f = ->(x) { x**3 - x + 1.0 }
f_prime = ->(x) { 3.0 * x**2 - 1 }

# Run the algorithm and assemble the output
output = newton(1.0, f, f_prime, 20)

#------------------------------------------------------------------------
# Past this is just output code. Feel free to ignore it

output << Result.new('did not converge') unless output.last.is_a?(Result)

puts ['n'.rjust(3), 'x_n'.rjust(12), 'f(x_n)'.rjust(12), "f'(x_n)".rjust(12)].join
puts '-' * 39
output.each { |step| puts step }