Borwein integrals

This notebook (which covers roughly the same ground as this video) is about the Borwein integrals. You can find more information on Wikipedia or in the original paper (titled Some remarkable properties of sinc and related integrals) by Borwein and Borwein. (To read this or other research papers in non-free journals, you should log in to MUSE, select StarPlus under My Services, then enter the title in the search bar.)

Consider the functions $$ \begin{array}{rl} f_0(x) &= \sin(x)/x \\ f_1(x) &= 3\sin(x)\sin(x/3)/x^2 \\ f_2(x) &= 3\cdot 5\sin(x)\sin(x/3)\sin(x/5)/x^3 \\ f_3(x) &= 3\cdot 5\cdot 7\sin(x)\sin(x/3)\sin(x/5)\sin(x/7)/x^4 \\ \dotsb & \dotsb \\ f_7(x) &= 3\cdot 5\cdot 7\cdot 9\cdot 11\cdot 13\cdot 15 \sin(x)\sin(x/3)\sin(x/5)\sin(x/7)\sin(x/9)\sin(x/11)\sin(x/13)\sin(x/15)/x^8, \end{array} $$ and the integrals $$ I_k = \int_{x=0}^\infty f_k(x)\,dx. $$ It is a remarkable fact that $$ I_0=I_1=I_2=I_3=I_4=I_5=I_6=\pi/2, $$ but $$ I_7 = \frac{467807924713440738696537864469}{935615849440640907310521750000} \pi. $$ This is almost exactly the same as $\pi/2$, but not quite: we have $I_7\approx\pi/2-2.31\times 10^{-11}$. (There is a brief explanation of this phenomenon in Wikipedia, and details in the original paper.) In order to check this by numerical integration, we need to be able to evaluate integrals to very high precision. Some very sophisticated methods have been developed for this, and are implemented in libraries such as scipy (specifically, in the function scipy.integrate.quad() and related functions). Our aim here is

• To investigate some simpler integration methods, as an introduction to the ideas needed when developing more sophisticated methods.
• To learn how to read and understand the scipy documentation well enough to supply appropriate options to the integration routines and so achieve the required accuracy.

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

First, for future reference we record the claimed value of $I_7$ and check the difference from $\pi/2$.

I7 = 467807924713440738696537864469/935615849440640907310521750000 * np.pi
print(f"The difference between I7 and pi/2 is {np.pi/2 - I7}")

The difference between I7 and pi/2 is 2.3100410473375632e-11

We now define the functions $f_0(x),\dotsc,f_7(x)$. These involve constant multipliers which we call $m_0,\dotsc,m_7$.

m0 = 1
m1 = 3
m2 = 3*5
m3 = 3*5*7
m4 = 3*5*7*9
m5 = 3*5*7*9*11
m6 = 3*5*7*9*11*13
m7 = 3*5*7*9*11*13*15

We now define the functions $f_k(x)$ themselves. We could use the usual syntax like

def f1(x):
    return m1 * np.sin(x) * np.sin(x/3) * x ** (-2)

but we choose to use the keyword lambda instead. This is an alternative way to define one-line functions that is convenient in a number of contexts; here we just use it to make the code more compact.

f0 = lambda x: m0 * np.sin(x) / x 
f1 = lambda x: m1 * np.sin(x) * np.sin(x/3) / x ** 2
f2 = lambda x: m2 * np.sin(x) * np.sin(x/3) * np.sin(x/5) / x ** 3
f3 = lambda x: m3 * np.sin(x) * np.sin(x/3) * np.sin(x/5) * np.sin(x/7) / x ** 4
f4 = lambda x: m4 * np.sin(x) * np.sin(x/3) * np.sin(x/5) * np.sin(x/7) * np.sin(x/9) / x ** 5
f5 = lambda x: m5 * np.sin(x) * np.sin(x/3) * np.sin(x/5) * np.sin(x/7) * np.sin(x/9) * np.sin(x/11) /  x ** 6
f6 = lambda x: m6 * np.sin(x) * np.sin(x/3) * np.sin(x/5) * np.sin(x/7) * np.sin(x/9) * np.sin(x/11) * np.sin(x/13) / x ** 7
f7 = lambda x: m7 * np.sin(x) * np.sin(x/3) * np.sin(x/5) * np.sin(x/7) * np.sin(x/9) * np.sin(x/11) * np.sin(x/13) *  np.sin(x/15) / x ** 8

