Cleve’s Corner: Cleve Moler on Mathematics and Computing
The MATLAB Blog
Guy on Simulink
MATLAB Community
Artificial Intelligence
Developer Zone
Stuart’s MATLAB Videos
Behind the Headlines
File Exchange Pick of the Week
Hans on IoT
Student Lounge
MATLAB ユーザーコミュニティー
Startups, Accelerators, & Entrepreneurs
Autonomous Systems
Quantitative Finance
MATLAB Graphics and App Building
Fun With The Pascal Triangle
The Wikipedia article on Pascal's Triangle has hundreds of properties of the triangle and there are dozens of other Web pages devoted to it. Here are a few facts that I find most interesting.
Blaise Pascal (1623-1662) was a 17th century French mathematician, physicist, inventor and theologian. His Traité du triangle arithmétique (Treatise on Arithmetical Triangle) was published posthumously in 1665. But this was not the first publication about the triangle. Various versions appear in Indian, Chinese, Persian, Italian and other manuscripts centuries before Pascal.
Contents
Blaise Pascal
Blaise Pascal (1623-1662) was a 17th century French mathematician, physicist, inventor and theologian. His Traité du triangle arithmétique (Treatise on Arithmetical Triangle) was published posthumously in 1665. But this was not the first publication about the triangle. Various versions appear in Indian, Chinese, Persian, Italian and other manuscripts centuries before Pascal.
Binomial Coefficients
The binomial coefficient usually denoted by ${n} \choose {k}$ is the number of ways of picking $k$ unordered outcomes from $n$ possibilities. These coefficients appear in the expansion of the binomial $(x+1)^n$. For example, when $n = 7$ syms x
n = 7;
x7 = expand((x+1)^n)
x7 = x^7 + 7*x^6 + 21*x^5 + 35*x^4 + 35*x^3 + 21*x^2 + 7*x + 1Formally, the binomial coefficients are given by $${{n} \choose {k}} = \frac {n!} {k! (n-k)!}$$ But premature floating point overflow of the factorials makes this an unsatisfactory basis for computation. A better way employs the recursion $$ {{n} \choose {k}} = {{n-1} \choose {k}} + {{n-1} \choose {k-1}}$$ This is used by the MATLAB function nchoosek(n,k).
Pascal Matrices
MATLAB offers two Pascal matrices. One is symmetric, positive definite and has the binomial coefficients on the antidiagonals.P = pascal(7)
P =
1 1 1 1 1 1 1
1 2 3 4 5 6 7
1 3 6 10 15 21 28
1 4 10 20 35 56 84
1 5 15 35 70 126 210
1 6 21 56 126 252 462
1 7 28 84 210 462 924
The other is lower triangular, with the binomial coefficients in the rows. (We will see why the even numbered columns have minus signs in a moment.)
L = pascal(7,1)
L =
1 0 0 0 0 0 0
1 -1 0 0 0 0 0
1 -2 1 0 0 0 0
1 -3 3 -1 0 0 0
1 -4 6 -4 1 0 0
1 -5 10 -10 5 -1 0
1 -6 15 -20 15 -6 1
The individual elements are
P(i,j) = P(j,i) = nchoosek(i+j-2,j-1)And (temporarily ignoring the minus signs) for i $\ge$ j
L(i,j) = nchoosek(i-1,j-1)The first fun fact is that L is the (lower) Cholesky factor of P.
L = chol(P)'
L =
1 0 0 0 0 0 0
1 1 0 0 0 0 0
1 2 1 0 0 0 0
1 3 3 1 0 0 0
1 4 6 4 1 0 0
1 5 10 10 5 1 0
1 6 15 20 15 6 1
So we can reconstruct P from L.
P = L*L'
P =
1 1 1 1 1 1 1
1 2 3 4 5 6 7
1 3 6 10 15 21 28
1 4 10 20 35 56 84
1 5 15 35 70 126 210
1 6 21 56 126 252 462
1 7 28 84 210 462 924
Pascal Triangle
The traditional Pascal triangle is obtained by rotating P clockwise 45 degrees, or by sliding the rows of L to the right in half increments. Each element of the resulting triangle is the sum of the two above it.triprint(L)
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
1 5 10 10 5 1
1 6 15 20 15 6 1
Square Root of Identity
When the even numbered columns of L are given minus signs the matrix becomes a square root of the identity.L = pascal(n,1) L_squared = L^2
L =
1 0 0 0 0 0 0
1 -1 0 0 0 0 0
1 -2 1 0 0 0 0
1 -3 3 -1 0 0 0
1 -4 6 -4 1 0 0
1 -5 10 -10 5 -1 0
1 -6 15 -20 15 -6 1
L_squared =
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1
Here is an exercise for you. What is sqrt(eye(n))? Why isn't it L?
Cube Root of Identity
When I first saw this, I was amazed. Rotate L counterclockwise. The result is a cube root of the identity.X = rot90(L,-1) X_cubed = X^3
X =
1 1 1 1 1 1 1
-6 -5 -4 -3 -2 -1 0
15 10 6 3 1 0 0
-20 -10 -4 -1 0 0 0
15 5 1 0 0 0 0
-6 -1 0 0 0 0 0
1 0 0 0 0 0 0
X_cubed =
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1
Sierpinski
Which binomial coefficients are odd? It's a fledgling fractal. odd = @(x) mod(x,2)==1;
n = 56;
L = abs(pascal(n,1));
spy(odd(L))
title('odd(L)')
Fibonacci
The sums of the antidiagonals of L are the Fibonacci numbers.n = 12; A = fliplr(abs(pascal(n,1))) for k = 1:n F(k) = sum(diag(A,n-k)); end F
A =
0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 1 1
0 0 0 0 0 0 0 0 0 1 2 1
0 0 0 0 0 0 0 0 1 3 3 1
0 0 0 0 0 0 0 1 4 6 4 1
0 0 0 0 0 0 1 5 10 10 5 1
0 0 0 0 0 1 6 15 20 15 6 1
0 0 0 0 1 7 21 35 35 21 7 1
0 0 0 1 8 28 56 70 56 28 8 1
0 0 1 9 36 84 126 126 84 36 9 1
0 1 10 45 120 210 252 210 120 45 10 1
1 11 55 165 330 462 462 330 165 55 11 1
F =
1 1 2 3 5 8 13 21 34 55 89 144
pi
The elements in the third column of lower triangular Pascal matrix are the triangle numbers. The n-th triangle number is the number of bowling pins in the n-th row of an array of bowling pins. $$t_n = {{n+1} \choose {2}}$$L = pascal(12,1); t = L(3:end,3)'
t =
1 3 6 10 15 21 28 36 45 55
Here's an unusual series relating the triangle numbers to $\pi$. The signs go + + - - + + - - .
pi - 2 = 1 + 1/3 - 1/6 - 1/10 + 1/15 + 1/21 - 1/28 - 1/36 + 1/45 + 1/55 - ...
type pi_pascal
function pie = pi_pascal(n)
tk = 1;
s = 1;
for k = 2:n
tk = tk + k;
if mod(k+1,4) > 1
s = s + 1/tk;
else
s = s - 1/tk;
end
end
pie = 2 + s;
Ten million terms gives $\pi$ to 14 decimal places.
format long
pie = pi_pascal(10000000)
err = pi - pie
pie =
3.141592653589817
err =
-2.398081733190338e-14
Matrix Exponential
Finally, I love this one. The solution to the (potentially infinite) set of ordinary differential equations $\dot{x_1} = x_1$ $\dot{x_j} = x_j + (j-1) x_{j-1}$ is $x_j = e^t (t + 1)^{j-1}$ This means that the matrix exponential of the simple diagonal matrixD = diag(1:7,-1)
D =
0 0 0 0 0 0 0 0
1 0 0 0 0 0 0 0
0 2 0 0 0 0 0 0
0 0 3 0 0 0 0 0
0 0 0 4 0 0 0 0
0 0 0 0 5 0 0 0
0 0 0 0 0 6 0 0
0 0 0 0 0 0 7 0
is
expm_D = round(expm(D))
expm_D =
1 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0
1 2 1 0 0 0 0 0
1 3 3 1 0 0 0 0
1 4 6 4 1 0 0 0
1 5 10 10 5 1 0 0
1 6 15 20 15 6 1 0
1 7 21 35 35 21 7 1
评论
要发表评论，请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。