automatic_differentiation(n) 1.0 "flibs"

NAME

automatic_differentiation - Implement the "automatic differentation" technique

TABLE OF CONTENTS

    TABLE OF CONTENTS
    SYNOPSIS
    DESCRIPTION
    EXAMPLE
    DATA TYPES AND ROUTINES
    POSSIBLE EXTENSIONS
    COPYRIGHT

SYNOPSIS

use automatic_differentiation
type(AUTODERIV)
x = derivvar( value )
call find_root( f, xinit, tolerance, root, found )
call find_root_iter( fvalue, root, tolerance, found )

DESCRIPTION

The automatic_differentiation module defines a derived type that enables you (via overloading common mathematical functions and operators) to:

without having to program that first derivative.

The module uses real numbers of the kind "wp" as defined in the auxiliary module select_precision.

(I was pointed to this technique by Simon Geard, a couple of years ago, in the context of one of the many language shootouts on the Internet.)

EXAMPLE

In numerical methods like Newton-Raphson for finding a root of the equation "f(x) = 0", you need the first derivative:

 
    x(k+1)  =  x(k)  - f(x(k)) / f'(x(k))

Rather than implementing the function and its first derivatives as separate functions, using this module you can dispense with manually determining the first derivative:

 
program root
    use automatic_differentation

    type(AUTODERIV) :: x
    type(AUTODERIV) :: res
    integer         :: i

    !
    ! First estimate
    !
    x = derivvar( 1.0_wp )

    do i = 1,10
        res = f(x)
        x = x - res.v / res.dv
    enddo

    write(*,*) 'Root: ', x.v
contains

type(AUTODERIV) function f(x)
    type(AUTODERIV) :: x

    f = x + cos(x)

end function f
end program

DATA TYPES AND ROUTINES

The module defines a single data type, AUTODERIV, and one specific function, derivvar().

use automatic_differentiation
The name of the module

type(AUTODERIV)
The type has two fields:

real(wp) v
The value of the variable/function (or any intermediate result)

real(wp) dv
The first derivative


x = derivvar( value )
Use this function to properly initialise the program variable x as a "mathematical" variable. (It sets the derivative of x to 1, because otherwise it would be regarded as a constant).

real(wp) value
The value of the "mathematical" variable.
call find_root( f, xinit, tolerance, root, found )
This subroutine is a simple implementation of the Newton-Raphson method to find roots. The function f takes one argument, x.

type(AUTODERIV) function f(x)
The function must have the interface:

 
    interface
        function f(x)
            use automatic_differentiation_type
            type(AUTODERIV), intent(in) :: x
            type(AUTODERIV)             :: f
        end function
    end interface

Its return value is the function value at point x and its first derivative.

type(AUTODERIV), intent(in) xinit
Initial estimate of the root - the end result may depend on this choice.

real(wp) tolerance
Maximum allowable error in the absolute value of the root. If the difference between the old and the new estimate is smaller than this value, the iteration stops.

type(AUTODERIV), intent(out) root
Final estimate of the root.

logical found
Indicator whether a root was found or not. There is a maximum number of iterations and the root must not grow too large.
call find_root_iter( fvalue, root, tolerance, found )
This subroutine performs a single iteration in the Newton-Raphson method to find roots. It can be used to implement a more general version of find_root, for instance if the interface to the function contains one or more asjustable parameters.

type(AUTODERIV) fvalue
The value of the function (plus its first derivative)

type(AUTODERIV), intent(in) root
Current estimate of the root

real(wp) tolerance
Maximum allowable error in the absolute value of the root. If the difference between the old and the new estimate is smaller than this value, the indicator found is set to true.

logical found
Indicator whether a root was found or not.

POSSIBLE EXTENSIONS

You can extend the module in at least two ways:

COPYRIGHT

Copyright © 2006 Arjen Markus <arjenmarkus@sourceforge.net>