It would be more satisfactory to give a single definition of $f_n(x)$ that works for all $n$. We can do this as follows:

def fn_basic(n, x):
    m = np.arange(1, 2*n+3, 2) # m = [1, 3, 5, ... , 2*n+1]
    # Note that x / m is the array [x/1, x/3, x/5, ... , x/(2*n+1)]
    # so np.sin(x / m) is the array [sin(x/1), sin(x/3), sin(x/5), ... , sin(x/(2*n+1))]
    # Also the .prod() method gives the product of all the elements in the array
    return m.prod() * np.sin(x / m).prod() / x ** (n+1)

There are two issues with the above definition:

• It only works when x is a scalar, not when x is an array. If x is an array then the expression x / m will typically give an error, because x and m have incompatible shapes and cannot be divided. If x happens to have the same shape as m then the expression x / m will not give an error, but the function will not return the correct result. As an exercise, you can analyse exactly what happens when n=1 and x=np.array([u,v]) say.
• With the above definition we have to enter fn_basic(7, x) to get $f_7(x)$. It would be better if we could enter f7=fn(7) and then f7(x).

These issues are fixed by the alternative version below. To explain how it works:

• The function np.expand_dims(x,-1) reshapes x by replacing each numerical entry a by [a]. For example, if x=[[1,2,3],[4,5,6]] (of shape (2,3)) then np.expand_dims(x,-1)=[[[1],[2],[3]],[[4],[5],[6]]] (of shape (2,3,1)).
• We then divide by the array m = np.arange(1, 2*n+3, 2)=[1,3,5,...,2*n+1], which has shape (2*n+1,). The division involves broadcasting, and the result is to replace each numerical entry a in the original array by [a, a/3, a/5, ..., a/(2*n+1)], giving a new array u. In the above example where x has shape (2,3), the array u will have shape (2,3,2*n+1).
• Now np.sin(u) / u is the array obtained from x by replacing each numerical entry a by the vector [np.sin(a)/a,...,np.sin(a/(2*n+1))/(a/(2*n+1))].
• We now take (np.sin(u) / u).prod(axis=-1) to take the product of each of these innermost vectors, giving an array of the same shape as x, in which each numerical value a has been replaced by $f_n(a)$.
• Inside the function fn, we define a function f and we return it. This means that we can enter fn(9)(x) to get $f_9(x)$, or equivalently we can enter f9=fn(9) and then f9(x). Although we will not explore the issues properly here, we will mention that this definition works correctly because of technical details about the scope of the variable n, which you might need to understand when doing more complicated things along similar lines.

def fn(n):
    def f(x):
        u = np.expand_dims(x,-1) / np.arange(1, 2*n+3, 2)
        return (np.sin(u) / u).prod(axis=-1)
    return f

When x is a scalar, the functions f7(x) and fn(7)(x) and fn_basic(7, x) all work correctly and give the same answer.

x = 42
print(f7(x), fn_basic(7, x), fn(7)(x))

-8.469186619278489e-10 -8.469186619278489e-10 -8.469186619278489e-10

When x is an array, the functions f7(x) and fn(7)(x) work correctly and give the same answer, and we can also get the same answer less efficiently by using fn_basic in a nested list comprehension, but fn_basic(7, x) gives an error.

x = np.array([[1,2,3],[4,5,6]])
print(f7(x))
print(fn(7)(x))
print(np.array([[fn_basic(7, a) for a in r] for r in x])) # same as above, but not vectorized
# print(fn_basic(7, x)) # gives an error, as explained above

[[ 0.81347988  0.39670628  0.03447635]
 [-0.10772798 -0.07761746 -0.01201976]]
