Recherches sur l'intégration des équations différentielles aux différences finies, et sur leur usage dans la théorie des hasards.

P.S. Laplace


I imagine a solid composed of a number n of faces perfectly equal, and which I designate by the numbers 1, 2, 3 .. n ; I wish to have the probability that, by a number x of casts, I will bring about these n faces in sequence in the order 1, 2, 3 .. n .

Let y(x) be the probability. Laplace obtains the recurrence relation y(x) = y(x-1)-y(x-n)/(n^n)+1/(n^n) subject to the initial conditions that y(x) = 0 for x < n and y(n) = 1/(n^n) .

> restart:

> eqn12:=y(x)=y(x-1)-y(x-n)/n^n+1/n^n;

eqn12 := y(x) = y(x-1)-y(x-n)/(n^n)+1/(n^n)

Example 1. Laplace solves this in the case of n = 2 only. We must have y(0) = 0 and y(1) = 1/4 .

> eqn12_2:=subs(n=2,eqn12);

eqn12_2 := y(x) = y(x-1)-1/4*y(x-2)+1/4

> g2:=rsolve({eqn12_2,y(1)=0,y(2)=1/4},y(x),'makeproc'):

> g2(x);


When do we have an even wager that the sequence 12 will appear? When x = 3 .

> g2(3);


Example 2. This problem is more difficult with n = 3 . The initial conditions are that y(1) = 0 , y(2) = 0 , and y(3) = 1/27 . It is for reason of difficulty, I suspect, that Laplace stopped with n = 2 .

> eqn12_3:=subs(n=3,eqn12);

eqn12_3 := y(x) = y(x-1)-1/27*y(x-3)+1/27

> g3:=rsolve({eqn12_3,y(1)=0,y(2)=0,y(3)=1/27},y(x),'makeproc'):

> factor(normal(g3(x)));


This solution can be simplified and the use of limits be avoided. Suppose alpha is a root of z^3-27*z+27 . Then z-alpha divides z^3-27*z+27 and the factorization of this polynomial is


Thus the other two roots can be expressed in terms of the first as (alpha^2-18)/3 and (18-3*alpha-alpha^2)/3 Moreover, we may rewrite (z-alpha)/(z^3-27*z+27) as

9/((3*x+18-alpha^2)*(3*x-18+3*alpha+alpha^2)) .

We may then substitute this new expression for (z-alpha)/(z^3-27*z+27) in the solution given above.

The problem is completed if we then sum over the three roots of z^3-27*z+27 , substituting them successively for alpha . Since the polynomial is cubic, we can find the roots exactly. These roots are

> solve(z^3-27*z+27);


This prospect is not fun.

An alternative is to compute the values recursively . We do this in general with the procedure G(n,x) .

> G:=proc(n,x::nonnegint) options remember;

> if x<n then 0 elif x=n then 1/n^n else evalf(G(n,x-1)-G(n,x-3)/n^n+1/n^n); fi; end:

For example , let's consider a die with 6 sides. In how many tosses do we have an even wager that the faces will appear in consecutive order?

> evalf(G(6,32342));


> evalf(G(6,32343));