Simulataneously calculate a function and it's derivitive using Dual Numbers

dual_number_derivative.rb

require 'matrix'

##  Note: In ruby the notation x**n means x raised to the nth power
##
##  In linear algebra, Dual numbers are in the form:
##  a + bε
##  
##  Where ε is what is called nilpotent, you can think of
##  ε as an infinitesimal (sort of kinda), like how they used to do Calculus in the
##  olden days.
##
##  Nilpotent means there is some number n where x**n = 0
##  
##  Dual numbers are similar to Complex numbers in the form a + bi where
##  i * i = -1.  A complex number is a 2 dimensional number, and 
##  so are dual numbers, but the ε dimension is infinitesimally small.
##
##  Where i * i = -1,  instead, in dual numbers ε * ε = 0
##  
##  Complex Numbers as Matrices:
##  a + bi    = | a, b|
##              |-b, a|
## 
##  So 0 + 1i = | 0, 1|
##              |-1, 0|
##
##  So i * i  = | 0, 1| * | 0, 1| = | 0 * 0 + 1 * -1,  0 * 1 + 1 * 0|
##              |-1, 0|   |-1, 0|   |-1 * 0 + 0 * -1, -1 * 1 + 0 * 0|
##
##            = |-1,  0|
##              | 0, -1|
##
##  Therefore: i * i = -1
##
##  Dual Numbers as Matrices:
##
##  If we have a + bε, and instead of setting element [0,1] to -b, we
##  set it to 0, we get a dual number
##  
##  a + bε    = |a, b|
##              |0, a|
##
##  So 0 + 1ε = |0, 1|
##              |0, 0|
##
##  So ε * ε  = |0, 1| * |0, 1| = |0 * 0 + 1 * 0, 0 * 1 + 1 * 0|
##              |0, 0|   |0, 0|   |0 * 0 + 0 * 0, 0 * 1 + 0 * 0|
##
##            = |0, 0|
##              |0, 0|
##
##  So ε * ε = 0 

####
##  Here is a class that implements Dual numbers, with addition,
##  multiplication, and exponentiation
class Dual
    attr_accessor :matrix

    def initialize(real, dual)
        @matrix = Matrix[[real, dual], [0, real]]
    end

    def +(other)
        result = @matrix + other.matrix
        Dual.new(result[0, 0], result[0, 1])
    end

    def *(other)
        result = @matrix * other.matrix
        Dual.new(result[0, 0], result[0, 1])
    end

    def **(exponent)
        result = @matrix**exponent
        Dual.new(result[0, 0], result[0, 1])
    end

    def to_s
        "#{real} + #{dual}ε"
    end

    def inspect
        to_s
    end

    def real
        @matrix[0, 0]
    end

    def dual
        @matrix[0, 1]
    end

end


####
##  Now let's have f(x) = x**2
def f(x)
    x**2
end


####
##  Let's print a few values of this function
print "f(x) = "
p (0..15).map{|x| f(x) }


## This prints:
## f(x) = [0, 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169, 196, 225]


####
##  Now lets run function f on a dual number, f(x + ε)
print "f(x + ε) = "
p (0..15).map{|x| f(Dual.new(x, 1)) }


##  This prints
##  f(x + ε) = [0 + 0ε, 1 + 2ε, 4 + 4ε, 9 + 6ε, 16 + 8ε, 25 + 10ε, 36 + 12ε, 
##              49 + 14ε, 64 + 16ε, 81 + 18ε, 100 + 20ε, 121 + 22ε, 144 + 24ε, 
##              169 + 26ε, 196 + 28ε, 225 + 30ε]


##  The real part of the dual number in f(x + ε) is exactly the same as f(x)
##  Let's just print out the dual component instead

print "Dual{f(x + ε)} = "
p (0..15).map{|x| f(Dual.new(x, 1)).dual }


##  This prints
##  Dual{f(x + ε)} = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30]
##
##  It looks like Dual{f(x + ε)} = f'(x) = 2 * x
##
##  In Calculus, the derivitive of the function f(x) = x**2, is f'(x) = 2 * x
##
##  So, what dual numbers are good for... is that they simultaneously calculate
##  a function f(x) and it's derivative f'(x), and this works for any function.
##
##  Tada!

