Fredholm Integral Equation of the 2nd Kind with a Singular Difference Kernel

Question

After reviewing the literature, I could not find an analytic solution to the equation
$$\int^{1}_{0}dx\frac{f(x)}{|x-y|^{2/3}}=cf(y)$$
for $f(y)$, where $c$ is a constant and $y\in[0,1]$.

I'm thus trying to find a numerical solution for it using Mathematica.

I've (unsuccessfully) tried using the approaches shown in "iteratively solve integral equation" and "Fredholm integral equation of the second kind with kernel containing Bessel and Struve functions" but was unable to apply the methods discussed there to my problem.

Thanks in advance

Answer 1

This is an eigenvalue problem.

Let's apply a Galerkin scheme: We fix a finite dimensional space of functions, pick a basis $u_0, u_2,\dotsc,u_{n}$ and define the matrices

$$A_{ij} = \int_0^1 \!\!\!\!\int_0^1 u_i(x) \, k(x,y) \, u_j(y) \, \operatorname{d}\! x \, \operatorname{d}\! y$$ and $$M_{ij} = \int_0^1 \!\!\!\!\int_0^1 u_i(x) \, u_j(y) \, \operatorname{d}\! x \, \operatorname{d}\! y,$$ where $k(x,y) = \frac{1}{|x-y|^{2/3}}$. Afterwards, we we solve the generalized eigensystem

$$ (A - \lambda) \, f = 0.$$

Depending on the basis you chose, you will get different methods. I guess trigonometric polynomials or Chebyshev polynomials should work fine with superior accuracy. I'll make it quick and dirty and use continuous piecewise-linear functions with respect to a grid $0 = t_0 < t_1 < \dots < t_n = 1$.

This is how we can assemble the matrices:

n = 1000;
xlist = Subdivide[0., 1., n];

M = SparseArray[{
    Band[{1, 1}] -> 2./(3. n),
    Band[{2, 1}] -> 1./(6. n),
    Band[{1, 2}] -> 1./(6. n)
    }, {n + 1, n + 1}, 0.];
M[[1, 1]] = 1./(3. n);
M[[n + 1, n + 1]] = 1./(3. n);

kernel = 1./(DistanceMatrix[xlist, xlist]^(2/3) + IdentityMatrix[n + 1, SparseArray]) 
  - IdentityMatrix[n + 1, SparseArray];
A = M.kernel.M;

We could use

{λ,U}=Eigensystem[{A,M}];

to solve the eigenvalue problem, but it takes rather long. Faster is this:

{λ, U} = Eigensystem[kernel.M];

A plot of the eigenfunctions for the eight highest eigenvalues can be obtained with

ListLinePlot[U[[1 ;; 8]]]

Answer 2

As mentioned by Henrik, this is an eigenvalue problem. Since Mathematica doesn't have a built-in eigenvalue problem solver for integral equation, we need to discretize the equation to matrix form by ourselves and solve the resulting system with Eigensystem. Henrik has discretize the integral based on Galerkin scheme. Here I'd like to add a solution discretizing the integral with Gauss-Legendre quadrature formula:

eps = 10^-2;
points = 400;
integrand[x_] = (1/(eps + Abs[x - y]^(2/3))) f[x];
domain = {0, 1};
{nodes, weights} = Most[NIntegrate`GaussRuleData[points, MachinePrecision]];
midgrid = Rescale[nodes, {0, 1}, domain];
grid = Flatten[{domain[[1]], midgrid, domain[[-1]]}];
int = -Subtract @@ domain weights.Map[integrand, midgrid];
{b, m} = CoefficientArrays[int, f /@ grid];
mat = Table[m, {y, grid}];
{val, vec} = Eigensystem[mat, 8];
ListLinePlot[vec[[;; 8]], PlotRange -> All]

Error check:

With[{i = 1}, 
 With[{f = ListInterpolation[vec[[i]], domain], c = val[[i]]}, 
   ListLinePlot[
    Table[NIntegrate[f[x]/Abs[x - y]^(2/3), {x, 0, 1}, 
       Method -> {Automatic, "SymbolicProcessing" -> 0}] - c f[y], {y, 0, 1, 1/100}], 
    PlotRange -> All]] // Quiet]

Remark

1. Introduction for NIntegrate`GaussRuleData can be found here.

2. The error of solution is relatively large, but it decreases as points increases.

3. The singularity is circumvented by adding a small value eps to the denominator of integrand, which is rough. I believe this is the root of inaccuracy of solution. Sadly I can't think out a better way to deal with the singularity.

4. Notice the definition of mat is somewhat opportunistic. It needs to be more sophisticated if the integrand involves f[y].