Numerical Integration (Trapezoidal / Simpson)

! -----------------------------------------------
! Trapezoidal Rule
! ∫[a,b] f(x) dx ≈ h/2 * (f(a) + 2*f(x1) + 2*f(x2) + ... + f(b))
! h = (b - a) / n  (n: number of subdivisions)
! -----------------------------------------------

real(kind=8) :: a, b, h, result
integer      :: n, i

n = 1000          ! Number of subdivisions (higher = more accurate)
h = (b - a) / n

result = 0.5_8 * (f(a) + f(b))
do i = 1, n - 1
    result = result + f(a + i * h)
end do
result = result * h

! -----------------------------------------------
! Simpson's Rule
! n must be even
! ∫[a,b] f(x) dx ≈ h/3 * (f(x0) + 4*f(x1) + 2*f(x2) + 4*f(x3) + ... + f(xn))
! -----------------------------------------------

n = 1000          ! Number of subdivisions (must be even)
h = (b - a) / n

result = f(a) + f(b)
do i = 1, n - 1
    if (mod(i, 2) == 1) then
        result = result + 4.0_8 * f(a + i * h)   ! Odd index: coefficient 4
    else
        result = result + 2.0_8 * f(a + i * h)   ! Even index: coefficient 2
    end if
end do
result = result * h / 3.0_8

jjk_integration.f90

! jjk_integration.f90 — Numerical integration using the trapezoidal and Simpson's rules
! Integrates a function modeled after a cursed-energy curve from Jujutsu Kaisen
!
! Compile and run:
!   gfortran jjk_integration.f90 -o jjk_integration && ./jjk_integration

program jjk_integration
    implicit none

    integer, parameter :: DP = selected_real_kind(15, 307)

    ! -----------------------------------------------
    ! Integration range and number of subdivisions
    ! Integrand: f(x) = x^2 * sin(x)  [0, π]
    ! Analytical solution: [-x^2*cos(x) + 2*x*sin(x) + 2*cos(x)] evaluated on [0,π] = π^2 - 4
    ! -----------------------------------------------
    real(DP), parameter :: PI     = 3.14159265358979323846_DP
    real(DP), parameter :: A      = 0.0_DP   ! Lower bound
    real(DP), parameter :: B      = PI       ! Upper bound
    real(DP), parameter :: EXACT  = PI**2 - 4.0_DP  ! Analytical solution

    integer  :: n
    real(DP) :: trap_result, simp_result
    real(DP) :: trap_err, simp_err

    ! -----------------------------------------------
    ! Compare accuracy at each number of subdivisions
    ! -----------------------------------------------
    print *, "===== Jujutsu Kaisen Cursed Energy Integration ====="
    print *, "  f(x) = x^2 * sin(x)  [0, π]"
    print "(a, f18.15)", "  Analytical solution (π^2 - 4) = ", EXACT
    print *, ""
    print "(a)", "  Subdivisions   Trapezoidal          Error (Trap)      Simpson's            Error (Simp)"

    do n = 10, 1000, 10
        ! Integration using the trapezoidal rule
        trap_result = trapezoidal(A, B, n)
        trap_err    = abs(trap_result - EXACT)

        ! Integration using Simpson's rule (n must be even)
        simp_result = simpsons(A, B, n)
        simp_err    = abs(simp_result - EXACT)

        if (n == 10 .or. n == 100 .or. n == 1000) then
            print "(i8, f20.15, es14.4, f20.15, es14.4)", &
                n, trap_result, trap_err, simp_result, simp_err
        end if
    end do

    print *, ""

    ! -----------------------------------------------
    ! Gojo's Infinite Cursed Energy (convergence check)
    ! Doubles the subdivisions until convergence
    ! -----------------------------------------------
    print *, "===== Gojo's Infinite Cursed Energy — Convergence Check ====="
    print *, "  Doubling subdivisions with Simpson's rule until convergence."
    print *, ""

    block
        real(DP) :: prev, curr, tol
        integer  :: nn

        tol  = 1.0e-12_DP   ! Tolerance
        nn   = 2
        prev = simpsons(A, B, nn)

        do
            nn   = nn * 2
            curr = simpsons(A, B, nn)

            print "(a, i6, a, es10.3, a, f18.15)", &
                "  n=", nn, "  diff=", abs(curr - prev), "  result=", curr

            if (abs(curr - prev) < tol) exit
            prev = curr
            if (nn > 100000) exit   ! Safety guard
        end do

        print *, ""
        print "(a, f18.15)", "  Converged value : ", curr
        print "(a, f18.15)", "  Analytical      : ", EXACT
        print "(a, es10.3)", "  Final error     : ", abs(curr - EXACT)
    end block

    print *, ""
    print *, "  Cursed energy integration for Itadori, Fushiguro, and Kugisaki complete."

