Series, sums, products

Series

Series is a number defined by

S=(lim_(N->∞))(Subscript[S, N])

where (Subscript[S, N]) is a sum

(Subscript[S, N])=(∑_(n=0)^N)(Subscript[a, n])

and (Subscript[a, n]) are series terms. Thus if the series is converging for N->∞, it is given by

S=(∑_(n=0)^∞)(Subscript[a, n])

Power series considered earlier is a generalization of the series that are functions of their arguments. Another generalization of the series are series with summation over several indices. Example of series in physics are Coulomb or dipolar interaction of atoms arranged in a crystal lattice. Also series may arise in calculating integrals and solving differential equations.

In some cases series can be calculated analytically. An example are Taylor series, e.g.

(∑_(n=0)^∞)(1/(n!))=e

that is a particular case of the general result

(∑_(n=0)^∞)((x^n)/(n!))=(e^x)

or

(∑_(n=0)^∞)(1/(2^n))=(1/(1-(1/2)))=2

that is a particular case of geometrical progression

Mathematica can calculate series

(∑_(n=0)^∞)(1/(2^n))

(∑_(n=0)^∞)((1-n)/(2^n))
(∑_(n=0)^∞)((1+n)/(2^n))

(∑_(n=0)^∞)((1+n)/(2^(a n)))

(∑_(n=0)^∞)((Sin[n])/(2^n))

(∑_(n=1)^∞)(1/(n^2))

(∑_(n=1)^∞)(1/(n^4))

If Mathematica cannot do it analytically, it can do it numerically

N[(∑_(n=1)^∞)((Log[n]Sin[n])/(2^n))]

Some series are expressed through special functions

(∑_(n=1)^∞)(1/(n^3))

and, more generally,

(∑_(n=1)^∞)(1/(n^x))

- Riemann zeta-function.

Plot[Zeta[x],{x,-3,7}]

Sums

Sums also can be calculated by Mathematica

(∑_(n=0)^m)(1/(2^n))
Simplify[%]
(∑_(n=0)^m)(1/(2^n))//Simplify

(∑_(n=0)^m)(n^3)

In general, the analytical expression for a series is simpler than that for the sum

(∑_(n=0)^∞)(1/(n!))
(∑_(n=0)^m)(1/(n!))

Products

Products

P=(∏_(n=0)^N)(Subscript[b, n])

can be reduced to sums by taking their logarithm:

Log[P]=Log[(∏_(n=0)^N)(Subscript[b, n])]=(∑_(n=0)^N)Log[(Subscript[b, n])],

thus

P=Exp[(∑_(n=0)^N)Log[(Subscript[b, n])]]

Mathematica can calculate products

(∏_(n=1)^N)n

(∏_(n=1)^N)Exp[n]

etc.

Analytical calculation of sums and series

Different methods can be used

Differentiation over parameter

To calculate

(∑_(n=0)^∞)(n/(2^n))

consider, for |x|<1,

f[x]=(∑_(n=0)^∞)(x^n)=(1/(1-x))

One has

f'[x]=(∑_(n=0)^∞)(nx^(n-1))=(1/x)(∑_(n=0)^∞)(nx^n)=(1/(((1-x))^2))

From here follows

g[x]=(∑_(n=0)^∞)(nx^n)=(x/(((1-x))^2))

Now

(∑_(n=0)^∞)(n/(2^n))=g[1/2]=2.

Series with higher powers of n can be calculated via high-order derivatives of the generating function f[x].

Using discrete derivatives

The sum of the terms having the form

(Subscript[a, n])=(Subscript[b, n+1])-(Subscript[b, n])

can be easily calculated:

(∑_(n=0)^N)(Subscript[a, n])=(∑_(n=0)^N)((Subscript[b, n+1])-(Subscript[b, n]))=(∑_(n=0)^N)(Subscript[b, n+1])-(∑_(n=0)^N)(Subscript[b, n])=(∑_(n=1)^(N+1))(Subscript[b, n])-(∑_(n=0)^N)(Subscript[b, n])=(Subscript[b, N+1])-(Subscript[b, 0])

Here (Subscript[a, n]) can be interpreted as a discrete derivative of (Subscript[b, n]) and the sum is similar to an integral of this derivative that is given by the difference of (Subscript[b, n]) between the upper and lower summation limits. One of the methods of summing the arithmetic progression

(∑_(n=0)^N)n

uses this idea. Using the discrete derivative

(((n+1))^2)-(n^2)=2n+1

one can write

