# It starts with an innocent power series

In this post I play around with some blogdown features. For that a crazy proof is resurrected. The sum of a funny, innocent, yet infinite power series is derived. It is very likely, that someone else did this before. Anyway, it is sufficient to brag around with some math and to include an R ggplot.

$\displaystyle\sum_{n=0}^{\infty}\frac{n^a}{a^n}$

I just had wondered if there was a more handsome form of its sum value. This thing is a little bit different from an ordinary power series. The variable $x$ normally does not like to live in the exponent of the coefficient, except in this sweet expression..

The equation we will get: $\sum_{n=0}^{\infty}\frac{n^a}{a^n}=\sum_{k=1}^{a}k!\cdot S_2(a,k)\cdot\frac{a}{(a-1)^{k+1}}$ To get there, a Lemma is required.

Lemma 1: For the $n$-times differentiable function $A(x)$, $n>0$, and the differential operator ${\mathrm{D}}=\frac{\mathrm{d}}{\mathrm{d}x}$ there is the equation: $(x{\mathrm{D}})^n A(x) = \sum_{k=1}^n S_2(n,k)\cdot x^k\cdot {\mathrm{D}}^k A(x)$ where $S_2(n,k)$ is the Stirling Number of the second kind. The Stirling Number of the second kind can be recursively defined as:

$S_2(n,k) = S_2(n-1,k-1) + k\cdot S_2(n-1,k)$

while $S_2(n,0)=0$ for $n>0$ and $S_2(n,m)=0$ for $m>n$.

Proof: We prove the claim by induction on $n$. For $n=1$ it is easy to see that

$x{\mathrm{D}}A(x) = \sum_{k=1}^{1}S_2(1,k)\cdot x^k {\mathrm{D}}^k A(x)\,,$ where $S_2(1,1)=1$. This starts the induction. The inductive hypothesis for $n-1$ states:

$(x{\mathrm{D}})^{n-1} A(x) = \sum_{k=1}^{n-1} S_2(n-1,k)\cdot x^k\cdot {\mathrm{D}}^k A(x)$

Hence it holds: \begin{aligned} (x{\mathrm{D}})^n A(x) &= x{\mathrm{D}}\left( (x{\mathrm{D}})^{n-1} A(x) \right) \\ &= x{\mathrm{D}}\left( \sum_{k=1}^{n-1} S_2(n-1,k)\cdot x^k\cdot {\mathrm{D}}^k A(x) \right)\\ &=\sum_{k=1}^{n-1}k\cdot S_2(n-1,k)x^k{\mathrm{D}}^kA(x) + \sum_{k=1}^{n-1}S(n-1,k)x^{k+1}{\mathrm{D}}^{k+1}A(x)\\ &=\sum_{k=1}^{n-1}k\cdot S_2(n-1,k)x^k{\mathrm{D}}^kA(x) + \sum_{k=2}^{n}S(n-1,k-1)x^{k}{\mathrm{D}}^{k}A(x)\\\end{aligned}

The indices can be extended as follows, because that Stirling Numbers are zeros: \begin{aligned} (x{\mathrm{D}})^n A(x) &= \sum_{k=1}^{n}k\cdot S_2(n-1,k)x^k{\mathrm{D}}^kA(x) + \sum_{k=1}^{n}S(n-1,k-1)x^{k}{\mathrm{D}}^{k}A(x)\\ &=\sum_{k=1}^{n}\left( k\cdot S_2(n-1,k)+S_2(n-1,k-1) \right) x^k {\mathrm{D}}^k A(x)\\ &=\sum_{k=1}^{n} S_2(n,k) x^k {\mathrm{D}}^k A(x)\end{aligned}

Theorem: For a given $a>1$ the sum of the following power series can be calculated: $\sum_{n=0}^{\infty}\frac{n^a}{a^n}=\sum_{n=0}^{\infty}n^a\left(\frac{1}{a}\right)^n=\sum_{k=1}^{a}k!\cdot S_2(a,k)\cdot\frac{a}{(a-1)^{k+1}}$

