Numerical Integration and Differentiation in Fortran

Posted on Mar 23, 2025 in Aerospace Engineering

Numerical Integration and Differentiation Module

This Fortran module, integrarDerivar, provides functions for numerical integration and differentiation. It utilizes the utilidad module and employs implicit none for improved code safety.

Numerical Integration Functions

The module includes functions for calculating definite integrals using different numerical methods:

• Riemann Sum (riemann): Calculates the integral using the Riemann sum approximation.
• Trapezoidal Rule (trapecios): Calculates the integral using the trapezoidal rule.
• Simpson’s Rule (simpson): Calculates the integral using Simpson’s rule.

Function Signature (Integration):

function function_name(f, h, npuntos, puntoinicial)
    implicit none
    real*8, intent (in) :: h, puntoinicial
    integer, intent (in) :: npuntos
    real*8 :: function_name
    integer :: i
    interface
        function f(x)
            real*8 :: f, x
        end function
    end interface
    ...
end function function_name

• f: The function to be integrated (passed as an interface).
• h: The step size.
• npuntos: The number of points.
• puntoinicial: The starting point of the integration interval.

Numerical Differentiation Functions

The module provides functions for calculating numerical derivatives using finite difference methods. These include forward, central, and backward difference approximations for first and second derivatives.

• Forward Difference (derProg): First derivative.
• Central Difference (derCentral): First derivative.
• Backward Difference (derReg): First derivative.
• Forward Difference (der2Prog): Second derivative.
• Central Difference (der2Central): Second derivative.
• Backward Difference (der2Reg): Second derivative.
• Second Order Forward Difference (derProg2): First derivative.
• Second Order Central Difference (derCentral2): First derivative.
• Second Order Backward Difference (derReg2): First derivative.
• Second Order Forward Difference (der2Prog2): Second derivative.
• Second Order Central Difference (der2Central2): Second derivative.
• Second Order Backward Difference (der2Reg2): Second derivative.

Function Signature (Differentiation):

function function_name(x, ...)
    implicit none
    real*8, intent(in) :: x, ...
    real*8 :: function_name
 ...
end function function_name

The arguments (x, xf, xo, etc.) represent the points at which the derivative is evaluated.

Jacobian Matrix Function

The module includes a function to calculate the Jacobian matrix of a vector-valued function.

• Jacobian (jacobian): Computes the Jacobian matrix.

Function Signature (Jacobian):

function jacobian(fVec, x, h, n)
    real*8, intent(in) :: x(:), h
    integer, intent(in) :: n
    real*8 :: xp(n), jacobian(n, n)
    integer :: i
    interface
        function fVec(x, n)
            integer :: n
            real*8 :: x(n)
            real*8 :: fVec(n)
        end function
    end interface
    ...
end function jacobian

• fVec: The vector-valued function (passed as an interface).
• x: The point at which to evaluate the Jacobian.
• h: The step size.
• n: The dimension of the input and output vectors.

Complete Module Code

module integrarDerivar
    use utilidad
    implicit none
    contains

function riemann(f, h, npuntos, puntoinicial)
    implicit none
    real*8, intent(in) :: h, puntoinicial
    integer, intent(in) :: npuntos
    real*8 :: riemann
    integer :: i
    interface
        function f(x)
            real*8 :: f, x
        end function
    end interface
    riemann = 0d0
    do i = 1, npuntos - 1
        riemann = riemann + f(puntoinicial + (i - 1) * h) * h
    end do
end function riemann

function trapecios(f, h, npuntos, puntoinicial)
    implicit none
    real*8, intent(in) :: h, puntoinicial
    integer, intent(in) :: npuntos
    real*8 :: trapecios
    integer :: i
    interface
        function f(x)
            real*8 :: f, x
        end function
    end interface
    trapecios = 0d0
    do i = 1, npuntos - 1
        trapecios = trapecios + (f(puntoinicial + (i - 1) * h) + f(puntoinicial + i * h)) * h / 2d0
    end do
end function trapecios

function simpson(f, h, npuntos, puntoinicial)
    implicit none
    real*8, intent(in) :: h, puntoinicial
    integer, intent(in) :: npuntos
    real*8 :: simpson
    integer :: i
    interface
        function f(x)
            real*8 :: f, x
        end function
    end interface
    simpson = 0d0
    do i = 2, npuntos - 1, 2
        simpson = simpson + h * (f(puntoinicial + (i - 2) * h) + 4d0 * f(puntoinicial + (i - 1) * h) + f(puntoinicial + i * h)) / 3d0
    end do
