# Numerical solution of a singular integral equation

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.

• Either you need to provide the function u[t], or you need to provide an equation that u[t] satisfies. May 19 '19 at 19:39
• In fact, I need to find this u numerically. May 19 '19 at 19:43
• The above is not valid Mathematica code? In Mathematica a function u(t) is written as u[t] also for u(s) it is u[s] May 19 '19 at 19:44
• Perhaps you meant you need to solve u[t] - Integrate[Abs[t-s]^(-1/2) u[s], {s, 0, 1}] == 0 numerically. If there is no equation, then there is no way to solve for u[t]. May 19 '19 at 19:45
• If an approximation is fine, you could write $u(t)=a_0+a_1t+a_2t^2+\mathcal O(t^3)$ and solve for $a_i$. I get $a_0=0$, $a_1=1$, $a_2=0$, etc, which confirms that $u(t)=t$ is a solution. May 19 '19 at 22:11

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.

• Dear Roman, I couldn't find any word that deserve saying Thank you for the idea you provided. The solution works perfectly for me. We just need to work to speed it up. Also I have another schemes that I need to discuss it with you. Please send me your email. I need someone to work with me on a paid project if you interested. May 20 '19 at 19:37
• You're welcome. Email me at Uncompress["1:eJxTTMoPChZnYGAoys9NzNMrTs7IzUxNcSjNy0xKLNZLzgAAnn0Kkg=="] May 20 '19 at 19:40
• @Roman How did you find the interpolation u[t]? May 21 '19 at 12:06
• @AlexTrounev every piece of the Piecewise function is a two-point Lagrange polynomial. Concretely, if you set f[x_] = y1*(x - x2)/(x1 - x2) + y2*(x - x1)/(x2 - x1), then f[x1] == y1 and f[x2] == y2. In this way the function f passes through the two points {x1,y1} and {x2,y2}. May 21 '19 at 12:27
• @jsxs 1. Correct. All I wanted here is a list of symbols that has the right length. Unique[]& is a pure function with zero arguments (which is perfectly legal), and so the supplied argument in Map is discarded. I agree it's a bit opaque; I could have done uvalues = Table[Unique[], {n+1}] instead for more clarity. May 23 '19 at 11:53

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?

• tried other starting point perhaps? e.g., #& is the exact solution, so a fixed point presumably (is the system stable?) May 20 '19 at 15:54
• I tried this Initial function too but fixpoint isn't stable May 20 '19 at 17:26

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]]