Proof: As power series are differentiable on the interior of the domain of convergence one can write, if $|x|<1$: \begin{aligned} \sum_{n=0}^{\infty}n^ax^n &= (x{\mathrm{D}})^a \sum_{n=0}^{\infty}x^n\\ &= (x{\mathrm{D}})^a \frac{1}{1-x}\end{aligned}

By Lemma [1] the result can be computed as: \begin{aligned} (x{\mathrm{D}})^a \frac{1}{1-x} &= \sum_{k=1}^a S_2(a,k)\cdot x^k\cdot {\mathrm{D}}^k \frac{1}{1-x}\\ &=\sum_{k=1}^a S_2(a,k)\cdot k!\cdot \frac{x^k}{(1-x)^{k+1}}\end{aligned} With $x=\frac{1}{a}$ the result can be written as follows: \begin{aligned} \sum_{n=0}^{\infty}\frac{n^a}{a^n} &= \sum_{k=1}^a k!\cdot S_2(a,k)\cdot \frac{a^{-k}}{\left(1-\frac{1}{a}\right)^{k+1}}\\ &=\sum_{k=1}^{a}k!\cdot S_2(a,k)\cdot\frac{a}{(a-1)^{k+1}}\end{aligned}

Example: For $a=2$ the sum is $\sum_{n=0}^{\infty}\frac{n^2}{2^n}=1!\frac{4}{2}+2!\frac{8}{4}=6$. With $a=3$ you will obtain $\frac{66}{16}$.

In Figure 1 you can see the behaviour of the sum approximation.

library(ggplot2)
library(plyr)
library(dplyr)
library(scales)
library(gmp)
sumit <- function(len,a) {
tsum <- 0
for(n in seq(0,len-1)) {
tsum <- tsum + n**a / a**n
}
return(tsum)
}
sumdir <- function(a) {
tsum <- 0
for(ka in seq(1,a)) {
tsum <- tsum + factorial(ka) * Stirling2(a,ka,method="direct")*a/(a-1)**(ka+1)
}
return(as.numeric(tsum))
}
sumit_to_frame <- function(lengths,a) {
dfp <- data.frame(x=lengths, y=sapply(lengths,sumit,a=a), a=as.integer(a))
dfp$c <- as.numeric(sumdir(a)) return(dfp) } lens <- seq(1,25) dfp <- sumit_to_frame(lens, 2) %>% bind_rows( sumit_to_frame(lens, 3) ) %>% bind_rows( sumit_to_frame(lens, 4) ) %>% bind_rows( sumit_to_frame(lens, 5) ) dfp$a <- as.factor(dfp\$a)

my_theme <-  theme_bw() + theme(axis.title.x = element_text(size=18),
axis.title.y = element_text(size=18),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
plot.margin = unit(c(8,10,1,1), "pt") # required otherwise labels are clipped in pdf output
)
my_theme <- my_theme + theme(#legend.title = element_blank(),
legend.title = element_text(size=16, face="bold"),
legend.text = element_text( size = 16),
legend.position="right",
legend.direction="vertical",
legend.box ="horizontal",
legend.box.just ="right",
legend.background = element_rect(colour = 'white', fill = 'white', size = 0., linetype='dashed'),
legend.key = element_rect(colour = 'white', fill = 'white', size = 0., linetype='dashed'),
legend.key.width = unit(1.1, "cm"),
panel.grid.major = element_line(color="#444444", size=0.5, linetype=3),
panel.grid.minor = element_line(color="#555555", size=0.5, linetype=3)
)

gg <- ggplot(dfp, aes(x,colour=a)) + geom_line(aes(y=y),size=1.2) + geom_line(aes(y=c),linetype="dashed",size=1.0) + geom_point(aes(y=y),size=2.5)
gg <- gg + my_theme
gg <- gg + xlab("iterations N")
gg <- gg + ylab(expression(sum(frac(n**a,a**n),n==0,N)))
gg
Share