Computing inverse discrete Stieltjes transform?

Question

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 thef(i)`. (3) You can isolate poles by rectangular subdivision, and use that to find all in some given rectangular range. — Daniel Lichtblau, Jun 27 '23 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 — Yaroslav Bulatov, Jun 27 '23 at 18:36

score 3 · Accepted Answer · answered Jun 27 '23 at 23: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 — Yaroslav Bulatov, Jun 28 '23 at 00: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. — Daniel Lichtblau, Jun 28 '23 at 13:12

Computing inverse discrete Stieltjes transform?

1 Answers1