Code: Select all
' ************************************************************************
' Numerical integration by Tanh-Sinh, Sinh-Sinh and Exp-Sinh formulas
' ************************************************************************
' Ref. : R. A. van Engelen, https://www.genivia.com/files/qthsh.pdf
' ************************************************************************
#include "crt.bi"
#define isfinite(x) (abs(x) <> INFINITY)
function exp_sinh_opt_d(func as function (x as double) as double, _
a as double, eps as double, d as double) as double
dim as double fl, fr, h, h2, lfl, lfr, lr, r, s
dim as long ev, i, j
ev = 2 : i = 1 : j = 32 ' j=32 is optimal to search for r
h2 = func(a + d / 2) - func(a + d * 2) * 4
if (isfinite(h2) and abs(h2) > 1e-5) then ' if |h2| > 2^-16
s = 0
lr = 2
do
' find max j such that fl and fr are finite
j /= 2
r = 1 shl (i + j)
fl = func(a + d / r)
fr = func(a + d * r) * r * r
ev += 2
h = fl - fr
loop while (j > 1 and not(isfinite(h)))
if (j > 1 and isfinite(h) and sgn(h) <> sgn(h2)) then
lfl = fl ' last fl=func(a+d/r)
lfr = fr ' last fr=func(a+d*r)*r*r
do ' bisect in 4 iterations
j /= 2
r = 1 shl (i + j)
fl = func(a + d / r)
fr = func(a + d * r) * r * r
ev += 2
h = fl - fr
if (isfinite(h)) then
s += abs(h) ' sum |h| to remove noisy cases
if (sgn(h) = sgn(h2)) then
i += j ' search right half
else ' search left half
lfl = fl ' record last fl=f(a+d/r)
lfr = fr ' record last fr=f(a+d*r)*r*r
lr = r ' record last r
end if
end if
loop while (j > 1)
if (s > eps) then ' if sum of |h| > eps
h = lfl - lfr ' use last fl and fr before the sign change
r = lr ' use last r before the sign change
' if last difference was nonzero, back up r by one step
if (h <> 0) then r /= 2
if (abs(lfl) < abs(lfr)) then
d /= r ' move d closer to the finite endpoint
else
d *= r ' move d closer to the infinite endpoint
end if
end if
end if
end if
return d
end function
function quad (func as function(x as double) as double, a as double, b as double, _
n as long = 6, eps as double = 1.0e-9, byref err_est as double = 0) as double
dim as double c, d, eh, fp, fm, h, p, q, r, s, sign, t, tol, u, v, w, x, y
dim as long k
dim as long mode ' Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
d = 1 : sign = 1 : h = 2 : tol = 10 * eps
if (b < a) then swap a, b : sign = -1
if (isfinite(a) and isfinite(b)) then
c = (a + b) / 2
d = (b - a) / 2
v = c
elseif (isfinite(a)) then
mode = 1 ' Exp-Sinh
d = exp_sinh_opt_d(func, a, eps, d)
c = a
v = a + d
elseif (isfinite(b)) then
mode = 1 ' Exp-Sinh
d = exp_sinh_opt_d(func, b, eps, -d)
sign = -sign
c = b
v = b + d
else
mode = 2 ' Sinh-Sinh
v = 0
end if
s = func(v)
do
p = 0 : fp = 0 : fm = 0
h /= 2
eh = exp(h)
t = eh
if (k > 0) then eh *= eh
if (mode = 0) then ' Tanh-Sinh
do
u = exp(1 / t - t) ' = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
r = 2 * u / (1 + u) ' = 1 - tanh(sinh(j*h))
w = (t + 1 / t) * r / (1 + u) ' = cosh(j*h)/cosh(sinh(j*h))^2
x = d * r
if (a + x > a) then ' if too close to a then reuse previous fp
y = func(a + x)
if (isfinite(y)) then fp = y ' if f(x) is finite, add to local sum
end if
if (b - x < b) then ' if too close to b then reuse previous fm
y = func(b - x)
if (isfinite(y)) then fm = y ' if f(x) is finite, add to local sum
end if
q = w * (fp + fm)
p += q
t *= eh
loop while (abs(q) > eps * abs(p))
else
t /= 2
do
r = exp(t - 0.25 / t) ' = exp(sinh(j*h))
w = r
q = 0
if (mode = 1) then ' Exp-Sinh
x = c + d / r
if (x = c) then exit do ' if x hit the finite endpoint then break
y = func(x)
if isfinite(y) then q += y / w ' if f(x) is finite, add to local sum
else ' Sinh-Sinh
r = (r - 1 / r) / 2 ' = sinh(sinh(j*h))
w = (w + 1 / w) / 2 ' = cosh(sinh(j*h))
x = c - d * r
y = func(x)
if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
end if
x = c + d * r
y = func(x)
if isfinite(y) then q += y * w ' if f(x) is finite, add to local sum
q *= t + 0.25 / t ' q *= cosh(j*h)
p += q
t *= eh
loop while (abs(q) > eps * abs(p))
end if
v = s - p
s += p
k += 1
loop while (abs(v) > tol * abs(s) and k <= n)
' if the estimated relative error is desired, then return it
if err_est <> 0 then err_est = abs(v) / (abs(s) + eps)
' result with estimated relative error err_est
return sign * d * s * h
end function
' ************************************************************************
' example functions to integrate, exact integrals are 1, 5 and 2
function f1(x as double) as double
return acos(x)
end function
function f2(x as double) as double
return exp(- x / 5)
end function
function f3(x as double) as double
return 1 / cosh(x)^2
end function
printf(!"integrate(acos(x), x=0..1) = %.15g\n", quad(@f1, 0, 1))
printf(!"integrate(exp(-x/5), x=0..+inf) = %.15g\n", quad(@f2, 0, INFINITY))
printf(!"integrate(1/cosh(x)^2, x=-inf..+inf) = %.15g\n", quad(@f3, -INFINITY, INFINITY))
sleep