## Numerical integration

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
jdebord
Posts: 551
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: 222
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: 551
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: 7998
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
``````
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

First of all, thank you for the beautiful program.

The following example works for n = 0 to 97 and then again for n = 168 to 170.
For n > 170, FreeBASIC-Double can no longer display the result.

Code: Select all

``````dim shared as double n
dim as double q, g, err_est

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

for n = 0 to 180
print "n ="; n
g = tgamma(n+1)
print "n! =  "; g
err_est = 1
q = quad(@f5, 0, INFINITY, , , err_est)
print "quad ="; q; " "; csng(err_est); " "; csng((g - q) / g)
print string(25, "-")
getkey
next
``````
If INFINITY is replaced by 1000, it works better:

Code: Select all

``````q = quad(@f5, 0, 1000, , , err_est)
``````
I am surprised that the values for INFINITY and n = 98 to 167 are too small. Is there a bug with INFINITY?

I have experimented with gas32.
With gcc and gas64 it behaves differently and with gas64 it crashes with b = 1000 at n = 103.

err_est is very useful.
SARG
Posts: 1790
Joined: May 27, 2005 7:15
Location: FRANCE

### Re: Numerical integration

hhr wrote: Jun 20, 2024 19:30 with gas64 it crashes with b = 1000 at n = 103.
I'll fix the problem in gas64 but in your example what is tgamma ?
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

tgamma is the gamma function in crt, which is declared in crt\math.bi.

https://en.wikipedia.org/wiki/Gamma_function

In the hope that I can be of some help, I suspect that the program is stuck in the loop that ends with line 141:

Code: Select all

``````loop while (abs(q) > eps * abs(p))
``````
If I switch to any other data type, it works, for example:

Code: Select all

``````loop while cbyte(abs(q) > eps * abs(p))
``````
or

Code: Select all

``````loop while cdbl(abs(q) > eps * abs(p))
``````
SARG
Posts: 1790
Joined: May 27, 2005 7:15
Location: FRANCE

### Re: Numerical integration

ok, thks.

I don't get a crash with b=1000. Using very last fbc 1.20 on Windows.

Code: Select all

``````n = 102
n! =   9.614466715035127e+161
-------------------------
n = 103
n! =   9.902900716486181e+163
-------------------------
n = 104
n! =   1.029901674514563e+166
-------------------------``````
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

Indeed, it works with fbc 1.20.0, I had used fbc 1.10.1.

Many thanks for gas64 and the effort.
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

With the latest fbc 1.20.0, gcc32, gcc64 and gas64 work exactly the same except for slight differences in the calculated numbers.
gas32 has a different behavior.

What distinguishes gas32 from the other compilers?
dodicat
Posts: 7998
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Numerical integration

32 bit -gen gas works out better for these areas.
I have added simpson's shipyard area (with 500000 iterations -which is more than any ship would get) and a string factorial. (for fun)
infinity at 10^4 for simpson
Choose your compiler option at line 1.

Code: Select all

``````'#cmdline "-gen gas64"
#include "crt.bi"
#define isfinite(x) (abs(x) <> INFINITY)
dim shared as double n
dim as double q, g, err_est

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

Function factorial(num As String) As String
static  As Ubyte _Mod(0 To 99),_Div(0 To 99),flag=0
if flag=0 then
For z As Integer=0 To 99:_Mod(z)=(z Mod 10+48):_Div(z)=z\10:next
flag=1
end if
Var fact="1",a="",b="",c=""
Dim As Ubyte n,carry,ai
For z As Integer=1 To Valint(num)
a=fact:b=Str(z):Var la =Len(a),lb =Len(b)
c =String(la+lb,"0")
For i As Integer =la-1 To 0 Step -1
carry=0:ai=a[i]-48
For j As Integer =lb-1 To 0 Step -1
n =ai*(b[j]-48)+(c[i+j+1]-48)+carry
carry =_Div(n):c[i+j+1]=_Mod(n)
Next j
c[i]+=carry
Next i
fact=Ltrim(c,"0")
Next z
Return fact
End Function

Function Simpson(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer) As Double
if abs(l)=infinity then l=sgn(l)*10^4
if abs(u)=infinity then u=sgn(u)*10^4
var n=Ordinates
If n Mod 2=0 Then n=n+1
Var h=(U-L)/n
Var Part1=fn(L+.5*h)
Var Part2=0.0
For i As Integer=1 To n-1
Part1+=fn(L+h*i+h/2)
Part2+=fn(L+h*i)
Next i
Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2)
End Function

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

