Numerical solution of a singular integral equation
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?
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
.
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]]