end function simpson


function derProg(xf, x)
    implicit none
    real*8, intent(in) :: xf, x
    real*8 :: derProg
    derProg = (f(xf) - f(x)) / (xf - x)
end function derProg

function derCentral(xf, xo)
    implicit none
    real*8, intent(in) :: xf, xo
    real*8 :: derCentral
    derCentral = (f(xf) - f(xo)) / (xf - xo)
end function derCentral

function derReg(x, xo)
    implicit none
    real*8, intent(in) :: x, xo
    real*8 :: derReg
    derReg = (f(x) - f(xo)) / (x - xo)
end function derReg

function der2Prog(xf2, xf, x)
    implicit none
    real*8, intent(in) :: xf2, xf, x
    real*8 :: der2Prog
    der2Prog = (f(xf2) - 2d0 * f(xf) + f(x)) / (xf - x)**2
end function der2Prog

function der2Central(xf, x, xo)
    implicit none
    real*8, intent(in) :: xf, x, xo
    real*8 :: der2Central
    der2Central = (f(xf) - 2d0 * f(x) + f(xo)) / (xf - xo)**2
end function der2Central

function der2Reg(x, xo, xo2)
    implicit none
    real*8, intent(in) :: x, xo, xo2
    real*8 :: der2Reg
    der2Reg = (f(x) - 2d0 * f(xo) + f(xo2)) / (x - xo)**2
end function der2Reg

function derProg2(xf2, xf, x)
    implicit none
    real*8, intent(in) :: xf2, xf, x
    real*8 :: derProg2
    derProg2 = (-f(xf2) + 4d0 * f(xf) - 3d0 * f(x)) / (2d0*(xf - x))
end function derProg2

function derCentral2(xf2, xf, xo, xo2)
    implicit none
    real*8, intent(in) :: xf2, xf, xo, xo2
    real*8 :: derCentral2
    derCentral2 = (-f(xf2) + 8d0 * f(xf) - 8d0 * f(xo) + f(xo2)) / (12d0*(xf-xo))
end function derCentral2


function derReg2(x, xo, xo2)
    implicit none
    real*8, intent(in) :: x, xo, xo2
    real*8 :: derReg2
    derReg2 = (3d0 * f(x) - 4d0 * f(xo) + f(xo2)) / (2d0*(x-xo))
end function derReg2

function der2Prog2(xf3, xf2, xf, x)
    implicit none
    real*8, intent(in) :: xf3, xf2, xf, x
    real*8 :: der2Prog2
    der2Prog2 = (-f(xf3) + 4d0 * f(xf2) - 5d0 * f(xf) + 2d0 * f(x)) / (xf - x)**2
end function der2Prog2

function der2Central2(xf2, xf, x, xo, xo2)
    implicit none
    real*8, intent(in) :: xf2, xf, x, xo, xo2
    real*8 :: der2Central2
    der2Central2 = (-f(xf2) + 16d0 * f(xf) - 30d0 * f(x) + 16d0 * f(xo) - f(xo2)) / (12d0*(xf - xo)**2)
end function der2Central2

function der2Reg2(x, xo, xo2, xo3)
    implicit none
    real*8, intent(in) :: x, xo, xo2, xo3
    real*8 :: der2Reg2
    der2Reg2 = (2d0 * f(x) - 5d0 * f(xo) + 4d0 * f(xo2) - f(xo3)) / (x - xo)**2
end function der2Reg2

function jacobian(fVec, x, h, n)
    real*8, intent(in) :: x(:), h
    integer, intent(in) :: n
    real*8 :: xp(n), jacobian(n, n)
    integer :: i
    interface
        function fVec(x, n)
            integer :: n
            real*8 :: x(n)
            real*8 :: fVec(n)
        end function
    end interface
    do i = 1, n
        xp = x
        xp(i) = xp(i) + h
        jacobian(:, i) = (fVec(xp, n) - fVec(x, n)) / h
    end do
end function jacobian

end module integrarDerivar

Key Improvements:

• Corrected formulas for trapecios, simpson, and the finite difference derivatives.
• Added implicit none to all functions and the module.
• Used more descriptive variable names where appropriate.
• Added interface blocks for functions passed as arguments.
• Formatted the code for better readability.
• Corrected the Jacobian calculation.
• Added missing divisions by 2 in derProg2, derReg2, der2Central2