MATLAB Answers

How to Get Answers in Non-PolyGamma Form?

2 views (last 30 days)
Mudtimud
Mudtimud on 22 Feb 2011
Answered: John D'Errico on 27 Mar 2020
>> syms k n; symsum(1/k - 1/(k + 1), 1, n)
The ideal answer should be:
ans = -1/(n+1)+1
However, I keep getting my answer in this form:
ans =
psi(n + 1) - psi(n + 2) + 1
which i don't understand since polygamma is currently beyond me. I can, however, achieve the simplified answer with:
>> simplify(ans)
ans =
1 - 1/(n + 1)
But not when:
>> symsum(1/k^2, 1, n)
ans =
pi^2/6 - psi(n + 1, 1)
>> simplify(ans)
ans =
pi^2/6 - psi(n + 1, 1)
How can i get my answer in a non-polygamma form?

  0 Comments

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 27 Mar 2020
(I'm sorry that my answer is coming 9 years too late to be of any significant value to you personally. But someone else might find this question of interest at some point in time.)
You presume an answer exists, in non-polygamma form. That need not always be true. Not all things you can write down and pose for solution need always have a simple answer. Sometimes the only simple answer is just in the form of a special function. And sometimes, the answer is not even simply writable using special functions. You get lucky here, in that a special function does exist to yield any solution at all.
In the first example, it does indeed have a trivial solution, because the series always compresses down. That is, we can expand the finite sum:
symsum(1/k - 1/(k + 1), 1, n)
into a pair of sums:
symsum(1/k, 1, n) - symsum(1/(k + 1), 1, n)
Then, inside the second sum, a simple change of variables is sufficient, to yield
symsum(1/k, 1, n) - symsum(1/k, 2, n+1)
But, now we can consider that some of those terms in the two sums are identical except for sign, and thus cancel out. The result now collapses into two terms only, this the first term from the first sum, and the last term from the second sum, yielding the result:
1/1 - 1/(n+1)
This is a common trick in these things. The intermediate result that symsum returned was unfortunately not as simple as you would like, but this is often true when using symbolic tools. They can generate a result, but it need not always be in as pretty a form as you would prefer. That often requires further artifice to convince the tools to reduce the result to something in a more elegant form. The problem is, a computer has no idea what would look elegant to you, perhaps as opposed to someone else. To paraphrase an already oft-used phrase: elegance is in the eyes of the beholder.
Sadly, you do not have that simple artifice for the direct sum of inverse squares. We get no such easy compression, although we do get a result that makes some mathematical sense.
The result you were given was:
symsum(1/k^2, 1, n) == pi^2/6 - psi(1,n + 1)
(By the way, when you wrote your question, you apparently made a typo, in that you swapped the arguments of psi. And, while your question is now 9 years old as I write this, I doubt they changed the functionality of psi in MATLAB even in that span of time.)
syms k n
>> symsum(1/k^2,[1,n])
ans =
pi^2/6 - psi(1, n + 1)
Anyway, what does this mean? Well the first term is easy to recognize, at least it is if we understand that the sum of squares of integers from 1 to inf is pi^2/6. That goes back to Euler, who derived that (very pretty) result for the rest of us to appreciate. In fact, just look online and you will find all sorts of pages explaining the result.
What then does the second term imply? It just tells us that if you want to compute your desired FINITE sum, the simplest way to represent it is to subtract the sum of the squares of integers from (n+1) to infinity, from a corresponding sum that extends from 1 to infinity. Mush as in your first problem, now we have two INFINITE but overlapping sums. Again, we can choose to cancel the corresponding terms, except that here the terms that cancel out are the terms in the finite sum. So the difference between the two infinite series is just the desired finite sum.
symsum(1/k^2, 1, n) == symsum(1/k^2, 1, inf) - symsum(1/k^2, n+1,inf)
Why does that help us? It does so because we can write each of those infinite sums in a known form. As we saw above, the first of those infinite series is just the Eulerian sum, which reduces to pi^2/6, as I said above.
The second term is only simply written using the polygamma function.
Can I explain Polygamma very easily? Hmm. I think the best answer is to just do some reading about the polygamma function, perhaps starting with that page. One thing you might consider is the gamma function itself is a very pretty thing, with the property that
gamma(n+1) == n*gamma(n)
In fact, we often use the gamma function as a factorial replacement, because the factorial function itself has a very similar property, where factorial(n+1) = (n+1)*factorial(n). Then we can take one more step, and find that for integer values of n, factorial(n)==gamma(n+1). There are lots of very pretty things you can do now, but I digress.
The whole point of my digression is to convince you to look at the gamma function, and the derivatives of the log of the gamma function. Again, do some reading, as this is wonderful a rabbit hole wherein I would be forced to write an entire book. And even though that book would be a wonderful read, it has already been written too many times by others more capable than am I.
The whole point of this is to convince that you can indeed write that second infinite sum using psi. For those that need some convincing, we can do a simple enough test. A wonderful feature about tools like MATLAB is they allow you to explore a result to convince yourself if it makes sense. Thus is psi(n+1,1) really the same as the infinite sum: symsum(1/k^2,n+1,inf)?
psi(1,4)
ans =
0.283822955737115
>> sum(1./(4:100).^2)
ans =
0.273872789073782
>> sum(1./(4:1000).^2)
ans =
0.282823455570449
>> sum(1./(4:10000).^2)
ans =
0.283722960736952
>> sum(1./(4:100000).^2)
ans =
0.283812955787118
>> sum(1./(4:1000000).^2)
ans =
0.283821955737617
>> sum(1./(4:10000000).^2)
ans =
0.283822855737108
>> psi(1,4)
ans =
0.283822955737115
As you should see, as I increase the series to have more and more terms, we really do seem to approach the same limit. Thus, with roughly 1e7 terms, we have convergence to the 7th decimal place, compared to psi(1,4). Unfortunately, the series does not converge that quickly, so I would need to take a huge number of terms to get double precision convergence.
All of this begs the question however, is is possible for a simple form for psi(1,n+1) to ever be written? And that we can claim does have an emphatic answer: No. Why not? It reduces to the transcendental nature of pi. That is, we know that for ANY finite value of n, the sum:
symsum(1/k^2,1,n)
does exist in a rational form. That is, for several values of n, we can write the sums, then expanding it in a purely numeric form. For a finite sum, a symsum is just a sum of syms. (I hope you can say that three times fast. Much like "How much wood can a woodchuck chuck...")
Again, lets do a few examples. (Sorry, but I never know just how simply I need to explain something, because for an online answer, I must reach for the lowest common denominator of those individuals who may one day read this.)
S = @(k) sum(1./sym(1:k).^2);
S(1)
ans =
1
>> S(2)
ans =
5/4
>> S(3)
ans =
49/36
>> S(4)
ans =
205/144
>> S(5)
ans =
5269/3600
And those are the exactly correct rational numbers, summed using exact arithmetic. They are not approximate floating point results. But you need to recognize that any finite sum of integer fractions is in the end a simple rational fraction, therefore a purely rational number.
We have found this relationship:
symsum(1/k^2, 1, n) == pi^2/6 - psi(1,n + 1)
it follows that
pi^2/6 == psi(1,n + 1) + symsum(1/k^2, 1, n)
and then we would have:
pi == 6*sqrt(psi(1,n + 1) + symsum(1/k^2, 1, n))
Essentially, we have just written a simple formula for pi. And since the finite part of that sum is clearly a rational number, this also leads us to the direct conclusion that psi(1,n+1) must itself be trascendental, since pi is itself transcendental. Effectively, it leads us to the simple conclusion that we won't ever be able to expand psi(1,n+1) in any nice algebraic form, because if we could, then we would also have a simple algebraic formula to compute pi.
What does this tell us? While we know that the desired finite sum, if evaluated for any given known integer n is representable as a rational number, it does look like we are now at a dead end. The simplest way to write the desired finite sum as a function of a general integer n is exactly what symsum yields in terms of pi and the function psi.

  0 Comments

Sign in to comment.

Sign in to answer this question.