contains

    ! -----------------------------------------------
    ! Integrand: f(x) = x^2 * sin(x)
    ! -----------------------------------------------
    pure real(DP) function f(x)
        real(DP), intent(in) :: x
        f = x**2 * sin(x)
    end function f

    ! -----------------------------------------------
    ! Trapezoidal rule: integrates over [a, b] with n subdivisions
    ! -----------------------------------------------
    real(DP) function trapezoidal(a, b, n)
        real(DP), intent(in) :: a, b
        integer,  intent(in) :: n
        real(DP) :: h, s
        integer  :: i

        h = (b - a) / n
        s = 0.5_DP * (f(a) + f(b))   ! Endpoints use coefficient 0.5

        do i = 1, n - 1
            s = s + f(a + i * h)
        end do
        trapezoidal = s * h
    end function trapezoidal

    ! -----------------------------------------------
    ! Simpson's rule: integrates over [a, b] with n subdivisions
    ! n must be even
    ! -----------------------------------------------
    real(DP) function simpsons(a, b, n)
        real(DP), intent(in) :: a, b
        integer,  intent(in) :: n
        real(DP) :: h, s
        integer  :: i, nn

        nn = n
        if (mod(nn, 2) /= 0) nn = nn + 1   ! Increment to even if odd

        h = (b - a) / nn
        s = f(a) + f(b)   ! Endpoints use coefficient 1

        do i = 1, nn - 1
            if (mod(i, 2) == 1) then
                s = s + 4.0_DP * f(a + i * h)   ! Odd index: coefficient 4
            else
                s = s + 2.0_DP * f(a + i * h)   ! Even index: coefficient 2
            end if
        end do
        simpsons = s * h / 3.0_DP
    end function simpsons

end program jjk_integration

gfortran jjk_integration.f90 -o jjk_integration && ./jjk_integration
===== Jujutsu Kaisen Cursed Energy Integration =====
  f(x) = x^2 * sin(x)  [0, π]
  Analytical solution (π^2 - 4) =  5.869604401089358

  Subdivisions   Trapezoidal          Error (Trap)      Simpson's            Error (Simp)
      10   5.869199919916741  4.0481E-04   5.869604401167741  7.8383E-11
     100   5.869600349476093  4.0516E-06   5.869604401089358  7.1054E-15
    1000   5.869604400684077  4.0528E-10   5.869604401089358  0.0000E+00

===== Gojo's Infinite Cursed Energy — Convergence Check =====
  Doubling subdivisions with Simpson's rule until convergence.

  n=     4  diff= 1.268E-08  result=  5.869604401097723
  n=     8  diff= 5.405E-11  result=  5.869604401089358
  n=    16  diff= 0.000E+00  result=  5.869604401089358

  Converged value :  5.869604401089358
  Analytical      :  5.869604401089358
  Final error     :  0.000E+00

  Cursed energy integration for Itadori, Fushiguro, and Kugisaki complete.

The trapezoidal rule approximates each subinterval with a trapezoid, yielding an error of O(h²) where h is the step size. Simpson's rule combines two adjacent subintervals into a single parabolic approximation, reducing the error to O(h⁴) and significantly improving accuracy with the same number of subdivisions. However, Simpson's rule requires the number of subdivisions to be even. A common approach for convergence testing is to compare the absolute difference between results for n and 2n subdivisions against a tolerance threshold. For even higher accuracy, consider Gaussian quadrature. Defining the integrand as an internal subprogram inside a contains section allows it to be referenced naturally from the main program and improves maintainability. For details on floating-point precision limits, see Floating-Point Precision and Rounding Errors.