# Computing inverse discrete Stieltjes transform?

Is there a way to invert the Stieltjes transform of a discrete distribution in Mathematica? IE, given $$g(u)=\sum_{i=1}^\infty f(i) \frac{1}{u-f(i)}$$, obtain $$f(i)$$

ClearAll["Global*"];
f = x^-2;
g = Sum[ f/(u - f), {x, 1, \[Infinity]}] (* (Sqrt[u] - \[Pi] Cot[\[Pi]/Sqrt[u]])/(2 Sqrt[u]) *)


For a continuous distribution, the recipe is given by JimB here

Edit

it seems possible to invert by following the recipe here if there were a way to get all the singularities of $$g$$

ClearAll["Global*"];
f = x^-2;
g = Sum[f/(u - f), {x,
1, \[Infinity]}] (*(Sqrt[u]-\[Pi] Cot[\[Pi]/Sqrt[u]])/(2 Sqrt[u])*)
singularities = Table[1/i^2, {i, 1, 10}];
Print["f=", f];
Print["g=", g];
Table[Residue[g, {u, u0}]/u0, {u0, singularities}]

• (1) You have infinitely many poles so they go to (complex) infinity and you won't be able to get all of them numerically. (2) As noted, you can get your values f(i) if you know the singularities of g. But that's kind of obvious already, because the singularities are exactly the f(i). (3) You can isolate poles by rectangular subdivision, and use that to find all in some given rectangular range. Commented Jun 27, 2023 at 18:21
• @DanielLichtblau I was curious if one could go from $\frac{\sqrt{u}-\pi \cot \left(\frac{\pi }{\sqrt{u}}\right)}{2 \sqrt{u}}$ to $x^{-2}$ using symbolic manipulation, analogous how JimB could do it here in the continuous case Commented Jun 27, 2023 at 18:36

This is (at best) a partial answer. A bit of calculus shows you can substitute x for Sqrt[u] to get x - Pi*Cot[Pi/x] (the change of variable differential conveniently cancels the 2*x denominator). The poles are clearly at reciprocals of integers. We will compute some residues.

ee = 2*x*(x - Pi*Cot[Pi/x])/(2*x);
etab = Chop[Table[
rep = 1/j + .001*Exp[2*Pi*I*t];
ff = ee*D[rep, t] /. x -> rep;
1/(2*Pi*I)*NIntegrate[ff, {t, 0, 1}, PrecisionGoal -> 12],
{j, 1, 20}]]

(* Out[89]= {1., 0.25, 0.111111, 0.0625, 0.04, 0.0277778, 0.0204082, \
0.015625, 0.0123457, 0.01, 0.00826446, 0.00694444, 0.00591716, \
0.00510204, 0.00444444, 0.00390625, 0.00346021, 0.00308642, \
0.00277008, 0.0025} *)


Now rationalize them to see what we have.

Rationalize[etab]

(* Out[92]= {1, 1/4, 1/9, 1/16, 1/25, 1/36, 1/49, 1/64, 1/81, 1/100, \
1/121, 1/144, 1/169, 1/196, 1/225, 1/256, 1/289, 1/324, 1/361, 1/400} *)


This does give the values for x^(-2). Not sure it was quite what you want but maybe it will give some ideas?

• I'm curious why you use NIntegrate instead of built-in Residue` function Commented Jun 28, 2023 at 0:19
• Didn’t think of it and was short on time. It occurs to me that it’s easier just to note that the poles of you function of u occur at reciprocals of squared integers. But it’s nice to see the CIF working as expected. Commented Jun 28, 2023 at 13:12