Numerical integration in Fortran 90

Question

In Fortran 90, I want to numerically integrate a mathematical function with one variable within a given limit. For example, integrating f(x) = x**2 from 0 to 10. The function I have is more complicated than this one and I have to run it several times changing the integration limits. I found out on internet that the 'QUADPACK' library might help me with this. But how can I install this library so that I can call this in my code? Provide some details as I won't be able to follow advanced instructions quickly.

https://en.wikipedia.org/wiki/Trapezoidal_rule or https://en.wikipedia.org/wiki/Simpson's_rule they are both pretty simple to implement your self — jk., Apr 01 '16 at 11:45
Please provide more information about the equation you're integrating, operating system you're using so that these instructions can be answered for your particular case. This would be more helpful. — Charles, Apr 08 '16 at 17:04
@Charlie- I would like to know this for any general equation. For example - integrating a*x**2 + b*x + c over the limit l1 to l2 where a,b, and c are some constants. I have an array of l1 and l2 i.e. I have to integrate this equation many times on an array of limits. — NewToFortran, Apr 10 '16 at 09:36
@Charlie - Is there a way to do it? In python, we can 'vectorize' an equation and feed it array limits. Is this possible in Fortran 90? Or do I have run a for loop for this? I was running Fortran 90 using Netbeans IDE and GCC compiler. I am soon shifting to Linux because I found that more material as well as library supposrt is available for Linux based systems when it comes to Fortran. So do you have a solution for this? Thank you! — NewToFortran, Apr 10 '16 at 09:36

Charles · Answer 1 · 2016-04-11T15:26:50.733

0

I've provided a simple program where midpoint method is used to integrate x^2. A more complicated formula may be entered, so long the mesh is fine enough (and the function is smooth), this should work..

    program integrate
    implicit none
    integer,parameter :: cp = selected_real_kind(14)
    integer,parameter :: N = 1000
    real(cp),dimension(N) :: f,xc
    real(cp),dimension(N+1) :: x
    real(cp) :: s,xmax,xmin,dx
    integer :: i
    xmin = 0.0_cp
    xmax = 10.0_cp
    dx = (xmax - xmin)/real(N,cp)
    x = (/(xmin + dx*(i-1),i=1,N+1)/)
    ! Define x at center
    do i=1,N
    xc(i) = x(i) + 0.5_cp*dx
    enddo
    ! Define f
    do i=1,N
    f(i) = xc(i)**2
    enddo
    ! Integrate (Midpoint method)
    s = 0.0_cp
    do i=1,N
    s = s + f(i)*dx
    enddo
    write(*,*) 'sum = ',s
    end program

edited Apr 11 '16 at 15:26

answered Apr 10 '16 at 20:06

Charles

947
1
15
39

I've used x**2.0_cp as a precaution to insure that selected_real_kind is used in this definition. I've had issues in the past using simply 2.0, or even real(2.0,cp). – Charles Apr 10 '16 at 21:23
@francescalus your example needs some parenthesis.. `-1.0_cp**2` is `-1.0` .. you are correct in the point though if you mean to square something use the integer 2. – agentp Apr 11 '16 at 10:09

score 0 · Answer 2 · answered Dec 11 '22 at 13:35

Here is one possible solution to integrate your function f(x) = x**2 from 0 to 10. This uses the Gaussian quadrature formula for 8 and 16 points.

program quadrature
  implicit none

    ! Declare variables
    integer, parameter :: n8 = 8, n16 = 16
    real(8) :: r, m, c
    real(8) :: a, b, result8, result16
    real(8), dimension (n8)   :: x8, w8
    real(8), dimension(n16)   :: x16, w16
    integer :: i

    ! Define the limits of integration
    a = 0.D0
    b = 10.D0

    ! Define the abscissas and weights for 8-point Gauss quadrature
    data x8 /-0.1834346424956498D0, 0.1834346424956498D0, -0.5255324099163290D0, 0.5255324099163290D0, &
             -0.7966664774136267D0, 0.7966664774136267D0, -0.9602898564975363D0, 0.9602898564975363D0/
    data w8 / 0.3626837833783620D0, 0.3626837833783620D0, 0.3137066458778873D0, 0.3137066458778873D0, &
              0.2223810344533745D0, 0.2223810344533745D0, 0.1012285362903763D0, 0.1012285362903763D0/ 

    ! Define the abscissas and weights for 16-point Gauss quadrature
    data x16 /-0.0950125098376374D0, 0.0950125098376374D0, -0.2816035507792589D0, 0.2816035507792589D0, &  
              -0.4580167776572274D0, 0.4580167776572274D0, -0.6178762444026438D0, 0.6178762444026438D0, &  
              -0.7554044083550030D0, 0.7554044083550030D0, -0.8656312023878318D0, 0.8656312023878318D0, &  
              -0.9445750230732326D0, 0.9445750230732326D0, -0.9894009349916499D0, 0.9894009349916499D0  /  
              
    data w16 /0.1894506104550685D0, 0.1894506104550685D0, 0.1826034150449236D0, 0.1826034150449236D0, &
              0.1691565193950025D0, 0.1691565193950025D0, 0.1495959888165767D0, 0.1495959888165767D0, &
              0.1246289712555339D0, 0.1246289712555339D0, 0.0951585116824928D0, 0.0951585116824928D0, &
              0.0622535239386479D0, 0.0622535239386479D0, 0.0271524594117541D0, 0.0271524594117541D0  /

    ! Compute the results using 8-point and 16-point Gauss quadrature
    r = 0.D0
    m = (b-a)/2.D0
    c = (b+a)/2.D0
    result8 = 0.D0
    result16 = 0.D0
    do i = 1, n8
        result8 = result8 + w8(i) * f(m*x8(i) + c)
    end do
    result8 = result8*m
    do i = 1, n16
        result16 =  result16 + w16(i) * f(m*x16(i) + c)
    end do
    result16 = result16*m
    ! Print the results
    print *, "Result using 8-point Gauss quadrature: ", result8
    print *, "Result using 16-point Gauss quadrature: ", result16

contains

    ! Function to be integrated
    real(8) function f(x)
        real(8), intent(in) :: x
        f = x**2.D0
    end function

end program

The data statement is really old stuff. Fortran 90+ has better options for initialization. Also, magic numbers like 8 in `real(8)` are not portable. For exponention, it is enough to just use `x**2`, there is no need for a real exponent. In compilers that do not optimize well, the real exponent could even be much slower. — Vladimir F Героям слава, Dec 19 '22 at 06:50
Ok I have a few questions then. What would be better for initialization? Also, what do you mean by not portable? What's wrong with using magic numbers? How much slower exactly is using x**2.D0, is it really noticeable? — siimurik, Dec 21 '22 at 19:24
For initialization one uses `x8 = [ 1., 2.,... ]` on the declaration line. For kinds see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter A good optimizing compiler will be able to optimize what you wrote in the exponentiation, but `exp(2.d0 * log(x))` is much slower than just `x*x`. And it is much easier for the compiler to see the latter in just `x**2`. — Vladimir F Героям слава, Dec 21 '22 at 20:21

Numerical integration in Fortran 90

2 Answers2