Numerical solution

In the most general form the integral equations we have to solve is written as

$$f(x)+\int\limits_{-Y}^Y K(x-y)f(y)\,dy = g(x),$$ (10.21)

where $f(x)$ is so far unknown solution, $K(x)$ is the kernel, $Y$ defines the integration limit, $g(x)$ defines the normalization (in LL case it is constant). The function $f(x)$ enters twice: once inside the integral and second time outside, this can be remedied inserting the $\delta$-function:

$$\int\limits_{-Y}^Y(\delta(x-y)+K(x-y))f(y)\,dy = g(x),$$ (10.22)

Now we do discretization with spacing $\Delta x$. The equation (A.22) now can be expressed in the matrix form:

$$(I+K\Lambda) \vec f \Delta x = \vec g,$$ (10.23)

here $I$ stands for a unity matrix and the diagonal matrix $\Lambda$ is defined by the integration method. Now the vector $f$ is obtained by multiplication of the inverse matrix on $\vec g$:

$$\vec f = \frac{1}{\Delta x} (I+K\Lambda)^{-1} \vec g,$$ (10.24)

For a uniform grid very good precision is achieved using the Simpson method. The matrix $\Lambda$ in this case is defined as $\Lambda = \mbox{diag}\{\frac13,\frac13,\frac43,\frac23,\frac43,..., \frac43,\frac23,\frac43,\frac13\}$. The residual term of the integration is very small and can be estimated as $I_{err} = \max f^{(4)}(x)\frac{(\Delta x)^5}{2880}$ and the error in the energy (which is defined by integrating the solution $f(x)$ with the weight proportional to $x^2$) is proportional to the spacing $\Delta x$ to the forth power.