# Computing numerically infinite sum of some double series

Let's consider the series: $$F(t) = \sum\limits_{n=0}^{\infty}\sum\limits_{k=0}^{\infty} \frac{(-b)^k(-a)^n\binom{n+k}{k} t^{2n+k(2-\alpha)}}{\Gamma(2n+k(2-\alpha)+2)}$$

where $$a,b$$ are positive reals, $$0 < \alpha < 1$$, and $$\Gamma(z)$$ is the gamma function.

I would like to numerically compute the sum $$F$$. I replace the doubly infinite sum by a double finite sum in the following way:

$$\sum\limits_{n=0}^{\infty}\sum\limits_{k=0}^{\infty}\left(\ldots\right) \approx \sum\limits_{i=0}^{N}\sum\limits_{{n+k = i}\atop{0\leq n,k\leq i}} \left(\ldots\right) \ \ \ (\star)$$

where bound $$N$$ is set in such a way, that $$\left|\frac{S_{N+1}-S_{N}}{S_{N+1}}\right|<\varepsilon$$ ($$S_j$$ - sum from right hand path of the ($$\star$$) for $$N=j$$; $$\varepsilon$$ - some small number, for example $$\varepsilon = 10^{-6}$$).

Unfortunately there is something wrong with my approach. For example for $$a=1$$, $$b=1$$, $$\alpha=0.5$$ and $$0\leq t \leq 30$$ I get following graph: Definitely, there is something wrong for $$t\geq 25$$.

My question: How to compute approximation of $$F(t)$$ correctly also for $$t\geq 25$$?

Code (in this way I define function $$F(t)$$):

F[t_, α_] := Block[{convergence, N1, S, SN, tab, re, i},
N1 = 0;
S = 1;
convergence = False;
While[! convergence && N1 <= 1000,
N1++;
(* ponizej k = i, n = N1-i *)
tab = Table[((-v0)^i (-w0^2)^(N1 - i) Binomial[N1, i] t^(
2 (N1 - i) + i (2 - α)))/
Gamma[2 (N1 - i) + i (2 - α) + 2], {i, 0, N1}];
SN = S + Sum[tab[[i]], {i, 1, Length[tab]}];
re = Abs[(S - SN)/S];
If[re < 10^-6, convergence = True];
S = SN;
];
Return[{S, N1, re}];
]

• Sorry for such strange formating but I could not add this post due to error "Your post appears to contain code that is not properly formatted as code. Please indent all code by 4 spaces using the code toolbar button or the CTRL+K keyboard shortcut. For more editing help, click the [?] toolbar icon." (it was some problem with LaTeX code) Sep 24, 2018 at 10:54
• You can probably do one summation symbolically, then your problem should be much easier. Sep 24, 2018 at 12:03

This is a very rough solution I just hacked out; I might try to refine it later.

I had previously talked about Cuyt's transformation of a double series as an aid in evaluating them numerically in this answer. Applied to this case, we have

Plot[With[{a = 1, b = 1, t = SetPrecision[tx, 25], α = 1/2},
NumericalMathNSequenceLimit[Accumulate[Table[Sum[
Function[{n, k}, (-a)^n (-b)^k Binomial[n + k, k] t^(2 n + k (2 - α))/
Gamma[2 n + k (2 - α) + 2]] @@ v,
{v, FrobeniusSolve[{1, 1}, k]}], {k, 0, 29}]]]], {tx, 0, 30},
MaxRecursion -> 0, PlotPoints -> 75, PlotRange -> All] where I use NumericalMathNSequenceLimit[] for performing the Shanks transformation on the result of Cuyt's transformation.

The method is a bit finicky numerically, which is why I needed the use of SetPrecision[]. I also very arbitrarily chose the cutoff 29 in Cuyt's transformation; for other parameter choices, less or more terms might be necessary. All of this will need to be experimented with.

• Thank you for very nice solution of my problem (unfortunately I use version <11.2 of Mathematica, so function NumericalMathNSequenceLimit[] in unavailable for me) Sep 24, 2018 at 20:23
• Replace that with SequenceLimit[]`, then. Sep 24, 2018 at 20:29
• It works, thank you :) Sep 26, 2018 at 19:37