I'd like to write a small library which uses the boundary element method to compute the charge distributions on bodies. This requires me to solve the following integrals $$\int_{0}^{1}{\log{\sqrt{[x_i-x(t)]^2+[y_i-y(t)]^2}}\,dt}$$ and $$\int_{0}^{1}{\frac{[x_i-(t)]\,\vec{e}_x+[y_i-y(t)]\,\vec{e}_y}{\left\{[x_i-x(t)]^2+[y_i-y(t)]^2\right\}^2}\cdot\vec{n}_i\,dt}.$$ Sometimes it holds that $x(t)=x_i\ t\in\{0,1\}$ or $y(t)=y_i,\ t\in\{0,1\}$. However, it can be assumed that $x(t)\ne x_i\ \forall\,t\in(0,1)$ and $y(t)\ne y_i\ \forall\,t\in(0,1)$.
I tried to solve these integrals with gsl_integration_qags. Unfortunately I get the following error message:
gsl: qags.c:563: ERROR: integral is divergent, or slowly convergent
I also tried to use gsl_integration_qagp and include the endpoint in the list of singular points. According to the GSL documentation, singular points given to gsl_integration_qagp must lie strictly within the integration interval and so I shouldn't be surprised to get an error message which was:
gsl: qagp.c:494: ERROR: bad integrand behavior found in the integration interval
Can you please recommend me a suitable C/C++ library (function) which is able to solve the integrals mentioned above. It would also be really helpful if someone could shed some light on what kind of integration schemes are suitable for such kind of integrals.
The first integral may be handled using analytic techniques. Notice that:
$$\log\sqrt{\Delta x^2+\Delta y^2}=\log\frac{\sqrt{\Delta x^2+\Delta y^2}}t+\log(t)$$
which allows the first logarithm to avoid a singularity at $t=0$. Similarly, if a singularity occurs at $t=1$, then divide by $1-t$. The leftover logarithm may then be computed as
$$\int_0^1\log(t)~\mathrm dt=\int_0^1\log(1-t)~\mathrm dt=-1$$
In the event that a singularity still remains, you can apply the above multiple times.
The second integral does not look like it can be handled so nicely. In general, singularities when integrating may be handled using adaptive step sizes. The idea is that you sample points sparsely where the integrand is well-behaved and you sample points densely where the integrand is ill-behaved (e.g. singularity). For example, instead of sampling points every $\Delta t=0.01$, you may choose to set
$\min|\Delta t|=0.00001$ at the singularity.
$\max|\Delta f(t)|=0.1$ on the singularity.
$\max|\Delta t|=0.01$ away from the singularity.
so the procedure would be
Set $\Delta t:=\pm0.01$ ($\pm$ for going forwards or backwards).
Compute $f(t+\Delta t)$.
Compute $|\Delta f(t)|$. If it is greater than $0.1$, set $\Delta t:=\pm0.1|\Delta t/\Delta f(t)|$.
If $\Delta t>0.00001$, continue with the chosen method. Otherwise set $\Delta t=0.00001$ and use one left/right Riemann sum (don't use the point $f(t+\Delta t)$ because it may be troublesome).
Note that you have to start at a point which is not a singularity. Note the above is just a simple example of the concept, a more complicated formulation is given on Wikipedia.