[[ 0.81347988  0.39670628  0.03447635]
 [-0.10772798 -0.07761746 -0.01201976]]
[[ 0.81347988  0.39670628  0.03447635]
 [-0.10772798 -0.07761746 -0.01201976]]

We can gain some insight by plotting the functions $f_6(x)$ and $f_7(x)$.

xs = np.linspace(0.01, 20, 1000)
fig, ax = plt.subplots(1,2)
fig.set_figwidth(15)
ax[0].plot(xs, fn(6)(xs), 'r-', label=r"$f_6(x)$")
ax[1].plot(xs, fn(7)(xs), 'b-', label=r"$f_7(x)$")
ax[0].spines['bottom'].set_position('zero') # put the bottom spine at y=0
ax[1].spines['bottom'].set_position('zero') # put the bottom spine at y=0
ax[0].legend()
ax[1].legend()

<matplotlib.legend.Legend at 0x26f9ce18dd0>

You can see that the two graphs look identical. Why is this? From the definitions we see that $$ f_7(x)=15f_6(x)\sin(x/15)/x. $$
Equivalently, if we write $g(t)=\sin(t)/t$ and $t=x/15$ then $f_7(x)=f_6(x)g(t)$. It is standard that when $t\approx 0$ we have $\sin(t)\approx t$ and so $g(t)\approx 1$. In the left hand half of the above graphs we have $0\leq x\leq 10$ so $0\leq t\leq 0.67$ so $0.92\leq g(t)\leq 1$ so $f_7(x)$ is close to $f_6(x)$. In the right hand half of the graph the expression $g(t)$ becomes smaller but both functions are already very small so you cannot see the difference between them. It is therefore not a surprise that $I_6$ and $I_7$ are close.

We next define a function I(f) that calls a scipy routine to calculate $\int_{x=0}^\infty f(x)\,dx$ numerically. Numerical evaluation of integrals is often called quadrature, so the appropriate scipy routine is called sp.integrate.quad(). We could enter help(sp.integrate.quad) to see the documentation for this, but it is easier to read it on the web at docs.scipy.org.

• The first argument f is the function to integrate
• The second and third arguments are the lower and upper limits of the region of integration. We want to integrate from $x=0$ to $x=\infty$, so we use the arguments 0 and np.inf. You should not generally assume that it is OK to use np.inf in Python functions, but in this case the documentation linked above tells us explicitly that np.inf is allowed.
• The integration routine is adaptive: it calculates f(a_i) at various points a_i, and tries to work out the integral based on those values. If it appears that f is changing in a relevant way in between the values that have already been calculated, then it will throw in some additional evaluation points and recompute. By default, it will only do this kind of refinement 50 times. We want an extremely accurate result, and our functions $f_n(x)$ have features that make integration difficult. We therefore add the argument limit=5000 to specify that the integral should be refined as many as 5000 times.
• We also add the argument epsabs=1e-14 to ask for an answer that is accurate with an error of less than $10^{-14}$. Of course the integration routine does not know the true answer and so does not have definitive knowledge of the accuracy of its approximation. Nonetheless, the standard algorithms include a method to estimate the probable size of the error, and that is what is relevant here.
• Similarly, we include the argument epsrel=1e-14 to say that the relative error should be less than $10^{-14}$, i.e. that the error should be less than $10^{-14}$ times the true answer. As the true answer is $\pi/2\approx 1.57$, this should not be very different from the effect of epsabs=1e-14. However, it seems that we get much better results by specifying both arguments. It would be an interesting exercise to understand in detail why this is, but probably quite a difficult exercise.

def I(f):
    """Compute the integral of f from 0 to infinity.

    Parameters:
    f (function): the function to integrate.

    Returns:
    A pair of numbers: the value of the integral and an estimate of the error.
    """
    return sp.integrate.quad(f, 0, np.inf, limit=5000, epsabs=1e-14, epsrel=1e-14)

def check_I(f, u):
    """Check the integral of f against the expected value u."""
    integral, err = I(f)
    print(f"Approximate integral           = {integral}")
    print(f"Error estimated by scipy       = {err}")
    print(f"Difference from expected value = {integral - u}\n")

