Newton-Raphson Method

Description of the problem

Today, I am going to look into one of the most fundamental numerical algorithms: Newton-Raphson method. It is so fundamental in fact that we teach it in calculus classes.

The main question is to find the set of solutions for an equation of the type \[ f(x) = 0 \] where $f(x)$ is a nice (nice in this context means differentiable) function. The Newton-Raphson algorithm tells us that the recursively defined sequence \[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \] with an initial guess $x_0$ for the solution will converge to one of the solutions provided that your initial guess was close enough.

I am not going to explain the mathematical reasoning behind the algorithm. This is an interesting question on its own, but I am here to write some lisp code.

An implementation

First, I will need a function which numerically differentiates a numerical function and returns the derivative function.

(defun differentiate (f &optional (epsilon 0.00001d0))
    (lambda (x) (/ (- (funcall f (+ x epsilon))
                      (funcall f (- x epsilon)))
                   (* 2.0d0 epsilon))))

Granted, I can write a simpler function which returns the number which is the derivative of a function at a point (which is all I need) but I like this better.

Now, let me test if this works as intended. First, I will create a random list of real numbers and then apply the derivative of $\sin(x)$ (which is $\cos(x)$) to this random list.

(defvar random-list 
        (loop for i from 0 to 9 collect (random 1.7)))

(format nil "~{~2,5F ~}~%" (mapcar (differentiate 'sin)  random-list))
.48007 .89578 .18994 .41848 .90823 -.02335 .16457 .73065 .90044 .49951 

(format nil "~{~2,5F ~}~%" (mapcar 'cos  random-list))
.48007 .89578 .18994 .41848 .90823 -.02335 .16457 .73065 .90044 .49951

Let me test it on a different function and on a different random list:

(setf random-list 
        (loop for i from 0 to 9 collect (- (random 2.1) (random 4.2))))

(defvar fn (lambda (x) (- (expt x 3.0d0) (* x x) (* 2.0d0 x) 10.0d0)))
(defvar dfn (lambda (x) (- (* 3.0d0 x x) (* 2.0d0 x) 2.0d0)))

(format nil "~{~2,5F ~}~%" (mapcar (differentiate fn)  random-list))
-2.02553 -.66356 52.53546 17.66807 13.66718 5.04636 14.92025 -1.01093 -.41888 5.03504 

(format nil "~{~2,5F ~}~%" (mapcar dfn  random-list))
-2.02553 -.66356 52.53546 17.66807 13.66718 5.04636 14.92025 -1.01093 -.41888 5.03504

Not bad… It is time for us to implement the recursive sequence:

(defun nr-iterate (f x)
   (- x (/ (funcall f x) (funcall (differentiate f) x))))

And now the main function which will return a solution within given tolerance.

(defun newton-raphson (f initial-guess &optional (tolerance 0.00001d0))
   (do* ((x initial-guess y)
         (y (nr-iterate f x) (nr-iterate f y)))
      ((< (abs (- x y)) tolerance) y)))

Let me test it: I will try to find a solution to the equation $e^{-x^2}-\ln(x^2+x+1)-x = 0$

(setf fn (lambda (x) (- (exp (- (* x x))) (log (+ (* x x) x 1)) x)))

(defvar solution (newton-raphson fn 1.0d0))
(format nil "~3,9F~%" solution)
.402628809

(format nil "~3,9F~%" (funcall fn solution))
-.000000000