for n = 0 to 180
print "n = "; n
g = tgamma(n+1)
var f=factorial(str(n))
print "n! =   "; g;" f = ";f
err_est = 1
q = quad(@f5, 0, INFINITY, , , err_est)
print "quad = "; q; " "; csng(err_est); " "; csng((g - q) / g)
var area=simpson(@f5,0,infinity,500000)
print "simpson";area;" "; " err not given"; " "; csng((val(f) - area) / val(f))
print string(25, "-")
getkey
next ``````
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

Simpson works well, it's just a bit slow.
I also find the factorial function interesting. I didn't realize that you can easily convert the string to Single or Double.

Code: Select all

``````#include "crt.bi"

Function factorial(num As String) As String 'by dodicat
Static  As Ubyte _Mod(0 To 99),_Div(0 To 99),flag=0
If flag=0 Then
For z As Integer=0 To 99
_Mod(z)=(z Mod 10+48)
_Div(z)=z\10
Next
flag=1
End If
Var fact="1",a="",b="",c=""
Dim As Ubyte n,carry,ai
For z As Integer=1 To Valint(num)
a=fact:b=Str(z):Var la =Len(a),lb =Len(b)
c =String(la+lb,"0")
For i As Integer =la-1 To 0 Step -1
carry=0:ai=a[i]-48
For j As Integer =lb-1 To 0 Step -1
n =ai*(b[j]-48)+(c[i+j+1]-48)+carry
carry =_Div(n):c[i+j+1]=_Mod(n)
Next j
c[i]+=carry
Next i
fact=Ltrim(c,"0")
Next z
Return fact
End Function

Dim As Long n
Dim As String f

For n=0 To 200
Print "n =";n
f = factorial(Str(n))
Print "f = ";f
Print "cdbl(f) ="; Cast(Double,f)
Print "tgamma = ";tgamma(n+1)
Print String(20,"-")
Getkey
Next

Sleep
``````
hhr
Posts: 222
Joined: Nov 29, 2019 10:41

### Re: Numerical integration

I wanted to replace Double with another data type and while searching I noticed this topic:

https://www.freebasic.net/forum/viewtopic.php?t=27314

I copied the bi files from the first three posts and rebuilt the example from jdebord's first post.
This gives me the same results for all four compilers.
This even works with Windows XP, but unfortunately not with Linux.

However, you have to refrain from including crt.bi.

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
' ************************************************************************

'#cmdline "-gen gcc"

Print __fb_signature__;" ";
Print __fb_backend__;" ";
#IFDEF __FB_64BIT__
Print "64Bit"
#ELSE
Print "32Bit"
#ENDIF
Print string(26,"-")

#include "clongdouble.bi"

#define INFINITY (clongdouble(1) / clongdouble(0))

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as clongdouble) as clongdouble, _
a as clongdouble, eps as clongdouble, d as clongdouble) as clongdouble

dim as clongdouble 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 clongdouble) as clongdouble, a as clongdouble, b as clongdouble, _
n as long = 6, eps as clongdouble = 1.0e-9, byref err_est as clongdouble = 0) as clongdouble

dim as clongdouble 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

' ************************************************************************

function factorial_d(n as long) as clongdouble
dim i as long
dim as clongdouble p = 1
for i  =  2 to n
p  *=  i
next
return p
end function

dim shared as long n
dim as clongdouble q, g, err_est

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

for n = 0 to 500
print "n ="; n
g = factorial_d(n)
print "n! =  "; g
err_est = 1
q = quad(@f5, 0, INFINITY, , , err_est)
print "quad ="; q; " "; csng(cdbl(err_est)); " "; csng(cdbl((g - q) / g))
print string(25, "-")
getkey
next

sleep
``````
Here are the examples from dodicat's first post:

Code: Select all

``````' example functions to integrate, exact integrals are 1, 5 and 2

function f1(x as clongdouble) as clongdouble
return acos(x)
end function

function f2(x as clongdouble) as clongdouble
return exp(- x / 5)
end function

function f3(x as clongdouble) as clongdouble
return 1 / cosh(x)^2
end function

'=================================================

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

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

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

'=================================================
#macro SimpsonArea(fn,Lo,Up,ordinates,sum)
scope
dim as clongdouble x,part1,part2
dim as clongdouble 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
'==================================================