check_I(f4, np.pi/2)
check_I(f5, np.pi/2)
check_I(f6, np.pi/2)
check_I(f7, I7)

Approximate integral           = 1.5707963267948974
Error estimated by scipy       = 2.3350579851661957e-13
Difference from expected value = 8.881784197001252e-16

Approximate integral           = 1.5707963267948963
Error estimated by scipy       = 8.711381198751283e-14
Difference from expected value = -2.220446049250313e-16

Approximate integral           = 1.5707963267948946
Error estimated by scipy       = 4.9057396950027255e-14
Difference from expected value = -1.9984014443252818e-15

Approximate integral           = 1.5707963267717977
Error estimated by scipy       = 3.6818168092483253e-14
Difference from expected value = 1.5543122344752192e-15

C:\Users\npstr\AppData\Local\Temp\ipykernel_18384\1988469615.py:10: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return sp.integrate.quad(f, 0, np.inf, limit=5000, epsabs=1e-14, epsrel=1e-14)

We now define two standard methods for numerical integration: the trapezium rule and Simpson's rule. We want to approximate $\int_{x=a}^bf(x)\,dx$ for some given function $f(x)$ and some interval $[a,b]$. We are also given an integer $n>0$, and we start by splitting $[a,b]$ into $n$ intervals $I_i=[u_i,u_{i+1}]$ of equal length, where $u_i=a+(i/n)(b-a)$ for $0\leq i\leq n$. We then put $y_i=f(u_i)$. If $n$ is large then the part of the graph between $(u_i,y_i)$ and $(u_{i+1},y_{i+1})$ will be close to a straight line segment. The area under that straight line segment is $(y_i+y_{i+1})(b-a)/(2n)$, and we can add up these contributions to get an approximation to the integral. In the case $n=5$, then answer is $$ \int_a^bf(x)\,dx\approx \frac{b-a}{2\times 5}\left((y_0+y_1)+(y_1+y_2)+(y_2+y_3)+(y_3+y_4)+(y_4+y_5)\right) = \frac{b-a}{5}\left(\frac{1}{2}y_0+y_1+y_2+y_3+y_4+\frac{1}{2}y_5\right). $$ The general case follows the same pattern. This approximation is called the trapezium rule. On the other hand, Simpson's rule is given by a slightly different formula, which works as follows when $n=6$: $$ \int_a^bf(x)\,dx\approx \frac{b-a}{3\times 6}\left(y_0+4y_1+2y_2+4y_3+2y_4+4y_5+y_6\right). $$ In the general case we have a denominator of $3n$ outside the brackets, and $y_0$ and $y_n$ appear with coefficient one, and the other values $y_i$ occur with coefficients alternating between $4$ and $2$. This method is only applicable when $n$ is even. We will not explain where this rule comes from, but we will mention the key property: if $f(x)$ is a polynomial of degree three or less, then Simpson's rule will give the exact answer for the integral (even with $n=2$), but the trapezium rule will only give the exact answer if the degree is zero or one.

def trapezium_rule(f, a, b, n):
    """Compute the integral of f from a to b using the trapezium rule.

    Parameters:
    f (function): the function to integrate.
    a (float): the lower limit of integration.
    b (float): the upper limit of integration.
    n (int): the number of subintervals to use.

    Returns:
    The value of the integral.
    """
    x = np.linspace(a, b, n+1) # n+1 points make n subintervals
    y = f(x)
    return (b-a) * (np.sum(y) - 0.5 * (y[0] + y[-1])) / n

def simpsons_rule(f, a, b, n):
    """Compute the integral of f from a to b using Simpson's rule.

    Parameters:
    f (function): the function to integrate.
    a (float): the lower limit of integration.
    b (float): the upper limit of integration.
    n (int): the number of subintervals to use.

    Returns:
    The value of the integral.
    """
    if not (isinstance(n, int) and n > 0 and n % 2 == 0):
        raise ValueError("n must be an even positive integer")
    x = np.linspace(a, b, n+1) # n+1 points make n subintervals
    y = f(x)
    h = (b-a) / n
    return h/3 * (y[0] + y[-1] + 2*np.sum(y[2:-1:2]) + 4*np.sum(y[1:-1:2]))

