automatic_differentiation(n) 1.0 "flibs"
automatic_differentiation  Implement the "automatic differentation" technique
TABLE OF CONTENTS
SYNOPSIS
DESCRIPTION
EXAMPLE
DATA TYPES AND ROUTINES
POSSIBLE EXTENSIONS
COPYRIGHT
The automatic_differentiation module defines a derived
type that enables you (via overloading common mathematical
functions and operators) to:

Define a mathematical function in the usual way

Evaluate that function and its first derivative at the same
time
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.)
In numerical methods like NewtonRaphson 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

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 NewtonRaphson
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 NewtonRaphson
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.
You can extend the module in at least two ways:

Determine the second derivative, the third and so on. This is
straightforward enough, but the formulae will get fairly complex after
the second derivative.

Determine the derivative of a function of two or more variables. This
will require a bit more: such functions have a vectorvalued
derivative and each independent variable will have to have a
vectorvalued derivative as well. For instance:
A function f(x,y) evaluated at (1,2), would take as the actual arguments
x = (1,1,0) (so no variation in the second direction) and y = (2,0,1)
(no variation in the first direction).
Copyright © 2006 Arjen Markus <arjenmarkus@sourceforge.net>