Numerical integration with tanh-sinh quadrature

From CodeCodex

Tanh-sinh quadrature is effective for calculating many integrands, even those with singularities at the endpoints. When iterative refinement is used, it is more efficient in function evaluations than standard Gaussian quadrature due to reuse of points.

Algorithm description[edit]

We approximate the integral with a sum:

 Integral [-1..1] f(x) ~= h * Sum[j=-n..n] { phi'(j*h) * f(phi(j*h)) }

where

 phi = tanh ((pi/2) * sinh(x))

and

 phi'(x) = (pi/2) * cosh(x) / ( (cosh((pi/2) * sinh(x)))^2 )

Start with h1 = 1/2 and choose a truncation value N=n*h with phi'(N) < eta, where eta is the required precision. Iteratively refine the integral by halving h and summing finer grids over this range. For each subsequent 'level' hk+1 = hk/2; we only need to sum the odd terms as the even terms are in the previous sum:

 Sum[k+1] = h_{k+1}*SumOdd[k+1] + (1/2)Sum[k]

Repeat the halving until a convergence condition is met:

 | Sum[k+1] - Sum[k] | < eta

Implementation[edit]

(* Created by Toby Kelsey 2008.03 - released to the public domain *)

let half_pi = 2.0 *. atan(1.0);;

let weight hj =
	let c = cosh (half_pi *. (sinh hj)) in
	half_pi *. (cosh hj) /. (c*.c);;

let abscissa hj = tanh (half_pi *. (sinh hj));;

(* Note phi'(x) = phi'(-x), and phi(-x) = -phi(x) *)
let term2 h j f =
	if j = 0 then
		half_pi *. (f 0.0)
	else (* terms for jh and -jh *)
	let r = h *. float_of_int j in
	let w = weight r and x = abscissa r in
	w *. ((f x) +. (f (~-.x)));;

(* Initial sum with h = 0.5. Return (n, h*sum) *)
let init_sum f eta =
	let trunc_fac = 5.0 in (* if say 5 levels then level truncation error < eta/5 *)
	let h = 0.5 in
	let rec psum n acc t = (
		let v = term2 h n f in
		let acc1 = acc +. v in
		let t1 = if eta > trunc_fac *. h *. abs_float v then (t+1) else 0 in
		if t1 > 1 then (* require 2 sequential small terms for safety *)
			n, h*.acc1
		else
			psum (n+1) acc1 t1 ) in
	psum 0 0.0 0;;

(* Sum over odd indices [n >= j >= -n] phi'(hj)*f(phi(hj)) *)
let sum_odd f h n =
	let rec psum i acc = 
		if i > n then acc else psum (i+2) (acc+.(term2 h i f)) in
	psum 1 0.0;;

(* main quadrature function *)
let ts_quad f a b e =
	let p = ((a +. b) *. 0.5) and q = ((b-.a) *. 0.5) in
	let trans x = p +. (x*.q) in (* [-1,1] => [a,b] *)
	let g x = f (trans x) in
	let n0, sum1 = init_sum g e in 
	let rec sum y h n acc = (
		let t = h *. (sum_odd y h n) in (* next level *)
		let acc1 = (acc*.0.5)+.t in (* new sum *)
		let delta = abs_float (acc -. acc1) in 
		if e > delta then (* check convergence *)
			acc1, delta (* heuristic assumes remaining error ~= delta *)
		else
			sum y (h*.0.5) (n*2) acc1 ) in (* finer grid *)
	let s, d = sum g 0.25 (2*n0) sum1 in
	(q*.s), d ;; (* 'q' factor for change of range *)


(* test section *)
open Printf;;

let ctr = ref 0;;

(* trail functions and anti-derivatives *)
let func1 x = incr ctr; sqrt x -. 1.5;;
let func2 x = incr ctr; x *. cos(x*.x);;
let f1i x = (2.0/.3.0) *. x *. sqrt x -. 1.5*.x;;
let f2i x = 0.5 *. sin(x*.x);;

(* test loop *)
let test f g a b =
let exact = g(b) -. g(a) in
for i = -1 downto -15 do
  ctr := 0;
  let eta = 10.0 ** float_of_int i in
  let v, e = ts_quad f a b eta in
  printf "target %6g nfunc %6d qval %.16g e_est %g e_acc %g\n"
    eta !ctr v e (abs_float (v-.exact));
done;;

(* run tests *)
let lo = 1.0 and hi = 6.0;;
printf "Testing func1\n";;
test func1 f1i lo hi;;
printf "Testing func2\n";;
test func2 f2i lo hi;;

Test output[edit]

Note the estimated error is the truncation error and does not include the accumulated rounding error.

Testing func1
target    0.1 nfunc     17 qval 1.631288877421028 e_est 1.99887e-05 e_acc 3.42705e-06
target   0.01 nfunc     21 qval 1.631292304226628 e_est 1.86661e-05 e_acc 2.39417e-10
target  0.001 nfunc     21 qval 1.631292304226628 e_est 1.86661e-05 e_acc 2.39417e-10
target 0.0001 nfunc     25 qval 1.631292304553765 e_est 1.8666e-05 e_acc 8.77203e-11
target  1e-05 nfunc     49 qval 1.631292304466042 e_est 3.50892e-11 e_acc 2.66454e-15
target  1e-06 nfunc     49 qval 1.631292304466042 e_est 3.50892e-11 e_acc 2.66454e-15
target  1e-07 nfunc     57 qval 1.631292304466045 e_est 3.50879e-11 e_acc 4.44089e-16
target  1e-08 nfunc     57 qval 1.631292304466045 e_est 3.50879e-11 e_acc 4.44089e-16
target  1e-09 nfunc     57 qval 1.631292304466045 e_est 3.50879e-11 e_acc 4.44089e-16
target  1e-10 nfunc     57 qval 1.631292304466045 e_est 3.50879e-11 e_acc 4.44089e-16
target  1e-11 nfunc    113 qval 1.631292304466045 e_est 0 e_acc 4.44089e-16
target  1e-12 nfunc    129 qval 1.631292304466045 e_est 0 e_acc 4.44089e-16
target  1e-13 nfunc    129 qval 1.631292304466045 e_est 0 e_acc 4.44089e-16
target  1e-14 nfunc    129 qval 1.631292304466045 e_est 0 e_acc 4.44089e-16
target  1e-15 nfunc    129 qval 1.631292304466045 e_est 0 e_acc 4.44089e-16
Testing func2
target    0.1 nfunc     65 qval -0.9166162598409749 e_est 0.00055452 e_acc 8.65928e-06
target   0.01 nfunc     81 qval -0.9166249158955178 e_est 0.000553239 e_acc 3.22999e-09
target  0.001 nfunc     81 qval -0.9166249158955178 e_est 0.000553239 e_acc 3.22999e-09
target 0.0001 nfunc    193 qval -0.9166249191254841 e_est 4.27436e-15 e_acc 2.20934e-14
target  1e-05 nfunc    193 qval -0.9166249191254841 e_est 4.27436e-15 e_acc 2.20934e-14
target  1e-06 nfunc    193 qval -0.9166249191254841 e_est 4.27436e-15 e_acc 2.20934e-14
target  1e-07 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-08 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-09 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-10 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-11 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-12 nfunc    225 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-13 nfunc    257 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-14 nfunc    257 qval -0.9166249191254983 e_est 1.66533e-15 e_acc 7.88258e-15
target  1e-15 nfunc   1025 qval -0.9166249191255038 e_est 4.44089e-16 e_acc 2.33147e-15