# 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. – Lutz Lehmann Dec 14 '18 at 21:57
• Could you add an image of the formula in the original paper for reference? – Lutz Lehmann 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. – Daniel 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. – Daniel 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. – Lutz Lehmann 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.