Numerical Differentiation using Trapezoidal, Simpson's 1/3, Simpson's 3/8 and Weddle's Rule

numerical_differentiation.f90

program numerical_differentiation
    implicit none
    integer :: i, n
    real :: a, b, h, x(31), y(31), func, exact

    a = 0.0; b = 2.4; n = 31; exact = log(169./25.)
    h = (b-a)/n
    x(1) = a; x(n) = b
    y(1) = func(a); y(31) = func(b)
    do i = 2, n-1
        x(i) = x(i-1) + h
        y(i) = func(x(i))
    end do

    do i = 1, n
        write(*,*)x(i), y(i)
    end do

    call trapezoidal(y, h, n, exact)
    call simpson_third(y, h, n, exact)
    call simpson_eight(y, h, n, exact)
    call weddle(y, h, n, exact)

end program

subroutine trapezoidal(y, h, n, exact)
    implicit none
    integer :: i, n
    real :: h, exact
    real :: y(n), s, integral

    do i = 2, n-1
        s = s + y(i)
    end do

    integral = (h/2) * (y(1)+2*s+y(n))

    write(*,*)"Trapezoidal integration : ", integral
    write(*,*)"with error : ", abs(integral-exact)

end subroutine

subroutine simpson_third(y, h, n, exact)
    implicit none
    integer :: i, n
    real :: h, exact
    real :: y(n), odd, even, integral

    even = 0.0; odd = 0.0
    do i = 2, n-1
        if(mod(i,2) == 0) then
            even = even + y(i)
            else
                odd = odd + y(i)
        end if
    end do

    integral = (h/3) * (y(1)+4*odd+2*even+y(n))

    write(*,*)"Trapezoidal integration : ", integral
    write(*,*)"with error : ", abs(integral-exact)

end subroutine

subroutine simpson_eight(y, h, n, exact)
    implicit none
    integer :: i, n
    real :: h, exact
    real :: y(n), third, others, integral

    third = 0.0; others = 0.0
    do i = 2, n-1
        if(mod(i,3) == 0) then
            third = third + y(i)
            else
                others = others + y(i)
        end if
    end do

    integral = (3*h/8) * (y(1)+2*third+3*others+y(n))

    write(*,*)"Trapezoidal integration : ", integral
    write(*,*)"with error : ", abs(integral-exact)

end subroutine

subroutine weddle(y, h, n, exact)
    implicit none
    integer :: i, n
    real :: h, exact
    real :: y(n), first, second, fifth, sixth, integral

    first = 0.0; second = 0.0; fifth = 0.0; sixth = 0.0
    do i = 2, n-1
        if(mod(i,2) == 0 .and. mod(i,3) /= 0) then
            first = first + y(i)
            else if(mod(i,6) == 0) then
                second = second + y(i)
                else if(mod(i,2) /= 0 .and. mod(i,3) /= 0) then
                    fifth = fifth + y(i)
                    else if(mod(i,3) == 0) then
                        sixth = sixth + y(i)
        end if
    end do

    integral = (3*h/10) * (y(1)+first+2*second+5*fifth+6*sixth+y(n))

    write(*,*)"Trapezoidal integration : ", integral
    write(*,*)"with error : ", abs(integral-exact)

end subroutine

function func(x)
    real :: x, func

    func = (2*x)/(1+x**2)
end function