How do we test for infinite series convergence or divergence in MATLAB

I am very new to MATLAB and still learning the ropes. I am looking for a generic method that tests a series for convergence or divergence for infinite series. This is what I tried (I am using MATLAB online):

syms n
term = 1/n^.5;
y = vpa(symsum(term, n, 1, inf));
if isinf(y)
    fprintf('divergent\n');
else
    fprintf('convergent\n');
end         

convergent

But this seems to give an incorrect answer! Can someone please explain how I can test for a generic series convergence in MATLAB?

[Dyuman Joshi]

But the sum of square root of Natural numbers, which corresponds to zeta(1/2), is convergent.

Did you mean to test for the sum of negative 5th power of Natural numbers?

[Someone Someone]

Its 1/(n^.5). If I am correct for any power less than or equal to 1 its divergent!

[Someone Someone]

Moreover, How can sum of square roots of natural numbers be convergent !!??!! With due respect, please don't comment with obviously wrong answers.

[Torsten]

sum_{n=1}^{n=Inf} 1/n^alpha 

converges for alpha > 1 and diverges for alpha <= 1.

[John D'Errico]

@Someone Someone - please be nice to people who are trying to help you. A comment is NOT an answer, merely a comment on your question in an attempt to resolve a point of confusion.

[Someone Someone]

John - Of course. I did not mean of offend anyone the slightest. And the reply was clearly polite too. Still no offense meant.

[Someone Someone]

@Torsten I know that one. That is why the example with a known series. The above code states 'convergent' as the answer (power is .5)! In case it was not clear I wanted to know how do we test for convergence for a generic series in MATLAB reliably (if such obvious series are given wrong answer)? Or am I doing something wrong here?

[Paul]

syms n
term = 1/n^.5;
s = symsum(term, n, 1, inf)

s = 

y = vpa(s)

y = 

Verify statement by @Dyuman Joshi

z = zeta(sym(1)/sym(2))

z = 

vpa(z)

ans = 

So vpa certainly seems to be doing something more than just coming up with a VPA value for s using ordinary VPA arithmetic. And if vpa "knows" the result (if we can call it that) is zeta(1/2), why doesn't symsum?

But this works as expected

s = symsum(n,n,1,inf)

s = 

∞

As opposed to being not evaluated and then being able to use vpa(s) to get zeta(-1) = -1/12

zeta(sym(-1))

ans = 

which seeemed to be all the rage on the internet a few years ago.

[Someone Someone]

@Paul I understand. But my primary concern is finding if a real number infinite series is convergent or not. vpa is doing something more than summation then it might/will not fit the purpose.

Is there any other way I can test if a series is convergent or not?

[Torsten]

@Paul

So it's a bug ?

[Paul]

@Someone Someone

Based on my very limited experience with symsum, I'm not so sure that I would rely on it coming up with a finite, closed-form expression as a test for convergence. That is, I'm pretty sure that I've seen symsum not find a closed-form expression when I know one exists. But then again, my use of symsum typically involves expressions that contain variables in addition to the summing variable, which is different than what you're trying to do. YMMV.

@Torsten

I don't know whether or not to think this is a bug. It's certainly unexpected and undocumented behavior of vpa IMO, and I'd like to know why vpa seems to be so clever (and if it will produce other "unexpected" results). vpasum also seems to be more clever than symsum

syms n
term = 1/n^.5;
s = vpasum(term, n, 1, inf)

s = 

Even though

vpasum(term,n,1,1e5)

ans = 

630.99675866237874808482661955721

To be clear, I don't like how vpa and vpasum work for this problem.

[Torsten]

Really strange.

syms n
s = symsum(1/n^sym('0.5'),n,1,Inf)

s = 

s = symsum(1/n^0.5,n,1,Inf)

s = 

syms alpha 
s = symsum(1/n^alpha,n,1,Inf)

s = 

assume(real(alpha)<=1)
s = symsum(1/n^alpha,n,1,Inf)

s = 

[Sam Chak]

Indeed, very peculiar!

syms n

Test #1: integer 0

s = symsum(1/n^sym('0'), n, 1, Inf)

s = 

∞

z = zeta(sym('0'))

z = 

Test #2: decimal form

s = symsum(1/n^sym('0.5'), n, 1, Inf)

s = 

z = zeta(sym('0.5'))

z = 

Test #3: fraction form

s = symsum(1/n^sym('1/2'), n, 1, Inf)

s = 

s = vpa(s)

s = 

z = zeta(sym('1/2'))

z = 

z = vpa(z)

z = 

Test #4: integer 1

s = symsum(1/n^sym('1'), n, 1, Inf)

s = 

∞

z = zeta(sym('1'))

z = 

∞

[Paul]

@Torsten

I only have a guess explanation for your examples. First define alpha.

syms n integer
alpha1 = sym(1)/sym(2)

alpha1 = 

alpha2 = sym('0.5')

alpha2 = 

0.5

We see that alpha1 and alpha2 aparently have two different internal represnetations, even though they are both of type sym. alpha1 looks like a symbolic reprensentation of "one half" while alpha2 looks like a VPA form of the decimal representation of the number 1/2. Using the former

S1 = symsum(1/n^alpha1,n,1,inf)

S1 = 

Using the latter

S2 = symsum(1/n^alpha2,n,1,inf)

S2 = 

My guess, and it's only a guess, is that because alpha2 is a VPA variable, the call to symsum actually results in a call to vpasum under the hood

V2 = vpasum(1/n^alpha2,n,1,inf)

V2 = 

And with alpha1 we get the same result, presumably because vpa(alpha1) == alpha2.

isAlways(vpa(alpha1) == alpha2)

ans = logical
   1

V1 = vpasum(1/n^alpha1,n,1,inf)

V1 = 

But the situation gets messier, i.e., inconsistent, if we treat alpha as a parameter and delay evaluation of S, either as a symfun with a parameter or as a sym expression followed by subs.

Define S as a symfun

syms alpha real
S(alpha) = symsum(1/n^alpha,n,1,inf)

S(alpha) = 

Note that S(alpha) is defined using piecewise and the otherwiseVal is not specified.

Evaluate S at alpha1.

S(alpha1)

ans = 

NaN

It returns NaN, which is consistent with its piecewise definition, but is inconsistent with S1.

Evaluate S at alpha2.

S(alpha2)

ans = 

NaN

Consistent with S(alpha1), which is good, but inconsistent with S2.

Using subs with sym expressions yields the same results

S = symsum(1/n^alpha,n,1,inf)

S = 

subs(S,[alpha1 alpha2])

ans = 

Define V(alpha) using vpasum

V(alpha) = vpasum(1/n^alpha,n,1,inf)

V(alpha) = 

Evaluate and show consistency with V2

V(alpha1)

ans = 

V(alpha2)

ans = 

So vpasum appears to be self-consistent in these cases, even if the result is sketchy IMO.

OTOH, symsum seems to be a mess.

The more I think about this, the more surprised I am that, IF zeta(1/2) is an acceptable answer for this problem, then symsum doesn't yield zeta(1/2) but vpasum does. IOW, if anything I'd expect symsum would have all of the fancy rules and tables for symbolic sums, whereas vpasum would just add up the terms until meeting criteria for either convergence or divergence.

Presumbaly this call using vpa is actually implemented using vpasum.

term = 1/n^sym(.5)

term = 

s = symsum(term,n,1,inf);
vpa(s)

ans = 

Answer 1

[John D'Errico]

I suppose with some effort I could find a series where convergence is not known. But it is early in the morning, and mental effort is just asking too much. :)

I would point out that the series you show, thus

syms n
term = 1/n^.5;
y = symsum(term, n, 1, inf)

y = 

In fact, does diverge This is easy enough to show, since we can see it dominates the similar series symsum(1/n,1,inf) term by term. And we do know the harmonic series diverges.

But then you tried:

vpa(y)

ans = 

which is a finite value. In fact, it is the same value that @Dyuman Joshi pointed out, as

zeta(1/2)

ans = -1.4604

It looks like VPA merely looks at that sum, and recognizes that zeta(s) is the same as that sum for general s, as long as real(s)>1. And so it returns zeta(1/2). Since zera(1/2) is finite, your test failed to produce the result you expected.

That is what happens, though I cannot really excuse VPA for having made that choice. Arguably, something in that chain should produce a result that indicates divergence. Feel free to send it into tech zupport and suggest that VPA might have known better.

However, is there some simple test one can apply that will always robustly indicate convergence or not, for some completely general series? I don't know of one.