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 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
|
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.
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 vector-valued
derivative and each independent variable will have to have a
vector-valued 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>