def scipy_rule(f, a, b, n):
   return sp.integrate.quad(f, a, b, limit=5000, epsabs=1e-14, epsrel=1e-14)[0]

def check_rule(rule, f, a, b, n):
    """Check the integral of f using the given rule."""
    int = rule(f, a, b, n)
    true_int = scipy_rule(f, a, b, n)
    print(f"Approximate integral = {int}")
    print(f"Accurate integral    = {true_int}")
    print(f"Difference           = {int - true_int}")
    print('')
    

Above we claimed that Simpson's rule gives the exact answer for cubic polynomials. We now check a case of this.

def p(x):
    return 4 + 6*x + 17*x**2 - 3*x**3

check_rule(simpsons_rule, p, 4, 7, 2)

Approximate integral = 83.25
Accurate integral    = 83.24999999999994
Difference           = 5.684341886080802e-14

C:\Users\npstr\AppData\Local\Temp\ipykernel_18384\2679450926.py:37: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return sp.integrate.quad(f, a, b, limit=5000, epsabs=1e-14, epsrel=1e-14)[0]

We next introduce a much less obvious integration rule called five-point Gaussian quadrature. This takes its simplest form when we are integrating from $-1$ to $1$, where we get $$ \int_{-1}^1 f(x)\,dx \approx w_2f(-x_2) + w_1f(-x_1) + w_0f(x_0) + w_1f(x_1) + w_2f(x_2), $$ where $$ \begin{array}{rlcrl} x_0 &= 0 && w_0 &= 128/225 \approx 0.5689\\ x_1 &= \frac{1}{3}\sqrt{5-2\sqrt{10/7}}\approx 0.5385 && w_1 &= (322 + 13\sqrt{70})/900 \approx 0.4786\\ x_2 &= \frac{1}{3}\sqrt{5+2\sqrt{10/7}}\approx 0.9062 && w_2 &= (322 - 13\sqrt{70})/900 \approx 0.2369. \end{array} $$ These numbers are carefully chosen to give the following key property: the rule gives the exact answer for all polynomials of degree less than or equal to six. The formula above is only valid for integration from $-1$ to $1$, but we can use a simple change of variables to reduce the general case to that one. This rule works very well for integration over a short interval where the graph does not have too many features. For longer intervals, we can split the interval into n equal parts for some given number n>0, use the five-point quadrature rule on each part, and then add up the results. This is implemented by the function split_gauss5_rule() defined below.

Wikipedia will tell you many more things about Gaussian quadrature if you are interested.

gauss5_points = np.array([
    -np.sqrt(5 + 2*np.sqrt(10/7))/3, 
    -np.sqrt(5 - 2*np.sqrt(10/7))/3,
     0,
     np.sqrt(5 - 2*np.sqrt(10/7))/3,
     np.sqrt(5 + 2*np.sqrt(10/7))/3,
])

gauss5_weights = np.array([
    (322 - 13*np.sqrt(70))/900,  
    (322 + 13*np.sqrt(70))/900,
    128/225,
    (322 + 13*np.sqrt(70))/900,
    (322 - 13*np.sqrt(70))/900  
])

def gauss5_rule(f, a, b, n = None):
    """Compute the integral of f from a to b using a 5-point Gauss rule.

    Parameters:
    f (function): the function to integrate.
    a (float): the lower limit of integration.
    b (float): the upper limit of integration.
    n (int): ignored, included for compatibility with other functions.

    Returns:
    The value of the integral.
    """
    x = 0.5 * (b-a) * gauss5_points + 0.5 * (b+a)
    return 0.5 * (b-a) * np.sum(gauss5_weights * f(x))

def split_gauss5_rule(f, a, b, n):
    return np.array([gauss5_rule(f, a+i*(b-a)/n, a+(i+1)*(b-a)/n) for i in range(n)]).sum()

