Numerical integration/Adaptive Simpson's method: Difference between revisions

=={{header|Fortran}}==

<syntaxhighlight lang="Fortran">

program adaptive_simpson_task

implicit none

integer, parameter :: dbl = kind (1.0d0)

interface

function function_dbl_to_dbl (x) result (y)

import dbl

real(dbl), intent(in) :: x

real(dbl) :: y

end function function_dbl_to_dbl

end interface

write (*, 1000) quad_asr (sine, 0.0_dbl, 1.0_dbl, 1E-9_dbl, 1000)

1000 format ("estimated definite integral of sin(x) ", &

& "for x from 0 to 1: ", F15.13)

contains

function sine (x) result (y)

real(dbl), intent(in) :: x

real(dbl) :: y

y = sin (x)

end function sine

function quad_asr (f, a, b, tol, depth) result (quadval)

procedure(function_dbl_to_dbl) :: f

real(dbl), intent(in) :: a, b, tol

integer, intent(in) :: depth

real(dbl) :: quadval

real(dbl) :: fa, fb, m, fm, whole

fa = f(a)

fb = f(b)

call quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, whole)

quadval = quad_asr_ (f, a, fa, b, fb, tol, whole, m, fm, depth)

end function quad_asr

recursive function quad_asr_ (f, a, fa, b, fb, tol, whole, &

& m, fm, depth) result (quadval)

procedure(function_dbl_to_dbl) :: f

real(dbl), intent(in) :: a, fa, b, fb, tol, whole, m, fm

integer, intent(in) :: depth

real(dbl) :: quadval

real(dbl) :: lm, flm, left

real(dbl) :: rm, frm, right

real(dbl) :: delta, tol_

call quad_asr_simpsons_ (f, a, fa, m, fm, lm, flm, left)

call quad_asr_simpsons_ (f, m, fm, b, fb, rm, frm, right)

delta = left + right - whole

tol_ = tol / 2

if (depth <= 0 .or. tol_ == tol .or. abs (delta) <= 15 * tol) then

quadval = left + right + (delta / 15)

else

quadval = quad_asr_ (f, a, fa, m, fm, tol_, left, &

& lm, flm, depth - 1)

quadval = quadval + quad_asr_ (f, m, fm, b, fb, tol_, &

& right, rm, frm, depth - 1)

end if

end function quad_asr_

subroutine quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, quadval)

procedure(function_dbl_to_dbl) :: f

real(dbl), intent(in) :: a, fa, b, fb

real(dbl), intent(out) :: m, fm, quadval

m = (a + b) / 2

fm = f(m)

quadval = ((b - a) / 6) * (fa + (4 * fm) + fb)

end subroutine quad_asr_simpsons_

end program adaptive_simpson_task

</syntaxhighlight>

{{out}}

<pre>$ gfortran -O -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out

estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318</pre>

On x86-64, if you leave out the optimizer you get a trampoline:

<pre>$ gfortran -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out

adaptive_simpson_task.f90:20:2:

20 | function sine (x) result (y)

| ^

Warning: trampoline generated for nested function ‘sine’ [-Wtrampolines]

/usr/lib/gcc/x86_64-pc-linux-gnu/12/../../../../x86_64-pc-linux-gnu/bin/ld: warning: /tmp/cc2mQOD2.o: requires executable stack (because the .note.GNU-stack section is executable)

estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318</pre>

Such a bit of executable code on the stack often is not allowed on systems hardened against hackers. Otherwise you probably needn't be concerned, even though warnings are being generated both by the compiler ''and'' the linker. :)