How to cope with the following singularity

Question

I have the following integral:

$\int_{1}^{Xd} \dfrac{(X^{z_i}-1)}{[X^2 \sum_{l=1}^{N}c_l(X^{z_l}-1)]^{1/2}}dX = \int_{1}^{Xd} h(X) dX$

where:

Xd is a real that can be either negative, positive or even 1.

$z_i$ is a integer (positive or negative)

$z_l$ also a integer ($z_i$ will be in the vector that contain all the $z_l$ of length N)

$c_l$ is a vector of positive real

You can observe that there is a singularity at the limit value 1. Where I got division of 0/0.

I am using Spyder idle (so python + libraries) and I have tried with:

import scipy.integrate as integrate

partB = integrate.quad(self.integrand_fun_Borkovec_1983_eqn_11, 1, Xd, arg = (zi,zl_vec, cl_vec))

def integrand_fun_Borkovec_1983_eqn_11 (self, x, zi,z_vec, cb):
    a = (x**zi)-1
    b= 0
    for i in range(0, len(z_vec)):
        b = b + cb[i]*((x**z_vec[i])-1)
    b = x*x*b
    return a/b

for the following values of Xd = 0.07544956110914136, zi = 1, zl_vec = [1,1,-1,-1] and cl_vec = [1.02882947e-07, 8.43609299e-05, 8.43566150e-05, 9.71978375e-08]

I get the following result in part B: (166085.34873353728, 6428.115530877825). The value 6428.115530877825 is an estimation of the error, as you can seen the estimated error is pretty high.

I have read in the paper that I am using Borkovec(1983) the following:

Expand the integrand to second order in X-1 for values of X near 1

What would mean (If I am not wrong):

$h(x) \approx h(X-1)+h'(X-1)+(1/2)h''(X-1)$

Hence, I will need to calculate h'(X-1) and h''(X-1). I wonder if there is an easier way to solve this singularity. I have tried using the points argument of the quad function, but the estimated error is still high.

partB = integrate.quad(self.integrand_fun_Borkovec_1983_eqn_11, 1, Xd, arg = (zi,zl_vec, cl_vec), points=[0.95, 1])

Thank you

Answer 1

Do not compute the expansion of $h$, compute the expansions of the simpler terms. For this purpose you can use a series approximation of \begin{align} X^z-1&=(1+(X-1))^z-1=\sum_{k=0}^\infty\binom{z}{k}(X-1)^k\\ &=z(X-1)+\frac{z(z-1)}2(X-1)^2+\frac{z(z-1)(z-2)}6(X-1)^3+... \end{align} where the terms in the expansion get rapidly small if $X-1$ is small enough.

def term(X,z):
    if abs(X-1)>1e-6: return X**z-1;
    return z*(X-1)*(1+(z-1)*(X-1)/2.0*(1+(z-2)*(X-1)/3.0*(1+(z-3)*(X-1)/4.0)));

The integral under consideration is $$ \int_1^{X_d}\frac{X^{z_i}-1}{\sqrt{X^2\sum_{l=1}^N c^B_l(X^{z_l}-1)}} $$ under the condition $\sum_{l=1}^N c^B_l z_l=0$. This has the consequence that close to $X=1$ the integrand is approximately $$ \frac{ z_i+\frac{z_i(z_i-1)}2(X-1)+... }{ X\sqrt{\left[\sum_{l=1}^N c^B_l \frac{z_l(z_l-1)}2\right] +\left[\sum_{l=1}^N c^B_l \frac{z_l(z_l-1)(z_l-2)}6\right](X-1)+...} } $$ which does no longer have a singularity at $X=1$. To compensate for the rounding error in the given coefficients, compute the defect $d=\sum_{l=1}^N c^B_l z_l$ and remove $d\,(X-1)$ from the sum in the denominator. This then clears all negative values that previously occurred and thus removes the necessity to "sanitize" the value of b.

def integrand_fun_Borkovec_1983_eqn_11 (x, zi,z_vec, cb):
    a = term(x,zi)
    b = sum(cB*term(x,z) for z,cB in zip(z_vec,cb))
    d = sum(cB*z for z,cB in zip(z_vec,cb))
    b = x*x*(b-d*(x-1))
    if b<0: print x,b
    return a/b**0.5

Now running the integration with the routines of scipy gives a value with a low error bound

import scipy.integrate as integrate
Xd = 0.07544956110914136
zi = 1
zl_vec = [1,1,-1,-1] 
cl_vec = [1.02882947e-07, 8.43609299e-05, 8.43566150e-05, 9.71978375e-08]
print "sum z[i]*cB[i] should be zero, but is",sum(cB*z for z,cB in zip(zl_vec,cl_vec))

partB = integrate.quad(integrand_fun_Borkovec_1983_eqn_11, 1, Xd, args = (zi,zl_vec, cl_vec))
print "partB =", partB

This gives as result

sum z[i]*cB[i] should be zero, but is 1.00000095e-08
partB = (157.85179247575064, 2.291846657454834e-06)

which is good enough relative to the accuracy of the provided data.