how to calculate this integral using scipy/numpy?

i have the folowing question enter image description here it's written badly, but for those values of alpha and beta (0.32, 12) i need to find the minimal k at which the integral from -inf to inf of that expression is less than pi.

but this code I tried didn't work as for the example I got the value of 4 and not 15 so I don't know how to continue.

this is what i got with some help of chatGPT, when the last print shows a positive value that means that i've met my condition and it happens at k=2.

import numpy as np
from scipy.integrate import quad

def a(i):
    return (12 * i - 11) ** 0.32

def sinc(t):
    if t == 0:
        return 1.0
    else:
        return np.sin(t) / t

def p(x, k):
    product = 1.0
    for i in range(1, k + 1):
        product *= sinc(x / a(i))
    return product

def P(k):
    integral, _ = quad(lambda x: p(x, k), -np.inf, np.inf)
    return integral

for value in range(1,5):
    print(f"P({value}) =", P(value))
    print(f'pi - P({value}) =', np.pi - P(value))

### solution was 2

Interestingly, most methods provided by scipy.integrate yields a number very close to pi but less than pi from k = 2 onwards. eg quad, fixed_quad, quadrature, romberg. Even using the diffferent gaus quadratures in scipy.special.root_* yields to a value less than pi by k=2.

import numpy as np
import scipy
    
def fn(x, k, alpha, beta):
  x = np.array(x).reshape(-1,1)
  w =  x / (beta * np.arange(1, k + 1) - (beta - 1)) ** alpha / np.pi
  return np.sinc(w).prod(1)

scipy.integrate.quad(fn, -np.inf, np.inf, (5, 1, 2))
(3.1415926520655124, 4.664127099015817e-08)

This value is already less than pi and k is only 5 so had to try other methods

In order to use the other methods, a vectorized function that is scaled to [-1,1] interval has to be defined:

def fn2(t, k, alpha, beta):
  t =  np.array(t).reshape(-1)
  results = np.zeros_like(t)
  index = np.abs(t) != 1
  t_ = t[index]
  denom = 1 - t_**2
  x = (t_/denom).reshape(-1,1)
  w =  x / (beta * np.arange(1, k + 1) - (beta - 1)) ** alpha / np.pi
  results[index] =  np.sinc(w).prod(1)*(1 + t_ **2) / denom**2
  return results


scipy.integrate.quadrature(fn2, -1,1,(10,1,2), maxiter=300, vec_func = True)
(3.1415900583875214, 3.7533582020188305e-08)

scipy.integrate.romberg(fn2, -1,1,(7,1,2), divmax = 20, vec_func = True)
array([3.14159266])

Still no success in the integration since the value is less than pi for K=10

Lets try julia:

## %%julia
using Pkg
using QuadGK

function fn(x, k, alpha, beta)
  w = x ./(beta .* range(1, k) .- (beta - 1)).^alpha
  x == 0 ? 1. : prod(sin.(w)./w)
end
quadgk(x -> fn(x, 5, 1., 2.), -Inf, Inf)
(3.1415926432027703, 4.3290522541976784e-8)

Also no success. since value is less than pi for k = 5

Mathematica?

beta = 2;
alpha = 1;
For[k = 2, k<10, k++, 
    Print["k=", k, " ", Integrate[Product[Sinc[x/(beta*i - beta+1)^alpha],{i,1,k}], {x, -Infinity, Infinity}]<Pi]]

Results:

k=1 False
k=2 False
k=3 False
k=4 False
k=5 False
k=6 False
k=7 False
k=8 True
k=9 True

So, Mathematica says that at k=8 the result is already less than pi.

What about R?

## %%R
fn <- function(x, k, alpha, beta){
  a <- (beta * seq_len(k) - beta + 1)**alpha
  ifelse(x==0, 1, apply(sin(w <- outer(x, a, '/'))/w, 1, prod))
}
integrate(fn, -Inf, Inf, k=10, alpha=1,beta=2)
3.141601 with absolute error < 0.00012

This seems to be successful. Compute for K such that the integral is less than pi.

##%%R
g <- function(k) {
  integrate(fn, -Inf, Inf, k=k, alpha=1,beta=2)$value - pi
}
uniroot(g, c(3,20))
$root
[1] 14.99996

$f.root
[1] 1.430762e-07

$iter
[1] 19

$init.it
[1] NA

$estim.prec
[1] 6.103516e-05

R gives the root as 15. You can use minimize/optim and both methods returns 15.

Since R uses Gauss-kronrod quadrature nodes and values, there should be a python implementation for that. but could not find any.

next try fixed_quad playing around with n. n = 200 seems to match the results from R:

scipy.integrate.fixed_quad(fn2, -1,1,(10,1,2), 200)
(3.141602220677316, None)

g = lambda k : scipy.integrate.fixed_quad(fn2, -1,1,(int(k),1,2), n = 200)[0] - np.pi
scipy.optimize.root_scalar(g, bracket=(4, 20))
      converged: True
           flag: converged
 function_calls: 53
     iterations: 52
           root: 15.000000000000028

Notice that the root = 15. Wow but that is trial and error with regards to the n, without which one cannot obtain the needed solution.

Other approximations?

nodes1, weights1 = scipy.special.roots_jacobi(1000,0,0)
fn2(nodes1,10,1,2)@weights1
3.1415923602907556
nodes2, weights2 = scipy.special.roots_legendre(1000)
fn2(nodes2,10,1,2)@weights2
3.1415923602907556

These values are still less than pi. So why Only R gives the results yet every other method fails?

Perhaps the author of that question used R. Every other approximation is close to pi, but mostly less than pi.

NB: I could not find the gauss-kronrod in python.

Conclusion: This is an ill-defined problem as the results of the integration need to be pi. thus an approximation could be smaller or larger than pi, and due to numerical approximations, different platforms would yield different results. eg the fixed_quad method used above will yield different results for different n. its not common to have a specific n that would yield the needed results. I had to optimize to find n=200.

For the case beta=2, alpha=1 the correct answer is k=8. (Kudos to Jared for his comment in another thread).

The confusion over the number 15 is that 2k-1 = 15 in that instance.

This was proved by Borwein and Borwein, "Some remarkable Properties of Sinc and Related Integrals", The Ramanujan Journal, 5, 73-89, 2001. You can find it at https://www.carmamaths.org/resources/db90/pdfs/db90-119.00.pdf Look at their equation 10 and the text around it. (Note that their integral runs from 0 to infinity, so is half the current one. To reiterate, it is k=8 and 2k-1=15 that is the interesting case).

For every k value below this Nate3384's integral (which runs from -infinity to infinity rather than 0 to infinity) would be exactly pi.

I include a screenshot from their paper. The paragraph after equation 10 is amusing (and salutary). I don't think any purely numerical method would have the accuracy to get this - you would need symbolic maths, such as Mathematica or, possibly, Sympy. So: over to the Sympy experts!

Screenshot from Borwein and Borwein's paper:

import numpy as np
from scipy.integrate import quad
from numpy import sinc

# Define the parameters
beta = 12
alpha = 32 / 100

# Define the series a_i
def a_i(i):
    return ((beta * i - beta + 1) ** alpha)

# Define the integrand
def integrand(x, k):
    product = 1
    for i in range(1, k + 1):
        product *= sinc(x / a_i(i))
    return product

# Function to find the minimum k such that the integral is less than π
def find_min_k():
    k = 1
    while True:
        integral_value, _ = quad(integrand, -np.inf, np.inf, args=(k,))
        if integral_value < np.pi:
            return k
        k += 1

# Find the result
result = find_min_k()
print(f"The minimum k such that the integral is less than π is: {result}")