NuMFor 9f2ab49 (2024-04-08)
Numerical (Modern) Fortran. Library for Simple Numerical computing
List of all members
qnthsh Interface Reference

Subroutine qnthsh implements integration by tanh-sinh method. More...

Detailed Description

Subroutine qnthsh implements integration by tanh-sinh method.

This routine is a simple non-adaptive automatic integrator using trapezoidal rule to integrate a the transformed integrand using the double exponential (also called tanh-sinh) method. The tahn-sinh scheme is based on the observation that the trapezoid rule converges very rapidly for functions on the entire real line that go to zero like exp( - exp(t) ). The change of variables $x = \tanh( \pi \sinh(t) /2)$ transforms an integral over [-1, 1] into an integral with integrand suited to the double exponential rule.

The transformed integral is infinite, but we truncate the domain of integration to [-a, a] ( $a \le 3$). The value of $a$ is estimated from the required precision, and the value 3 was chosen for two reasons: for $ t = 3$, the transformed $x$ values are nearly equal to 1 (up to 12 significant figures). Also, for $ t = 3$, the smallest weights are 12 orders of magnitude smaller than the largest weights.

The integration first applies the trapezoid rule to $[-a, a]$. Then it subsequently cuts the step size in half each time, comparing the results. Integration stops when subsequent iterations are close enough together or the maximum integration points have been used. By cutting $h$ in half, the previous integral can be reused; we only need evaluate the integrand at the newly added points.

References:

Parameters
[in]fa real or complex function (as defined in Integrable functions)
[in]a(real(dp)) are the limits of integration
[in]b(real(dp)) are the limits of integration
[out]IntValue(real or complex, depending on f) Estimated value of the integral
[in]epsabs(real, optional) absolute desired precision. Default: epsabs=1.e-7
[in]epsrel(real, optional) relative desired precision. Default: epsrel=1.e-5
[out]abserr(real, optional) absolute estimated error
[out]Neval(integer, optional) Number of function evaluations

Example:

real(dp) :: Integ1
intrinsic dsin
print "(A)", center(" Integrate sin(x) between 0 and pi ", 70, '-')
call qnthsh(dsin, zero, m_pi, integ1)
print "(A)", '\int \sin(x) dx = '//str(integ1)//" (Difference="//str(abs(2 - integ1))//")"
! that prints:
!------------------ Integrate sin(x) between 0 and pi -----------------
! \int \sin(x) dx = 2 (Difference=0)
The documentation for this interface was generated from the following file: