# How to cope with the following singularity

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

• What is $[c_l(X^{z_l}-1)]^{1/2}$ when the radicant is negative? This you get for $l=1,2$, and $X=0.1$ for example, which is in the interval from $X_d$ to $1$ that you integrate over. I notice that this square root does not appear in your code. Dec 14 '18 at 21:57
• Could you add an image of the formula in the original paper for reference? Dec 14 '18 at 22:21
• cl is a concentration, hence a positive real, and zl is the charge (e.g. for H+ it is 1, for Cl- it is -1, for Mg+2 it is 2, ...). Right now, I am not that worry about a negative radicand, that it is also a problem, but more about the fact that at X = 1, there is a singularity. Dec 17 '18 at 8:04
• "Solution of the Poisson-Boltzmann equation for surface excesses of ions in the diffuse layer at the oxide-electrolyte interface" Borkovec M., and Westall J., 1983. Dec 17 '18 at 8:09
• There remain some small errors due to rounding in the given data so that the condition $\sum z_lc^B_l=0$, which makes the integrand singularity-free, is violated in your data. But that can be cut away, the current (and hopefully final) version of my answer computes this defect and removes its effect. Dec 17 '18 at 10:42

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.