(((N+1))^2)-0=(∑_(n=0)^N)((((n+1))^2)-(n^2))=(∑_(n=0)^N)(2n+1)=2(∑_(n=0)^N)n+(∑_(n=0)^N)1=2(∑_(n=0)^N)n+N+1,

wherefrom one obtains

(∑_(n=0)^N)n=(1/2)((((N+1))^2)-N-1)=(1/2)((N^2)+N)=((N(N+1))/2)

For the arithmetic progression this method if not the fastest. However, this method works for sums with higher powers. For instance, one can calculate

(∑_(n=0)^N)(n^2)

using the discrete derivative

(((n+1))^3)-(n^3)=3(n^2)+3n+1

and the formula for (∑_(n=0)^N)n. One writes

(((N+1))^3)-0=(∑_(n=0)^N)((((n+1))^3)-(n^3))=(∑_(n=0)^N)(3(n^2)+3n+1)=3(∑_(n=0)^N)(n^2)+3(∑_(n=0)^N)n+(∑_(n=0)^N)1=3(∑_(n=0)^N)(n^2)+((3N(N+1))/2)+N+1,

wherefrom one obtains

(∑_(n=0)^N)(n^2)=(1/3)((((N+1))^3)-((3N(N+1))/2)-N-1)=(1/3)((N^3)+3(N^2)+3N+1-((3N(N+1))/2)-N-1)=(N/3)((N^2)+3N+3-((3(N+1))/2)-1)=(N/3)((N^2)+((3N)/2)+(1/2))=(N/6)(2(N^2)+3N+1)=(N/6)(N+1)(2N+1)

Sums of any powers of n can be calculated with this method but calculations become more and more difficult.

Convergence of series

The main question about the series is their convergence, i.e., the existence of the limit S=(lim_(N->∞))(Subscript[S, N]). There are several criteria of convergence.

Ratio criterium

A useful criterion in the so-called ratio criterion. With the asymptotic ratio defined as

ρ=(lim_(n->∞))((Subscript[a, n+1])/(Subscript[a, n]))

the criterium reads

({{ρ<1,, series converges}, {ρ=1,, series can converge or diverge}, {ρ>1,, series diverges}})

For instance, for the series

(∑_(n=0)^∞)((n!)/((2n)!))

one has

ρ=(lim_(n->∞))(((n+1)!)/((2(n+1))!))(((2n)!)/(n!))=(lim_(n->∞))((n+1)/((2n+2)(2n+1)))=0

and the series converges. Mathematica can express this series via the error function:

(∑_(n=0)^∞)((n!)/((2n)!))
N[%]

For the series representing the zeta function

(∑_(n=1)^∞)(1/(n^x))

one has

ρ=(lim_(n->∞))((n^x)/(((n+1))^x))=(lim_(n->∞))(1/(((1+(1/n)))^x))=(1/((lim_(n->∞))(((1+(1/n)))^x)))=1

and one cannot decide whether this series converges or diverges.

Integral criterium

Even more useful is the integral criterion that says that the series converges or diverges if the similar integral converges or diverges at the upper limit. This criterion is convenient because it is easier to analytically calculate integrals than series. In particular, for the series representing the zeta function above one obtains that the series converges for x>1.

Alternating series

Alternating series are series of the type

S=(∑_(n=0)^∞)(((-1))^n)(Subscript[a, n]),    (Subscript[a, n])>0.

Alternating series converge if (Subscript[a, n]) decreases steadily to zero in the limit n->∞. Example:

(∑_(n=1)^∞)((((-1))^n)/n)
N[%]

(∑_(n=1)^∞)((((-1))^n)/(Sqrt[n]))
N[%]

(∑_(n=1)^∞)((((-1))^n)/(n^(1/10)))
N[%]

Coulomb energy of an ionic lattice

General

Coulomb energy of a system of charges is given by

U=(k/2)(∑_(i,j)^N)(((Subscript[q, i])(Subscript[q, j]))/(|(Subscript[r, i])-(Subscript[r, j])|)),

where k is related to the system of units and the factor 1/2 compensates double counting of terms. The Coulomb energy corresponding to charge (Subscript[q, i]) is

(Subscript[U, i])=(k/2)(∑_j^N)(((Subscript[q, i])(Subscript[q, j]))/(|(Subscript[r, i])-(Subscript[r, j])|)),

so that

U=(∑_j^N)(Subscript[U, i]).

Consider alternating ±q charges on an infinite square lattice with lattice constant a. The energies corresponding to each charge in this infinite system are the same, so that one can consider the energy of the charge at the origin of the coordinate system, (Subscript[r, i])= 0. Using

