Discrete Calculus with Mathematica

({{, , , , }})

Devendra Kapadia and Charles Pooh

Wolfram Research, Inc.

Outline of Presentation

Overview of Discrete Calculus in Mathematica 7

Evaluation of symbolic sums and solutions of difference equations

Sequence-guessing and generating functions

Applications to chaotic dynamical systems and heat dissipation

Architecture of Sum and user-extensibility

Questions and Discussion

Overview of Discrete Calculus

What is Discrete Calculus?

Discrete Calculus is the study of sequences using finite differences, sums, products, and generating functions.

Practical problems typically require the solution of difference equations given some initial or boundary conditions.

Many areas such as Probability Theory, Dynamical Systems and the Theory of Algorithms use methods from Discrete Calculus.

Examples of Sequences

The Fibonacci numbers are defined by a recurrence relation.

Table[(Fibonacci[n]),{n,1,10}]

The harmonic numbers are defined by (∑_(k=1)^n)(1/k).

Table[(HarmonicNumber[n]),{n,1,10}]

DiscretePlot[(HarmonicNumber[n]),{n,1,30},PlotStyle->{PointSize[0.015],Red},Background->Yellow]

Sequences such as Factorial[n] and BarnesG[n] grow very quickly.

Table[n!,{n,0,10}]

Table[(BarnesG[n]),{n,1,10}]

Problems involving periodic sequences can be handled very efficiently once the period has been determined.

GraphicsRow[{DiscretePlot[sin(((π n)/6)),{n,30},PlotStyle->{PointSize[0.015],Red},Background->Yellow],
DiscretePlot[(Mod[n, 4]),{n,20},PlotStyle->{PointSize[0.015],Red},Background->Yellow]},ImageSize->700]

Hypergeometric functions often have simple expressions for certain discrete ranges of their parameters.

(Hypergeometric2F1[a, b, c, x])/.{a->-2}

Evaluation of Sums and Products

The sum of the squares of the first billion natural numbers can be computed in two ways:

(∑_(k=1)^(10^9))(k^2)

(∑_(k=1)^n)(k^2)

%/.{n->(10^9)}

Infinite sums played an important role in the development of quantum theory and its applications, such as proving the Stefan-Boltzmann law.

Assuming[(h/(k T))>0,(∫_0^∞)(((2 π (ν^2)) (∑_(n=0)^∞)h n ν (E^(-((h n ν)/(k T)))))/((c^3) (∑_(n=0)^∞)(E^(-((h n ν)/(k T))))))\[DifferentialD]ν]

The following multiple sum has a simple answer due to  a symmetry in the denominator:

(∑_m^∞)(∑_n^∞)(((m^2) n)/((2^m) (m (2^n)+(2^m) n)))

Hypergeometric term products such as the following can be evaluated in terms of Gamma and BarnesG:

(∏_(k=1)^n)(((7^k) (2k+3))/((1+k)!))

Factor[(DiscreteRatio2[(((7^k) (2 k+3))/((k+1)!)), k])]

(∏_(k=1)^n)(((7^k) (2 k+3))/((k+1)!))

Solutions of Difference Equations

This difference equation generates the Fibonacci sequence.

RSolve[{a(n)==a(n-1)+a(n-2),a(1)==1,a(2)==1},a,n]

Table[a(n)/.%[[1]],{n,1,10}]

The number of moves required to solve the Towers of Hanoi puzzle with n disks and three rods is given by:

RSolve[{T(n)==2 T(n-1)+1,T(1)==1},T,n]

We can use RecurrenceTable to generate the moves for different values of  n:

RecurrenceTable[{T(n)==2 T(n-1)+1,T(1)==1},T,{n,10}]

Sequence-guessing and sequence-generation

In the following, we find a closed form corresponding to a given sequence:

FindSequenceFunction[{1,3,4,7,11,18,29,47},n]

For this example, the generating function is rational in x:

FindGeneratingFunction[{1,3,4,7,11,18,29,47},x]

Series[%,{x,0,7}]

Symbolic Summation in Mathematica

Indefinite Sums

Consider the following indefinite sum :

(∑_k)((k^2)+3 k+5)

Then the difference of this indefinite sum is the original expression:

(DifferenceDelta2[%, k])

The indefinite sum of a rational function is composed of a rational part and a transcedental part based on PolyGamma :

(∑_k)(((k^4)+3 (k^3)-(k^2)+10)/((k^5)+11 (k^4)+35 (k^3)+55 (k^2)+45 k+21))

Simplify[(DifferenceDelta2[%, k])-(((k^4)+3 (k^3)-(k^2)+10)/((k^5)+11 (k^4)+35 (k^3)+55 (k^2)+45 k+21))]

