In the fall of 1972 President Nixon announced that the rate of increase of inflation was decreasing. This was the first time a sitting president used the third derivative to advance his case for reelection.

Hugo Rossi. Mathematics Is an Edifice, Not a Toolbox, Notices of the AMS, v. 43, no. 10, October 1996.

Welcome to the third part is the third part of the essential introduction to the FORTRAN programming series. In this part, we will introduce the use of functions and subroutines to structure the programs. Then we use them to make a program to calculate numerical derivatives of mathematical functions.

FUNCTIONS

We are already familiar with intrinsic functions. As the name say, these are function that come with the FORTRAN compiler as a an intrisic mathematical library and directly accessible to the program. In addition to this rich but limited set of functions, the programmer can define statments functions and user defined functions. The first one are single line defined functions that are useful to incpasulate a complex formula.

The are defined as

name(argument list ) = expression

Example

The followin function return the distance travelled by a body with an initial speed V0 and constant acceleration (A) after a time (t)

PROGRAM ACCEL_MOTION
    REAL A,V0,t,DIST
    DIST(V0,A,t)=0.5*A*t*t+V0*t

    A=30.
    V0=2.
    t=100.
    PRINT 20,DIST(A,V0,t)

    20 FORMAT("DIST=",F12.5)
END PROGRAM ACCEL_MOTION

To use in the program you just need to provide the value for the constants as

print *,DIST(A,V0,t)

If the unit are meters and seconds, the program will print

DIST=13000

The user-define functions are useful when the expression defining the function is complex or is calculated using a block of staments (e.g. piecewise defined mathematical function. Such functions are definied with the following (essential ) syntax.

SYNTAX

type FUNCTION  name(argument list ) 

      [ STATEMENTs BLOCK] 

      name=[ STATMENTS]
RETURN  

END FUNCTION name

Example

PROGRAM PIECEWISE_FUNCT
    REAL :: x

    x=0.8
    PRINT 20,FN(x)

    20 FORMAT("F(X)=",F12.5)
END  PROGRAM PIECEWISE_FUNCT

REAL FUNCTION FN(x)

    REAL :: x

    IF (x>0 .AND. x< 1.57) THEN
        FN=SIN(x)
    ELSE IF (x<0) THEN
        FN=0.0
    ELSE
        FN =1.0
    END IF

    RETURN

END FUNCTION FN

SUBROUTINES

The subroutine is a portion of the program that cannot be executed independently but only as part of another program. As for the user-defined function, the subroutine contains a specific algorithm or a stand-alone procedure. Differently from the user-defined functions, there is no type to declare at the beginning of the definition, and there is no value associated with the subroutine name. On the contrary, the argument list contains both input and output variables. Finally, the subroutine is called using the CALL command.
In modern FORTRAN compilers, the possibility of checking the intent of each variable passed to a function or a subroutine using the command INTENT has been added. However, by default, the compiler does not perform such checking that it can help to detect erroneous coding. Instead, the INDENT is done at the variable declaration by specifying with INTENT(IN) that the value of the dummy argument may be used, but not modified, within the procedure; INTENT(OUT) that the dummy argument may be set and then changed within the subroutine, and the values returned to the caller. Finally, with INTENT(INOUT), the dummy argument values may be used and modified within the procedure and returned to the caller.

EXAMPLE: CALCULATION OF NUMERICAL DERIVATIVES

The numerical calculation of derivative is a standard procedure used, for example, in optimization problems (finding maximum and minimum of functions) or in evaluating velocities or forces in simulations of physical models using classical mechanics. However, the study of numerical derivation offers the possibility of introducing some pitfall in numerical methods, and it is therefore used to start this series of lab sessions on numerical methods.

The derivative in the point x of a function is defined as:

(1)

To calculate the numerical derivative as in the equation (1), we need first to discretize the derivative as the incremental ratio of single-variable function f(x) on a grid of points equally spaced at the interval

(2)

The function f(x) is known on an equally spaced lattice of x values (see Figure 1) f(x_n) with x_n=nh and n=0,±1,±2,…, and h is the step size.

Depending on the order and the point on the grid that we use for the approximate formula, we can use three different simple formulas to perform numerical derivation:

• Forward difference

• Backward difference

• Center difference

Equation (2) corresponds to the forward formula. The error in estimating the exact (analytical) derivative is proportional to the step size (h). In particular, in the forward (as in the backward difference method), the error is on the square of h. The last formula, centred or symmetric difference approximation, gives a more accurate derivative. Information on the derivation of these formulas can be found in the numerical analysis textbook (see the References). The FORTRAN program implements the centre difference formula and uses FUNCTION and SUBROUTINE statements.

PROGRAM CALCDERIV

    !
    ! Program for the calculation of derivative 
    ! using the central formula
    !
    ! (c) Roccatano

    IMPLICIT NONE

    real :: x, h, df, PI
    real :: fn
    real, dimension(:), allocatable :: fx
    integer :: i,NP

    x=0.0

    ! Number of sampling points

    NP=100

    PI=4*ATAN2(1.0,1.0)
    allocate (fx(NP))

    ! Sampling step 

    h=2*PI/float(NP)

    !
    ! Calculate and print the derivative of the function FN
    !

    do i=1,NP
        CALL DERIVATE(x,h,DF)
        PRINT 20,x,FN(x),DF
        x=x+h
    end do

    20      FORMAT (E12.5,2x,F12.5,5X,E15.8)

END PROGRAM CALCDERIV

FUNCTION FN(x)

    ! Define the function FN     

    IMPLICIT NONE
    REAL :: FN
    REAL, INTENT (IN) :: x
    FN=cos(x)
    RETURN
END FUNCTION FN

SUBROUTINE DERIVATE(x,h,DF)

    !
    ! Calculate the numerical derivative of FN
    ! using the central formula
    !

    REAL, INTENT(IN) :: x,h
    REAL, INTENT(OUT) :: DF
    DF=(FN(X+H)-FN(X-H))/(2*H)
    RETURN
END SUBROUTINE DERIVATE

The program calculate the derivative of cos(x) using 100 points. The numerical data are plotted in Figure 2 using the program xmgrace.

REFERENCES

• Etter. Structured Fortran 77 for engineers and scientists, The benjamin/ Cummings Publishing Company.
• Burden, R.L., Faires, J.D. and Burden, A.M., 2015. Numerical analysis. Cengage learning.