dual_numbers_newton_raphson.rb

require 'matrix'

##  An update of my dual number library, used to find roots using
##  Newton-Raphson method, using Dual number's automatic differentiation

##  Note: In ruby the notation x**n means x raised to the nth power
##
##  In linear algebra, Dual numbers are in the form:
##  a + bε
##
##  Where ε is what is called nilpotent, you can think of
##  ε as an infinitesimal (sort of kinda), like how they used to do Calculus in the
##  olden days.
##
##  Nilpotent means there is some number n where x**n = 0
##
##  Dual numbers are similar to Complex numbers in the form a + bi where
##  i * i = -1.  A complex number is a 2 dimensional number, and
##  so are dual numbers, but the ε dimension is infinitesimally small.
##
##  Where i * i = -1,  instead, in dual numbers ε * ε = 0
##
##  Complex Numbers as Matrices:
##  a + bi    = | a, b|
##              |-b, a|
##
##  So 0 + 1i = | 0, 1|
##              |-1, 0|
##
##  So i * i  = | 0, 1| * | 0, 1| = | 0 * 0 + 1 * -1,  0 * 1 + 1 * 0|
##              |-1, 0|   |-1, 0|   |-1 * 0 + 0 * -1, -1 * 1 + 0 * 0|
##
##            = |-1,  0|
##              | 0, -1|
##
##  Therefore: i * i = -1
##
##  Dual Numbers as Matrices:
##
##  If we have a + bε, and instead of setting element [0,1] to -b, we
##  set it to 0, we get a dual number
##
##  a + bε    = |a, b|
##              |0, a|
##
##  So 0 + 1ε = |0, 1|
##              |0, 0|
##
##  So ε * ε  = |0, 1| * |0, 1| = |0 * 0 + 1 * 0, 0 * 1 + 1 * 0|
##              |0, 0|   |0, 0|   |0 * 0 + 0 * 0, 0 * 1 + 0 * 0|
##
##            = |0, 0|
##              |0, 0|
##
##  So ε * ε = 0

class Numeric
  def to_dual
    self.kind_of?(Dual) ? self : Dual.new(self, 0)
  end
end


####
##  Here is a class that implements Dual numbers, with addition,
##  multiplication, and exponentiation
class Dual < Numeric
    attr_accessor :matrix

    def initialize(real, dual)
        @matrix = Matrix[[real, dual], [0, real]]
    end

    def +(other)
        result = @matrix + other.to_dual.matrix
        Dual.new(result[0, 0], result[0, 1])
    end

    def -(other)
        result = @matrix - other.to_dual.matrix
        Dual.new(result[0, 0], result[0, 1])
    end

    def *(other)
        result = @matrix * other.to_dual.matrix
        Dual.new(result[0, 0], result[0, 1])
    end

    def /(other)
      fail('Think about how this works')
    end

    def **(exponent)
        result = @matrix**exponent
        Dual.new(result[0, 0], result[0, 1])
    end

    def ratio
      Rational(self.real, self.dual)
    end

    def to_s
        "#{real} + #{dual}ε"
    end

    def inspect
        to_s
    end

    def real
        @matrix[0, 0]
    end

    def dual
        @matrix[0, 1]
    end

end

##  Pick a starting point x = 0.11 + ε
x = Dual.new(0.11, 1)
##  Define the non-linear equation f(x) you want to find a root for
f = lambda{|x| 4.to_dual * x**5 + x**2 -3}

##  Iterate on this a few times, The dual numbers calculate f(x) and f'(f)
##  simulataneously in the real and dual parts.
##  So Newton-Raphson gives you a better approximation for the root
##  each iteration by x - f(x) / f'(x)
##  Dual#ratio here is defined as the real part over the dual part,
##  same as f(x) / f'(x).
15.times do
  f0 = f[x]
  x = x - f0.ratio
end

puts "Approximate root #{x.real}"