## Numerical integration

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
jdebord
Posts: 550
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Numerical integration

Here is my FreeBASIC version of an integration program in C by Robert van Engelen. It is posted here with permission of the original author.

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
``````
hhr
Posts: 216
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

Very nice, is there a way to estimate the error?

Code: Select all

``````function f4(x as double) as double
return 1-sqr(x*x-1)/x
end function

printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY))
printf(!"should be                   (pi-2)/2 = %.15g\n", M_PI/2 - 1)
``````

Code: Select all

``````dim shared as double n

function f5(x as double) as double
return (x^n)*exp(-x)
end function

for n = 0 to 4
print "n =";n
print "n! =  ";tgamma(n+1)
print string(25, "-")
next
``````
jdebord
Posts: 550
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: Numerical integration

The quad function has a sixth parameter err_est with a default value of 0. To return it, initialize it to a different value.

Code: Select all

``````dim as double err_est = 1

function f4(x as double) as double
return 1-sqr(x*x-1)/x
end function

printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY, , , err_est))
printf(!"should be                   (pi-2)/2 = %.15g\n", M_PI/2 - 1)

print "estimated error = ", err_est
``````
Result:

Code: Select all

``````integrate(1-sqr(x*x-1)/x, x=1..+inf) = 0.570796326697611
should be                   (pi-2)/2 = 0.570796326794897
estimated error =            2.962455091046144e-011
``````
dodicat
Posts: 7993
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Numerical integration

Using Simpson's rules (as in shipyards for waterplane areas) gets quite close.
...
COMPARE:

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

function f4(x as double) as double
return 1-sqr(x*x-1)/x
end function

function f5(x as double) as double
return x^2
end function

function f6(x as double) as double
return 1/(1+exp(x))
end function

'=================================================
#macro SimpsonArea(fn,Lo,Up,ordinates,sum)
scope
dim as double x,part1,part2
dim as double l=lo,u=up,n=iif(ordinates mod 2=0,ordinates+1,ordinates)
var z=(len(str(n))-1)/2
if abs(l)=infinity then l=sgn(l)*10^4
if abs(u)=infinity then u=sgn(u)*10^4
Var h=(U-L)/n
x=L+.5*h:Part1=fn
For i As long=1 To n-1
x=L+h*i+h/2
Part1+=fn
x=L+h*i
Part2+=fn
Next i
x=L:var fn1=fn
x=U:var fn2=fn
sum= (h/6)*(fn1+fn2+Part1*4+Part2*2)
end scope
#endmacro
'==================================================

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))
printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", quad(@f4, 1, INFINITY))
printf(!"integrate(1/(1+exp(x)), x=0..+inf)   = %.15g\n", quad(@f6, 0, INFINITY))
printf(!"integrate(x^2, x=-100.005..12.556)   = %.15g\n",quad(@f5, -100.005, 120.556))

print
'Simpson macro
dim as double Result
const ordinates=1000000 ' no less anyway.
SimpsonArea(acos(x),0,1,ordinates,Result)
printf(!"integrate(acos(x), x=0..1)           = %.15g\n", Result)

SimpsonArea(exp(-x/5),0,infinity,ordinates,Result)
printf(!"integrate(exp(-x/5), x=0..+inf)      = %.15g\n", Result)

SimpsonArea(1/cosh(x)^2,-infinity,infinity,ordinates,Result)
printf(!"integrate(1/cosh(x)^2, x=-inf..+inf) = %.15g\n", Result)

SimpsonArea(1-sqr(x*x-1)/x,1,infinity,ordinates,Result)
printf(!"integrate(1-sqr(x*x-1)/x, x=1..+inf) = %.15g\n", Result)

SimpsonArea(1/(1+exp(x)),0,infinity,ordinates,Result)
printf(!"integrate(1/(1+exp(x)), x=0..+inf)   = %.15g\n", Result)

SimpsonArea(x^2,-100.005,120.556,ordinates,Result)
printf(!"integrate(x^2, x=-100.005..120.556)  = %.15g\n", Result)
print
print "actual value for last integral (x^2)  ";120.556^3/3- -100.005^3/3
print

sleep
``````