(Subscript[r, j])=(Subscript[e, x])a (Subscript[n, x])+(Subscript[e, y])a (Subscript[n, y]),       |(Subscript[r, j])|=a(Sqrt[(Subscript[n, x]^2)+(Subscript[n, y]^2)]),

one obtains

(Subscript[U, 0])=((k (q^2))/(2a))M,

where

M=(Underoverscript[∑', ((Subscript[n, x]),(Subscript[n, y])=-∞), ∞])((((-1))^((Subscript[n, x])+(Subscript[n, y])))/(Sqrt[(Subscript[n, x]^2)+(Subscript[n, y]^2)]))

is the Madelung constant. Prime at the sum means that the non-existent term with (Subscript[n, x])=(Subscript[n, y])= 0 is omitted. The sum can be calculated numerically by at first summing in the range -N ≤ (Subscript[n, x]),(Subscript[n, y]) ≤ N and then taking the limit N->∞. Similar formulas can be written in 1d and 3d cases. In all cases M < 0, because alternating arrangements of positive and negative ions reduce the energy.

One-dimensional toy model

Let us at first consider the 1d model

FM1[n_]:=If[n==0,0,((((-1))^n)/(Abs[n]))]
M1N[NN_]:=(∑_(n=1)^NN)2.((((-1))^n)/n) 
M1=(∑_(n=1)^∞)2((((-1))^n)/n)    (* Exact result *)
N[M1]

M1NList=Table[{NN,M1N[NN]},{NN,10,400,5}]  (* Approximations with different N *)

Plot1Exact=Plot[M1,{NN,10,400},PlotStyle->{Black,Thick}];
Plot1Approx=ListPlot[M1NList,PlotRange->All];
Show[Plot1Approx,Plot1Exact]

Convergence to the exact result is slow. To improve the convergence, take the contribution of the last point n=N with the factor 1/2. This is reasonable because the last point in fact should be shared 50%-50% between the interval taken into account and the remaining region. So we define

M1NCorrected[NN_]:=(∑_(n=1)^NN)((((-1))^n)/n) If[n==NN,1,2.]

Now do everything again and compare

M1NList=Table[{NN,M1N[NN]},{NN,5,100,5}]
M1NCorrectedList=Table[{NN,M1NCorrected[NN]},{NN,5,100,5}]  (* Approximations with different N *)

Plot1Exact=Plot[M1,{NN,5,100},PlotStyle->{Black,Thick}];
Plot1Approx=ListPlot[M1NList,PlotStyle->{Blue,PointSize[0.015]},PlotRange->All];
Plot1ApproxCorr=ListPlot[M1NCorrectedList,PlotRange->All,PlotStyle->{Red,PointSize[0.02]},PlotRange->All];
Show[Plot1Approx,Plot1ApproxCorr,Plot1Exact]

One can see that the corrected procedure has much better convergence. Note that the total charge in the summation range is zero, that is, the electroneutrality holds:

Q1NCorrected[NN_]:=(∑_(n=0)^NN)(((-1))^n) If[n==0||n==NN,1.,2.]

Q1NCorrected[20]
Q1NCorrected[21]

Two-dimensional model

Here there is no analytical expression for the Madelung constant. Similarly to the 1d model we define

FM2[nx_,ny_]:=If[nx==ny==0,0,((((-1.))^(nx+ny))/(Sqrt[(nx^2)+(ny^2)]))]
M2N[NN_]:=(∑_(nx=-NN)^NN)(∑_(ny=-NN)^NN)FM2[nx,ny]
M2NCorrected[NN_]:=(∑_(nx=-NN)^NN)(∑_(ny=-NN)^NN)FM2[nx,ny]Which[Abs[nx]==Abs[ny]==NN,(1/4),(Abs[nx]==NN)&&(Abs[ny]≠NN)||(Abs[ny]==NN)&&(Abs[nx]≠NN),(1/2),True,1]

M2NList=Table[{NN,M2N[NN]},{NN,2,50}]
M2NCorrectedList=Table[{NN,M2NCorrected[NN]},{NN,1,50}]

Plot2Approx=ListPlot[M2NList,PlotStyle->{Blue,PointSize[0.015]},PlotRange->All];
Plot2ApproxCorr=ListPlot[M2NCorrectedList,PlotRange->All,PlotStyle->{Red,PointSize[0.02]},PlotRange->All];
Show[Plot2Approx,Plot2ApproxCorr]

One can see that in 2d convergence of the Madelung constant to the exact result in 2d (M = -1.61554) is even much faster than in 1d.