We now check our various rules on a highly oscillatory polynomial of degree six. Even with 100 points, the trapezium rule gives poor accuracy and Simpson's rule is mediocre. The five-point Gaussian rule gives the exact answer with no subdivision.

def p5(x):
    return x * (x - 1) * (x - 2) * (x - 3) * (x - 4) * (x - 5)

check_rule(trapezium_rule, p5, 1, 6, 100)
check_rule(simpsons_rule, p5, 1, 6, 100)
check_rule(gauss5_rule, p5, 1, 6, None)

Approximate integral = 221.49341332031273
Accurate integral    = 221.13095238095246
Difference           = 0.36246093936026114

Approximate integral = 221.13110859375016
Accurate integral    = 221.13095238095246
Difference           = 0.00015621279769106877

Approximate integral = 221.13095238095246
Accurate integral    = 221.13095238095246
Difference           = 0.0

C:\Users\npstr\AppData\Local\Temp\ipykernel_18384\3177476011.py:37: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return sp.integrate.quad(f, a, b, limit=5000, epsabs=1e-14, epsrel=1e-14)

We now try to calculate $\int_1^{10}\sin(x)\,dx$. Here the trapezium rule (with $n=100$) is OK, and Simpson's rule (with $n=100$) is better. The five-point Gaussian rule is quite poor, but if we split the interval into five parts first then we get a very accurate answer. Note that this only involves evaluating the function at $25$ points, whereas the other two methods used $101$ points.

check_rule(trapezium_rule, np.sin, 1, 10, 100)
check_rule(simpsons_rule, np.sin, 1, 10, 100)
check_rule(gauss5_rule, np.sin, 1, 10, None)
check_rule(split_gauss5_rule, np.sin, 1, 10, 5)

Approximate integral = 1.378442631886318
Accurate integral    = 1.3793738349445925
Difference           = -0.0009312030582744324

Approximate integral = 1.3793743382115935
Accurate integral    = 1.3793738349445925
Difference           = 5.032670009619267e-07

Approximate integral = 1.3736652384967674
Accurate integral    = 1.3793738349445925
Difference           = -0.005708596447825176

Approximate integral = 1.3793738351641145
Accurate integral    = 1.3793738349445925
Difference           = 2.1952195616847803e-10

C:\Users\npstr\AppData\Local\Temp\ipykernel_18384\3177476011.py:37: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return sp.integrate.quad(f, a, b, limit=5000, epsabs=1e-14, epsrel=1e-14)

Finally, we try to evaluate the Borwein integral $\int_0^\infty f_7(x)\,dx$ using the rules defined above. We cannot take $a=0$ because the definition of $f_7(x)$ involves division by $x$. (It is easy to see that $f_7(x)\to 1$ as $x\to 0$, but that is not visible to Python.) We therefore take $a$ to be a very small positive number instead, specifically $a=10^{-14}$. We also cannot take $b=\infty$, because all of our methods assume that $b$ is finite. We therefore take $b=1000$ in the expectation that the function is so small for $x\geq 1000$ that this cutoff will make a negligible difference.

print(f"trapezium rule  : {trapezium_rule(   fn(7), 1e-14, 1000, 1000) - I7}")
print(f"Simpson's rule  : {simpsons_rule(    fn(7), 1e-14, 1000, 1000) - I7}")
print(f"split Gauss rule: {split_gauss5_rule(fn(7), 1e-14, 1000, 1000) - I7}")
print(f"scipy rule      : {scipy_rule(       fn(7), 1e-14, 1000, None) - I7}")

trapezium rule  : -9.992007221626409e-15
Simpson's rule  : -9.992007221626409e-15
split Gauss rule: -1.021405182655144e-14
scipy rule      : -1.0658141036401503e-14

C:\Users\npstr\AppData\Local\Temp\ipykernel_18384\2679450926.py:37: IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  return sp.integrate.quad(f, a, b, limit=5000, epsabs=1e-14, epsrel=1e-14)[0]