Newton-Raphson method

Solving equations with the Newton-Raphson method

by Stephen R. Schmitt

Introduction

The Newton-Raphson method allows one to solve equations of the form f(x) = 0 by finding the values of x for which the equality is valid.

In 1669, Isaac Newton found an algorithm to solve for the roots (values for which the function equals zero) of a polynomial equation. In this method, one guesses a starting value and then repeatedly makes small changes that yield improved approximations to the solution. The process terminates once the desired precision is reached. Newton's original method did not use the derivative of f(x).

Later, in 1690, Joseph Raphson found an improvement to Newton's method which did use the derivative of f(x), f'(x). Each iterative step of the Newton-Raphson method is

x_n+1 := x_n - f(x_n)/f'(x_n)

If an estimated root of f(x) = 0 is x_n, then the line tangent to f(x_n) at x_n crosses the x-axis at a point, x_n+1, which is an improved estimate of the root. The slope of this tangent is given by the derivative. Repeated application of the iterative step improves the estimate.

Approximate derivative

Sometimes it is inconvenient to explicitly determine the derivative of a function for use in a computer program. An approximate derivative of f(x) at x can be computed from

        f(x+h) - f(x-h)
f'(x) ≈ ---------------
              2·h

where h is small.

The steps to apply Newton-Raphson method to find the root of an equation f(x) = 0 using an approximate derivative are

Pick an arbitrary starting estimate, x_n.
Compute an approximate value of f'(x_n) and the value of f(x_n).
Perform the iterative step calculation to compute x_n+1.
Repeat steps 2 and 3 until f(x_n)/f'(x_n) in close enough to zero.

Note: To avoids infinite loops, it is good practice to ensure that the number of iterations has not exceeded some arbitrary maximum.

Zeno source code

Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.

Download Page

const TOLERANCE : real := 0.00000001
const MAX_ITER  : int  := 50

program

    var root : real
    var num_roots, i : int

    put "solving: f(x) = exp(x) - 3*x^2"
    
    repeat
    
        put "initial guess"...
        get root
    
        if newton_approx( root ) then      
            put "root = ", root
        else
            put "could not solve"
        end if

    until not another      
    
end program

% f(x) = exp(x) - 3*x^2
function Fx( x : real ) : real

    return exp(x) - 3*x^2

end function

%
% Solve for a root of a function of a single variable by starting
% from an initial guess and using the function value and an
% approximate derivative to search for the solution
%
function newton_approx( var root : real ) : boolean
  
    var iter : int := 0
    var h, diff : real
    
    repeat
    
        h := 0.01 * root

        % root could be small or zero
        if abs( root ) < 1 then 
            h := 0.01
        end if
        
        % calculate guess refinement
        diff := 2*h*Fx(root) / (Fx(root + h) - Fx(root - h))
        
        % update guess
        root := root - diff
        iter := iter + 1

    until ( iter > MAX_ITER ) or ( abs( diff ) < TOLERANCE )
  
    if abs( diff ) <= TOLERANCE then
        return true
    end if

    return false

end function

% returns the absolute value of the argument
function abs( x : real ) : real

    if x < 0.0 then
        x := -x
    end if

    return x

end function

% ask user to continue or not
function another : boolean

    var ans : string
    put "another [Y|N]"...
    get ans
    return (ans[1] = 'Y') or (ans[1] = 'y')

end function

Sample output

solving: f(x) = exp(x) - 3*x^2
initial guess? 0
root = -0.458962
another [Y|N]? y
initial guess? 1
root = 0.910008
another [Y|N]? y
initial guess? 2
root = 0.910008
another [Y|N]? y
initial guess? 3
root = 3.73308
another [Y|N]? n

AbCd Classics - free on-line books