Indefinite sums of hypergeometric terms can often be computed by Gosper' s algorithm :

(∑_k)(((27^(-k)) (3 k)!)/(k! (k+1)! (k+2)!))

Homotopy methods are used for indefinite sums of exact expressions :

(∑_k)(-((u(k))^13)-((u(k))^2)-u(k)+((u(k+1))^13)+((u(k+1))^2)+u(k+1))

(DifferenceDelta2[%, k])

Finite Sums

Summation by parts is a very useful technique for problems involving suitable pairs of sequences :

∑ f[k] g[k] δk==g[k] (Subscript[Δ, k]^(-1))f[k]-∑(Subscript[Δ, k])g[k] (Subscript[E, k])(Subscript[Δ, k]^(-1))f[k] δk

We use this idea to compute sums involving special functions, such as the following :

(∑_(k=1)^m)((PolyGamma[k])/k)

Zeilberger's algorithm provides an important method for definite hypergeometric sums:

(∑_(k=0)^n)((k! (Binomial[n, k]))/((k+n+1)!))

The key idea in this case is to solve a recurrence for the sum viewed as a function of n :

Sum`SumZeilbergerRecurrence[(k!*Binomial[n, k])/(k + n + 1)!, {k, 0, n}, y]

Infinite Sums

The first step in evaluating infinite series is to test whether they are convergent, using SumConvergence :

SumConvergence[((x^n)/n),n]

This information is used to generate conditions in the result from Sum :

Sum[((x^n)/n),{n,1,∞},GenerateConditions->True]

The following quadratic Euler sum can be evaluated in terms of Zeta:

(∑_(k=1)^∞)(((((HarmonicNumber[k])))^2)/(k^7))

Sums involving Fibonacci and Stirling numbers are evaluated by converting them to rational - exponential or hypergeometric series :

(∑_(k=2)^∞)(((Fibonacci[2 k+7]) (StirlingS1[k+1, k-1]))/((4^k) (((k!))^2)))

The following remarkable infinite sum is due to Srinivas Ramanujan, and is done by table look-up.

(((Sqrt[8]) (∑_(n=0)^∞)(((26390 n+1103) (4 n)!)/((396^(4 n)) (((n!))^4))))/9801)

Divergent Series

Here is a divergent alternating series :

(∑_(n=1)^∞)(((-1))^n)

We can use Abel regularization to assign a value to this series:

(∑_(n=0)^∞)(((-1))^n) (x^n)

(lim_(x->1)) %

Sum[(((-1))^n),{n,0,∞},Regularization->"Abel"]

The following finite sum is computed using Borel regularization :

(∑_(k=0)^n)k! (((-1))^(n-k)) (Binomial[n, k])

Sum[k! (((-1))^(n-k)) (Binomial[n, k]),{k,0,∞},Regularization->"Borel"]

q-Summation

Indefinite sums of q - rational functions may involve the QPolyGamma function :

expr=(((q^(2 k-1))-(q^(5 (2 k-1))))/(1-(q^(6 (2 k-1)))));

(∑_k)expr

(DifferenceDelta2[%, k])

Together[%-expr]

We can also find q-rational solutions of with q-rational coefficients:

eqn=(-15 (q^(k+1))-17 (q^(k+2))+(q^(2 k+2))-8 (q^(k+3))-11 (q^(2 k+3))-45 (q^(k+4))-30 (q^(2 k+4))-45 (q^(k+5))-28 (q^(2 k+5))-2 (q^(3 k+5))-9 (q^(k+6))-40 (q^(2 k+6))-20 (q^(3 k+6))-51 (q^(2 k+7))-21 (q^(3 k+7))-39 (q^(2 k+8))-18 (q^(3 k+8))-4 (q^(4 k+8))-18 (q^(2 k+9))-36 (q^(3 k+9))-7 (q^(4 k+9))-27 (q^(3 k+10))-6 (q^(4 k+10))-15 (q^(3 k+11))-8 (q^(4 k+11))-(q^(5 k+11))-9 (q^(3 k+12))-8 (q^(4 k+12))-(q^(5 k+12))-6 (q^(4 k+13))-3 (q^(4 k+14))-(q^(5 k+14))-(q^(5 k+15))+3 (q^k)-18 (q^2)-15 q+3) y(k+1)+(18 (q^(k+1))+5 (q^(k+3))+8 (q^(2 k+3))+30 (q^(k+4))+12 (q^(2 k+4))+18 (q^(k+5))+21 (q^(2 k+5))+3 (q^(3 k+5))+11 (q^(2 k+6))+3 (q^(3 k+6))+23 (q^(2 k+7))+11 (q^(3 k+7))+33 (q^(2 k+8))+19 (q^(3 k+8))+(q^(4 k+8))+6 (q^(3 k+9))+3 (q^(4 k+9))+14 (q^(3 k+10))+3 (q^(4 k+10))+18 (q^(3 k+11))+4 (q^(4 k+11))+4 (q^(4 k+12))+(q^(5 k+12))+3 (q^(4 k+13))+3 (q^(4 k+14))+(q^(5 k+15))-3 (q^k)+18 q-3) y(k+2)+(-3 (q^(k+1))+17 (q^(k+2))-(q^(2 k+2))+3 (q^(k+3))+3 (q^(2 k+3))+15 (q^(k+4))+17 (q^(2 k+4))-(q^(3 k+4))+27 (q^(k+5))+14 (q^(2 k+5))+6 (q^(3 k+5))+9 (q^(k+6))+15 (q^(2 k+6))+3 (q^(3 k+6))+30 (q^(2 k+7))+11 (q^(3 k+7))+21 (q^(2 k+8))+21 (q^(3 k+8))+3 (q^(4 k+8))+9 (q^(2 k+9))+7 (q^(3 k+9))+4 (q^(4 k+9))+15 (q^(3 k+10))+3 (q^(4 k+10))+12 (q^(3 k+11))+4 (q^(4 k+11))+(q^(5 k+11))+4 (q^(4 k+12))+3 (q^(4 k+13))+(q^(5 k+14))+18 (q^2)-3 q) y(k);

sol=RSolve[eqn==0,y,k]

Simplify[eqn/.sol[[1]]]

Sequence guessing

1. What is the next element in the sequence 3, 7, 13, 21, 31, 43, 57?

2. Find a formula that generates 2, 5/2, 10/3, 17/4, 26/5, 37/6, 50/7?

3. Find a generating function having first coefficients 1, 1, 2, 2, 3, 4, 5, 6, 7, 8, 11, 12?

FindSequenceFunction

FindSequenceFunction[{(Subscript[a, 1,])(Subscript[a, 2]),(Subscript[a, 3]),...},n] attempts to find a simple function that yields the sequence (Subscript[a, n]).

FindSequenceFunction[{1,1,2,3,5,8,13},n]

Table[Fibonacci[n],{n,10}]

Functions

Periodic sequence

Polynomial

Rational function

Hypergeometric terms

Recurrence equations

FindSequenceFunction[{2,3/2,7/3,13/4,26/5,49/6,92/7,169/8,307/9,551/10,980/11,1729/12,3030/13,5279/14, 
 9151/15,15793/16,27150/17,46513/18},n]

Generalizations and Extensions

FindSequenceFunction[{1,2,3,5,17,305,34865,24918065,125436246065,5056710181206065},n]

Options

FunctionSpace

ValidationLength

FindGeneratingFunction

FindGeneratingFunction[{(Subscript[a, 1]),(Subscript[a, 2]),...},x] attempts to find a simple generating function in x whose (n^th) series coefficient is (Subscript[a, n]).

FindGeneratingFunction[{1,1,2,3,5,8,13},x]

Series[%, {x,0,7}]

Functions

Rational function

FindGeneratingFunction[{1,1,2,2,3,4,5,6,7,8,11,12},x]

Hypergeometric function

FindGeneratingFunction[{1,3/2,9/5,2,15/7,9/4,7/3,12/5},x]

Differential equations

Generalizations and Extensions

FindGeneratingFunction[{1,-t,1/2+1/2t^2,-t/2-t^3/6,1/8+t^2/4+t^4/24,-t/8-t^3/12-t^5/120},x]

FindGeneratingFunction[{0,1,3,12,55,273,1428,7752,43263,246675},x]

Options

FunctionSpace

ValidationLength

LinearRecurrence & FindLinearRecurrence

LinearRecurrence

LinearRecurrence[ker,init,n] gives the sequence of length n obtained by iterating the linear recurrence with kernel ker starting with initial values init.

LinearRecurrence[{1,1},{1,1},10]

LinearRecurrence[{a,b},{c,d},5]

LinearRecurrence[{1,1},{{1,2},{2,1}},10]

FindLinearRecurrence

FindLinearRecurrence[list] finds if possible the minimal linear recurrence that generates list.

FindLinearRecurrence[{0,1,2,12,45,200,832,3549,14994,63580}]

FindLinearRecurrence gives the shortest kernel ker which reproduce list from 
LinearRecurrence[ker, init, ...].

LinearRecurrence[{3,6,-3,-1},{0,1,2,12},15]

Sequence guessing

1. What is the next element in the sequence 3, 7, 13, 21, 31, 43, 57?

FindSequenceFunction[{3,7,13,21,31,43,57},8]

2. Find a formula that generates 2, 5/2, 10/3, 17/4, 26/5, 37/6, 50/7?

FindSequenceFunction[{2,5/2,10/3,17/4,26/5,37/6,50/7},n]

3. Find a generating function with first coefficients 1, 1, 2, 2, 3, 4, 5, 6, 7, 8, 11, 12?

FindGeneratingFunction[{1,1,2,2,3,4,5,6,7,8,11,12},x]

Iterative Solutions of Difference Equations

Scope of RecurrenceTable

RecurenceTable can find iterative solutions for initial value problems involving:

Ordinary difference equations (linear, nonlinear, single, system)

Evolutionary partial  difference equations  with two variables.

Difference-algebraic equations (linear constant coefficient systems)

The Logistic Map

The onset of chaos in the logistic map may be modelled as follows :

GraphicsGrid[Partition[Table[ListPlot[RecurrenceTable[y(n+1)==r y(n) (1-y(n))∧y(0)==(1/3),y,{n,1,600}],ImageSize->150,PlotStyle->Red,Background->Yellow],{r,3,3.8,0.04}],5]]

A Rabbit Fractal

Douady's rabbit fractal may be constructed in the following steps.

Choose a complex number c :

c=-0.122561+0.744862 I;

Set up a vector Initial condition with 250 points in each direction on the rectangle  [-1.3-1.3 I, 1.3 + 1.3I]:

zinit=Table[a+I b,{a,-1.3,1.3,(2.6/249)},{b,-1.3,1.3,(2.6/249)}];

Iterate starting from these initial conditions:

iterates=RecurrenceTable[{z(k+1)==c+((z(k))^2),z(0)==zinit},z,{k,0,20}];

Use ArrayPlot to show the fractal:

data=Map[Which[max((Abs[#1]))≤1,1,max((Abs[#1]))≤1.25,(1/5),max((Abs[#1]))≤1.5,(1/7),True,0]&,Flatten[iterates,(({{2}, {3}, {1}}))],{2}];

Dimensions[data]

ArrayPlot[(Transpose[data]),ColorRules->{1->Orange,(1/5)->Blue,0->Yellow,(1/7)->Red},ImageSize->400]

Smoothing of Initial Data in the Discretized Heat Equation

We can use RecurrenceTable to solve an initial-value problem for the discretized heat equation.

h=0.08;k=(1/320);

sol = RecurrenceTable[{v[i, j] == (k*(v[i - 1, j - 1] + v[i + 1, j - 1]))/h^2 + 
       (1 - (2*k)/h^2)*v[i, j - 1], v[i, 0] == 100*UnitBox[(i*h)/4]}, v, {i, -70, 70}, 
    {j, 0, 300}];

The following plots show the smoothing of the initial data as time progresses :

Quiet[Table[ListLinePlot[(Transpose[sol])[[i]],PlotRange->All,Filling->Axis,FillingStyle->Yellow],{i,{1,20,100,150,200,300}}]]

ListPlot3D[(Transpose[sol]),Mesh->None,PlotStyle->{Orange,Specularity[White,10]}]

Internals of Sum

Techniques of Summation

The following techniques are used for summation :

Table look - up for basic examples and for sums from standard references.

Algorithms for special classes such as polynomial, rational, hypergeometric term, PolyGamma sums, etc.

Special methods for divergent sums, Piecewise sums, holonomic - type sums, etc.

The Architecture of Sum

The Method Option

The following sum produces an answer in terms of the Binomial function by default :

(∑_(k=0)^n)((Binomial[n, k])^2)

This example is done using table look-up:

Sum[((Binomial[n, k])^2),{k,0,n},Method->"TableLookUp"]

We can confirm the result by Zeilberger's algorithm:

Sum[(((Binomial[n, k])^2)/(Binomial[2 n, n])),{k,0,n},Method->"HypergeometricTermZeilberger"]

Finally, Sum can be extended with any user-defined method :

SumSquares(expr_,{k_,min_,max_})/;IntegerQ[min]∧IntegerQ[max]∧max>min:=(∑_(k=min)^max)(expr^2)

SumSquares(___):=$Failed

Sum[k,{k,1,200},Method->"SumSquares"]

(∑_(k=1)^200)(k^2)

Summary

Sum, Product, RSolve, GeneratingFunction and other functions are well-equipped for solving symbolic problems of all types in Discrete Calculus.

New functions such as RecurrenceTable and FindSequenceFunction will allow users to conduct experiments in dynamical systems, combinatorics, etc.

Mathematica 7 provides a state-of-the art implementation of  Discrete Calculus functionality.