Numerical solution of a singular integral equation

Question

I am looking to approximate the solution u of the following equation using discretization method or any other idea. Is there any way on how to find a numerical solution for it:

u[t]-Integrate[Abs[t - s]^(-1/2)*u[s], {s, 0, 1}] == 1/3 (-2 Sqrt[1 - t]+3t-4 Sqrt[1-t]t-4t^(3/2)) where 0<t<1.

The solution is u[x]=x but I am assuming that I dont know the answer and we need to find approximation for it.

Answer 1

Not an answer, only an idea to solve the problem.

I tried to solve your integral equation iterativ using NestList:

sol = NestList[
Function[fu,
FunctionInterpolation[
 1/3 (-2 Sqrt[1 - t] + 3 t - 4 t Sqrt[1 - t] - 4 t^(3/2)) + 
  NIntegrate[fu[s]/Sqrt[Sqrt[(t - s)^2]] , {s, 0, 1}, 
   Method -> "LocalAdaptive" ], {t, 0, 1 }]
] , 0 &,  (* initial function *)5];

Unfortunately the Picarditeration doesn't converge in your case:

    Plot[Map[#[t] &, sol], {t, 0, 1}

Perhaps you have additional system knowhow to force a convergent iteration?

Answer 2

I will add another method that is not as accurate as method @Roman, but faster. It uses expression describing the integral Integrate[1/Sqrt[Abs[t-s]], {s, 0, 1}]

ker[s_, t_] := If[t > s, -2*Sqrt[t - s], 2*Sqrt[s - t]]

Then everything is as usual

np = 51; points = fun = Table[Null, {np}];
Table[points[[i]] = i/np, {i, np}];
sol = Unique[] & /@ points;

Do[fun[[i]] = 
   1/3 (-2 Sqrt[1 - t] + 3 t - 4 Sqrt[1 - t] t - 4 t^(3/2)) /. 
    t -> points[[i]], {i, np}];

sol1 = sol /. 
   First@Solve[
     Table[sol[[j]] - 
        Sum[.5*(sol[[i]] + 
            sol[[i + 1]])*(ker[points[[i + 1]], points[[j]]] - 
            ker[points[[i]], points[[j]]]), {i, 1, np - 1}] == 
       fun[[j]], {j, 1, np}], sol];

u = Transpose[{points, sol1}];

Show[Plot[t, {t, 0, 1}], ListPlot[u]]

Answer 3

Here's a general solution that works by interpolation. I'll present the method in a very slow way, and we can work on speeding it up later on if desired.

First, we make an ansatz for the function $u(t)$ on the interval $[0,1]$. Here I use a grid of $n+1$ equidistant points and a linear interpolation scheme:

n = 10;
tvalues = Subdivide[n];
uvalues = Unique[] & /@ tvalues;  (* we don't care what these variables are called *)
tupairs = Transpose[{tvalues, uvalues}];
u[t_] = Piecewise@BlockMap[{((t-#[[2,1]])#[[1,2]]-(t-#[[1,1]])#[[2,2]])/(#[[1, 1]]-#[[2, 1]]),
          #[[1,1]]<=t<=#[[2,1]]}&, tupairs, 2, 1]

Check that this interpolation scheme has indeed the values uvalues on the grid points tvalues:

u /@ tvalues == uvalues
(* True *)

Define the integral $\int_0^1 ds\,u(s)/\sqrt{\lvert t-s\rvert}$:

uint[t_] := Integrate[u[s]/Sqrt[Abs[t-s]], {s, 0, 1}]

Evaluate this integral on the same grid of tvalues: here is the slow part of this calculation, and could probably be sped up dramatically,

uintvalues = uint /@ tvalues
(* long output where every element is a linear combination of the uvalues *)

The right-hand side of the integral equation, evaluated on the same grid of tvalues:

f[t_] = 1/3 (-2 Sqrt[1 - t] + 3 t - 4 Sqrt[1 - t] t - 4 t^(3/2));
fvalues = f /@ tvalues
(* long output *)

Solve for the coefficients of $u(t)$: a linear system of equations for the grid values uvalues, found by setting the left and right sides of the integral equation equal at every grid point in tvalues,

solution = tupairs /.
  First@Solve[Thread[uvalues - uintvalues == fvalues] // N, uvalues]

{{0, 5.84947*10^-16}, {1/10, 0.1}, {1/5, 0.2}, {3/10, 0.3}, {2/5, 0.4}, {1/2, 0.5}, {3/5, 0.6}, {7/10, 0.7}, {4/5, 0.8}, {9/10, 0.9}, {1, 1.}}

This confirms your analytic solution $u(t)=t$ but is much more general.

You don't need the // N in the last step if you prefer an analytic solution; however, numerical solution is very much faster.

ListLinePlot[solution, PlotMarkers -> Automatic]

Update: much faster version

To speed up this algorithm, the main point is to speed up the calculation of the uintvalues from the uvalues. Instead of doing piecewise integrals, this calculation can be expressed as a matrix multiplication, uintvalues == X.uvalues, with the matrix X defined as

n = 10;
X = N[4/(3 Sqrt[n])]*
  SparseArray[{{1,1} -> 1.,
               {-1,-1} -> 1.,
               Band[{2,2}, {-2,-2}] -> 2.,
               Band[{2,1}, {-1,1}, {1,0}] ->
                 N@Table[(i-2)^(3/2)-(i-1)^(3/2)+3/2*(i-1)^(1/2), {i,2,n+1}],
               Band[{1,-1}, {-2,-1}, {1,0}] -> N@Reverse@Table[(i-2)^(3/2)-(i-1)^(3/2)+3/2*(i-1)^(1/2), {i,2,n+1}],
               Sequence @@ Table[Band[{1,a}, {1+n-a,n}] -> N[a^(3/2)-2*(a-1)^(3/2)+(a-2)^(3/2)], {a,2,n}],
               Sequence @@ Table[Band[{a+1,2}, {n+1,n+2-a}] -> N[a^(3/2)-2(a-1)^(3/2)+(a-2)^(3/2)], {a,2,n}]},
              {n+1, n+1}] // Normal;

(The coefficients follow from the Piecewise ansatz and analytic integration.)

With this matrix defined, the algorithm becomes simply

tvalues = Subdivide[n];
f[t_] = 1/3 (-2 Sqrt[1 - t] + 3 t - 4 Sqrt[1 - t] t - 4 t^(3/2));
fvalues = f /@ tvalues;
solution = Inverse[IdentityMatrix[n+1] - X].fvalues
ListLinePlot[Transpose[{tvalues, solution}]]

In this way, $n=1000$ grid points can be achieved in a few seconds, most of which is still spent in assembling the X-matrix. The next step would be to write down a faster way of assembling X.