print "integrate(acos(x), x=0..1)           = ", quad(@f1, 0, 1)
print "integrate(exp(-x/5), x=0..+inf)      = ", quad(@f2, 0, INFINITY)
print "integrate(1/cosh(x)^2, x=-inf..+inf) = ", quad(@f3, -INFINITY, INFINITY)
print "integrate(1-sqr(x*x-1)/x, x=1..+inf) = ", quad(@f4, 1, INFINITY)
print "integrate(1/(1+exp(x)), x=0..+inf)   = ", quad(@f6, 0, INFINITY)
print "integrate(x^2, x=-100.005..12.556)   = ", quad(@f5, -100.005, 120.556)

print
'Simpson macro
dim as clongdouble Result
const ordinates=1000000 ' no less anyway.
SimpsonArea(acos(x),0,1,ordinates,Result)
print "integrate(acos(x), x=0..1)           = ", Result

SimpsonArea(exp(-x/5),0,infinity,ordinates,Result)
print "integrate(exp(-x/5), x=0..+inf)      = ", Result

SimpsonArea(1/cosh(x)^2,-infinity,infinity,ordinates,Result)
print "integrate(1/cosh(x)^2, x=-inf..+inf) = ", Result

SimpsonArea(1-sqr(x*x-1)/x,1,infinity,ordinates,Result)
print "integrate(1-sqr(x*x-1)/x, x=1..+inf) = ", Result

SimpsonArea(1/(1+exp(x)),0,infinity,ordinates,Result)
print "integrate(1/(1+exp(x)), x=0..+inf)   = ", Result

SimpsonArea(x^2,-100.005,120.556,ordinates,Result)
print "integrate(x^2, x=-100.005..120.556)  = ", Result
print
print "actual value for last integral (x^2)  ";120.556^3/3- -100.005^3/3
print

sleep
``````
Arguments must be of type clongdouble:

Code: Select all

``````#include "clongdouble.bi"

dim as double a
dim as clongdouble b, c, arg = 1

a = 4*atn(1)
b = 4*atn(clongdouble(1))
c = 4*atn(arg)

print a
print b
print c
print clongdouble("3.1415926535897932384626433832795")

sleep
``````
Translated with the help from DeepL.com (free version)
dodicat
Posts: 7998
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: Numerical integration

Simpson is a local baker in these parts.
https://www.simpsonsbakery.co.uk/
So here is one of their pies, tied up.(for fun)

Code: Select all

``````
namespace PiBySimpson
const infinity=1.797693134862316e+308
Function Simpson(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer) As Double
if abs(l)=infinity then l=sgn(l)*10^4
if abs(u)=infinity then u=sgn(u)*10^4
var n=Ordinates
If n Mod 2=0 Then n=n+1
Var h=(U-L)/n
Var Part1=fn(L+.5*h)
Var Part2=0.0
For i As Integer=1 To n-1
Part1+=fn(L+h*i+h/2)
Part2+=fn(L+h*i)
Next i
Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2)
End Function
function PieSlice(x as double) as double
return exp(-x*x)
end function
dim as const double i=0
sub CreatePi() constructor
*cptr(double ptr,@i)=simpson(@PieSlice,-infinity,infinity,1000000)^2
end sub
end namespace
'======================================================
'PiBySimpson.i=9 not allowed(const)
dim as const double pi=PiBySimpson.i
const as double fbpi=acos(-1)
for n as long =1 to 5
var t=timer
print pi,
print timer-t,"pi"
t=timer
print fbpi,
print timer-t,"fbpi"
print
next
'fbpi=9 not allowed(const)
'pi=9 not allowed(const)

sleep
``````
jdebord
Posts: 551
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: Numerical integration

Here is a test using the MPFR library according to srvaldez (https://www.freebasic.net/forum/viewtopic.php?t=24110)

The library has been compiled for 100 digits precision.

The quad routine is called with 10 levels and a tolerance eps = 1E-90

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 "mpfr-fb.bi"

#define INFINITY 1/0

#define isfinite(x) (abs(x) <> INFINITY)

function exp_sinh_opt_d(func as function (x as mpfr) as mpfr, _
a as mpfr, eps as mpfr, d as mpfr) as mpfr

dim as mpfr 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 mpfr) as mpfr, a as mpfr, b as mpfr, _
n as long = 6, eps as mpfr = 1.0e-9, byref err_est as mpfr = 0) as mpfr

dim as mpfr 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)

' result with estimated relative error err_est
err_est = abs(v) / (abs(s) + eps)
return sign * d * s * h
end function

' *****************************************************************************

dim shared as long n
dim as mpfr g_ref, g_quad, err_est

function f(x as mpfr) as mpfr
return x^(n-1) * exp(-x)
end function

n = 200
g_ref = gamma(n)
? "n        : "; n
? "gamma(n) : "; g_ref
``````n        :  200