Using Numeric Inverse Stieltjes Transform from here:
ClearAll["`*"]; Remove["`*"];
f[s_] := (Sqrt[2] (\[Pi] - 2 ArcTan[Sqrt[s]/Sqrt[2]]))/(2 Sqrt[s] +
Sqrt[2] s ArcTan[Sqrt[2]/Sqrt[s]])
g1[x_] := 1/\[Pi]^2*NIntegrate[1/(s - x)*f[s], {s, -Infinity, x, Infinity},
Method -> "PrincipalValue", WorkingPrecision -> 20,
MaxRecursion -> 20](*Method One*)
A = ListLinePlot[Table[{x, Re[g1[x]]}, {x, -3 + 10^-15, 3, 1/50}],
PlotStyle -> Red] // Quiet
g2[x_] := -(2/Pi^2)*NIntegrate[(f[x - s] - f[x + s])/(2 s), {s, 10^-15,
Infinity}](*Method Two*)
B = ListLinePlot[Table[{x, Re[g2[x]]}, {x, -3 + 10^-15, 3, 1/50}],
PlotStyle -> {Black, Dashed}] // Quiet
Show[{A, B}](*Two Methods on one Plot*)
But that doesn't plot look like @JimB answer!.
EDITED 17.10.2024
From comments @Yaroslav Bulatov Harder example, to invert:
distY = TransformedDistribution[1/z^2,
z \[Distributed] ZipfDistribution[1]]; Expectation[1/(s - i),i\[Distributed] distY]
(*(3 (Sqrt[s] - \[Pi] Cot[\[Pi]/Sqrt[s]]))/(\[Pi]^2 Sqrt[s])*)
If I use Inverse Stieltjes Transform to:$\frac{3 \left(\sqrt{s}-\pi \cot \left(\frac{\pi }{\sqrt{s}}\right)\right)}{\pi ^2 \sqrt{s}}$ i should get:$\begin{cases}
\frac{6 \left(\psi ^{(1)}\left(\left\lceil \frac{1}{\sqrt{x}}\right\rceil \right)-\psi ^{(1)}\left(\left\lfloor
\frac{1}{\sqrt{x}}\right\rfloor +1\right)\right)}{\pi ^2} & 0<x\leq 1 \\
0 & \text{True}
\end{cases}$.
Let try:
PDF[TransformedDistribution[1/z^2,
z \[Distributed] ZipfDistribution[1]], x]
(*Piecewise[{{(6*(PolyGamma[1, Ceiling[1/Sqrt[x]]] -
PolyGamma[1, 1 + Floor[1/Sqrt[x]]]))/Pi^2,
Inequality[0, Less, x, LessEqual, 1]}}, 0]*)
Plot[%, {x, -2, 2}, PlotRange -> Full]
(*I get Zero for all range*)
Using: $$\sum _{n=1}^{\infty } \frac{6}{\pi ^2 \left(-1+n^2 s\right)}=\frac{3 \left(\sqrt{s}-\pi \cot \left(\frac{\pi }{\sqrt{s}}\right)\right)}{\pi ^2 \sqrt{s}}$$ Then:
f[s_] := 6/(\[Pi]^2 (-1 + n^2 s));
Sum[InverseLaplaceTransform[
InverseLaplaceTransform[f[s], s, t] /. t -> -t, t, x], {n, 1,
Infinity}](*Inverse Stieltjes Transform -method 1*)
(*Sum[(6*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}]*)
Sum[InverseLaplaceTransform[FullSimplify[
InverseMellinTransform[f[s], s, w, GenerateConditions -> False]
/. w -> Exp[t], Assumptions -> t > 0], t, x], {n, 1, Infinity}]
(*Inverse Stieltjes Transform -method 2*)
(* 0 *)
Sum[Simplify[InverseMellinTransform[
Simplify[
InverseMellinTransform[f[s], s, Exp[t],
GenerateConditions -> False], Assumptions -> t > 0], t, Exp[-x],
GenerateConditions -> False], Assumptions -> x > 0], {n, 1,
Infinity}](*Inverse Stieltjes Transform -method 3*)
(* 0 from InverseMellinTransform *)
Sum[FullSimplify[-1/Pi*InverseFourierTransform[-I Sign[z] FourierTransform[f[s], s, z,
FourierParameters -> {1, -1}], z, x,
FourierParameters -> {1, -1}], Assumptions -> x > 0], {n, 1, Infinity}]
(*Inverse Stieltjes Transform -method 4*)
(*Sum[(12*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}] two times bigger than method 1*)
Sum[Simplify[1/Pi*InverseFourierTransform[
I*(2*HeavisideTheta[z] - 1)*
FourierTransform[f[s], s, z, FourierParameters -> {1, -1},
Assumptions -> n > 0], z, x, FourierParameters -> {1, -1}],
Assumptions -> x > 0], {n, 1, Infinity}]
(*Inverse Stieltjes Transform -method 5*)
(*Sum[(12*DiracDelta[-(1/n^2) + x])/(n^2*Pi^2), {n, 1, Infinity}] two times bigger than method 1*)
Sum[1/\[Pi]^2*Integrate[1/(s - x)*f[s], {x, -Infinity, Infinity},
PrincipalValue -> True, Assumptions -> x > 0], {n, 1,
Infinity}](*Inverse Stieltjes Transform -method 6*)
(*Integral does not converge*)
Sum[-(2/Pi^2)*Limit[Integrate[(f[x - s] - f[x + s])/(2 s), {s, eps, Infinity},
Assumptions -> {eps > 0, x > 0}, GenerateConditions -> False],
eps -> 0, Assumptions -> x > 0, Direction -> -1], {n, 1,
Infinity}](*Inverse Stieltjes Transform -method 7*)
(*Sum[-((6*(Log[1/n^2 - x] - Log[-(1/n^2) + x]))/(Pi^4*(-1 + n^2*x))), {n, 1, Infinity}]*)
h[x_] := Sum[-((6 (Log[1/n^2 - x] - Log[-(1/n^2) + x]))/(\[Pi]^4 (-1 + n^2 x))), {n, 1, 100}]
ReImPlot[h[x], {x, -2, 2}](*For Real part h[x] is 0 *)
Computing with Numerics:
R[s_] := (3 (Sqrt[s] - \[Pi] Cot[\[Pi]/Sqrt[s]]))/(\[Pi]^2 Sqrt[s]);
g1[x_] := 1/\[Pi]^2*
NIntegrate[1/(s - x)*R[s], {s, -Infinity, x, Infinity},
Method -> "PrincipalValue", WorkingPrecision -> 20,
MaxRecursion -> 20](*Method One*)
dane = Table[{x, Re[g1[x]]}, {x, -2, 2, 1/20}] // Quiet;
ListLinePlot[dane, PlotStyle -> Red, PlotRange -> All] // Quiet
(*Random noise*)
g2[x_] := -(2/Pi^2)*NIntegrate[(R[x - s] - R[x + s])/(2 s), {s, 10^-15, Infinity},
WorkingPrecision -> 20](*Method Two*)
B = ListLinePlot[Table[{x, Re[g2[x]]}, {x, -2, 2, 1/20}]] // Quiet
(*Random noise*)
Could say that the solution is distribution nature that's why numeric's methods fail.