Limit of partial sums?

Question

So I'd like to calculate limit of partial sums, where `a_n` is `sp.Lambda(x, 1/x * sp.factorial2(2*x-1)/sp.factorial2(2*x))`. I tried this: sp.limit(sp.summation(a(i), (i,1,n)), n, sp.oo) But as I start it, it starts calculating but never ends, so I don't get an answer. Interesting fact is that if I try to calculate infinite sum, I'll get an answer: sp.summation(a(i), (i,1,sp.oo)) >>>2log(2) I can even create a sequence of partial sums: s = sp.SeqFormula(sp.summation(a(i), (i,1,n)), (n, 1, sp.oo)) Which, for example in Maple, evaluates the limit kinda. So is there a way to calculate such limit in sympy? Actually I'd like to get an answer similar to summation: `2log(2)`, but using the limit of partial sums.

Answer 1

Short version first: sympy can already give you this answer, just not by the route you are taking. The infinite sum closes. The limit of the partial sums does not, and I can show you exactly where it breaks.

Before anything else, your term simplifies a lot. Since (2n−1)!!/(2n)!! is the same as C(2n,n)/4ⁿ, your a_n is

a_n = C(2n,n) / (n · 4ⁿ)

which is a much friendlier object than a pair of double factorials. Either form works below, but it helps to know that is what you are looking at.

The infinite sum comes back immediately:

import sympy as sp
n = sp.symbols('n', positive=True, integer=True)
a = sp.factorial2(2*n-1)/(n*sp.factorial2(2*n))

sp.Sum(a, (n, 1, sp.oo)).doit()      # log(4)

log(4) is your 2log(2), so you already have the number.

Now the part you actually asked about. It is worth separating two steps here, because only one of them is broken. The symbolic partial sum works fine:

N = sp.symbols('N', positive=True, integer=True)
sp.Sum(a, (n, 1, N)).doit()

and it returns

2*log(2) - gamma(N + 3/2)*hyper((1, N + 1, N + 3/2), (N + 2, N + 2), 1)
           / (sqrt(pi)*(N + 1)*gamma(N + 2))

I checked that against simply adding the terms up at N = 1, 2, 5, 10 and 40, and it agrees to twenty digits. So sympy can find your partial sum.

What it cannot do is take the limit of that. Call sp.limit on it and you get NotImplementedError, with a message about not knowing how to compute the mrv of (1, N + 1, N + 3/2). Its limit algorithm does not handle a hypergeometric whose parameters themselves depend on the variable you are sending to infinity. That is the real wall, and it is not something you can configure around. Your original code was hanging rather than erroring because sp.summation inside sp.limit kept re-deriving that same partial sum.

The good news is that the closed form hands you the answer anyway. Look at its shape: it is literally 2*log(2) minus a tail. So the moment you know the tail goes to zero, the limit is 2log(2) by inspection, and you do not need sympy to finish that step.

If you would rather see the tail vanish than take it on faith, use the size of the terms. For large n, C(2n,n)/4ⁿ is about 1/√(πn), so

a_n ≈ 1 / (√π · n^3/2)

and summing that from N onward gives a tail of roughly 2/√(πN), which goes to zero. I checked it numerically and the agreement gets very tight:

    N        2log2 - S_N      2/sqrt(pi N)     ratio
   100       0.11250988       0.11283792      0.99709
   400       0.056377853      0.056418958     0.99927
 10000       0.011283463      0.011283792     0.99997
 40000       0.0056418547     0.0056418958    0.99999

That table also quietly explains why you should not chase this numerically. The tail shrinks like one over the square root of N, so to pin down even three decimal places you would need to add on the order of five million terms. Your instinct to do it symbolically was the right one, the partial sum just is not the symbolic path that closes.

So in practice: use sp.Sum(a, (n, 1, sp.oo)).doit() if you want the number, use sp.Sum(a, (n, 1, N)).doit() if you genuinely want the partial sum as a formula, and read 2*log(2) off the front of that formula if what you want is the limit with a reason attached rather than a black box result.