7

I'm trying to write a Python program to use Tanh-sinh quadrature to compute the value of \begin{equation} \int_{-1}^1 \frac{dx}{\sqrt{1-x^2}} \end{equation} but although the program converges to a sensible value with no errors in every case, it's not converging to the correct value (which is $\pi$ for this particular integral) and I can't find the problem.

Instead of asking for a desired level of accuracy, the program asks for the number of function evaluations wanted, to make comparisons of convergence with simpler integration methods easier. The number of evaluations needs to be an odd number as the approximation used is \begin{equation} \int_{-1}^1 f(x) dx = \sum_{k=-n}^n w_k f(x_k) \end{equation} Can anyone suggest what I might have done wrong?

import math

def func(x):
    # Function to be integrated, with singular points set = 0
    if x == 1 or x == -1 :
        return 0
    else:
        return 1 / math.sqrt(1 - x ** 2)

# Input number of evaluations
N = input("Please enter number of evaluations \n")
if N % 2 == 0:
    print "The number of evaluations must be odd"
else:
    print "N =", N  

# Set step size
h = 2.0 / (N - 1)
print "h =", h

# k ranges from -(N-1)/2 to +(N-1)/2
k = -1 * ((N - 1) / 2.0)
k_max  = ((N - 1) / 2.0)
sum = 0

# Loop across integration interval
while k < k_max + 1:

    # Compute abscissa
    x_k = math.tanh(math.pi * 0.5 * math.sinh(k * h))

    # Compute weight
    numerator = 0.5 * h * math.pi * math.cosh(k * h)
    denominator = math.pow(math.cosh(0.5 * math.pi * math.sinh(k * h)),2)
    w_k =  numerator / denominator

    sum += w_k * func(x_k)

    k += 1

print "Integral =", sum
ThomasNicholas
  • 171
  • 4

3 Answers3

15

There are a bunch of pitfalls when it comes to tanh-sinh quadrature, one being that the integrand needs to be evaluated very closely to the interval boundaries, at distances less than machine precision, e.g., 1.0 - 1.0e-20 in the original example. When this point is evaluated, it rounds to 1.0 at which f has a singularity, and anything can happen. You're getting the wrong result because you're setting that value to 0.

You'll first have to transform the function such that the singularities sit at 0 to avoid the round-off. In the case of 1 / sqrt(1 - x**2), this is 1 / numpy.sqrt(2*x - x**2) for both the left and the right singularity. With tanh_sinh (a project of mine), one then gets

import numpy as np
import tanh_sinh

def f(x):


return 1 / np.sqrt(1 - x ** 2)


val, error_estimate = tanh_sinh.integrate_lr(
    lambda x: 1 / np.sqrt(2 * x - x ** 2),  # = 1 / sqrt(1 - (x-1)2)
    lambda x: 1 / np.sqrt(2 * x - x  2),  # = 1 / sqrt(1 - (-(x-1))**2)
    2.0,  # length of the interval
    1.0e-10,
)
print(val, val - np.pi)

3.1415926533203944 -2.693987255497632e-10
Nico Schlömer
  • 3,126
  • 17
  • 36
5

Your approach to selecting $h$ based on $N$ isn't a very sophisticated one and probably isn't getting $h$ to be small enough. See these notes for a more reasonable way to adjust $h$ and estimate the error in the integral evaluation:

http://www.davidhbailey.com/dhbpapers/dhb-tanh-sinh.pdf

Brian Borchers
  • 18,719
  • 1
  • 39
  • 70
0

One more remark:

I would suggest you to read the original papers on the development and implementation of the DE quadrature (see here for some papers) and try to understand the implementation by Prof. Ooura that you can download here (make sure to download what he calls "the easy version", since that one is more or less readable). The straight-forward calculation of the weights $w_k$ as you do also leads to round-off errors.

GertVdE
  • 6,179
  • 1
  • 21
  • 36