http://math.stackexchange.com/q/2071775/86776

planet.rb

# planet.rb
#
# Planet Ship Intersection
#
# http://math.stackexchange.com/q/2071775/86776

require "matrix.rb"

# problem instance
$x0 = 1.0
$y0 = 2.0
$u0 = Vector[$x0, $y0]

$a = 2.0
$b = 1.0

$v = 1.0
$w = 1.0

# Newton Raphson iteration:
#
# Scalar, one argument:
#  0 = f(x_n) + f'(x_n) (x_{n+1} - x_n)
#  x_{n+1} = x_n - f(x_n) / f'(x_n)
#
# General case:
#  0 = f(x_n) + J(x_n) (x_{n+1} - x_n)
#  J_{ij} = \partial f_i / \partial x_j

def f(u)
  phi = u[0]
  t = u[1]
  Vector[$x0 + $v * Math.cos(phi) * t - $a * Math.cos($w * t),
         $y0 + $v * Math.sin(phi) * t - $b * Math.sin($w * t)]
end

def jacobian(u)
  phi = u[0]
  t = u[1]
  Matrix.rows([[-$v * Math.sin(phi) * t,
                $v * Math.cos(phi) + $a * Math.sin($w * t) * $w],
               [$v * Math.cos(phi) * t,
                $v * Math.sin(phi) - $b * Math.cos($w * t) * $w]])
end


def newton(u0, eps, maxstep)
  puts "x0 = #{$x0}"
  puts "y0 = #{$y0}"
  puts "a  = #{$a}"
  puts "b  = #{$b}"
  puts "v  = #{$v}"
  puts "w  = #{$w}"
  k = 0
  u = u0
  fu = f(u)
  d = nil
  puts "#{k}: u=#{u} f(u)=#{fu}"
  loop do
    jac = jacobian(u)
    if not jac.regular?
      puts "Jacobian #{jac} is not regular at u=#{u}!"
      return nil
    end
    u2 = u - jac.inverse * fu
    u = u2
    fu = f(u)
    d = fu.magnitude
    k += 1
    puts "#{k}: u=#{u} f(u)=#{fu} d=#{d}"
    break if d < eps or k >= maxstep
  end
  if d < eps
    puts "d=#{d} < eps=#{eps}, success"
    return u
  end
  puts "maxstep = #{maxstep} reached, fail"
  nil
end