Uses of the Spreadsheet in Numerical Analysis

Jen-chung Chuan

Department of Mathematics

National Tsing Hua University

jcchuan@math.nthu.edu.tw

Introduction

When the personal computer arrived in the early 1980's, there were very few appropriate software programs available for computation. Lotus 1-2-3 first took advantage of the PC environment, by making computations so convenient that using the spreadsheet was the main reason to turn on the computer. Right now the PC is often performing more anti-intellectual computations than people are aware of. Yet the spreadsheet still provides the rare kind of computational environment that reflects the user's mathematical sophistication.

Numerical analysis is a branch of discipline that people perform computations in order to see the insight of mathematics. To seek the insight of the natural numbers, Euler started by examining the first 99 numbers and calculated the sum of the divisors. He then experimented with the results, looking for regularities, and gradually strengthened the belief in a conjecture about the sums. In the process of discovering the law of quadratic reciprocity, Gauss set up huge tables, of prime numbers, of quadratic residues and non-residues, and of the fractions 1/p for p = 1 to p = 1000 with their decimal expansions carried out to a complete period, and therefore sometimes to several hundred places. The spirit of Euler and Gauss still lives, if the appropriate technology is used. The spreadsheet offers one such choice.

In the undergraduate numerical analysis course, the spreadsheet may be applied to experiments in topics that deal with number patterns, especially patterns associated with the recursion. We begin by discussing one such possibility in detail.

Chebyshev Polynomials

From the formula of the cosine expansion
 cos 0x = 1
 cos x = cosx
 cos 2x = 2cos2x-1
 cos 3x = 4cos3x-3cosx
 cos 4x = 8cos4x-8cos2x+1
 cos 5x = 16cos5x-20cos3x+5cosx
 cos 6x = 32cos6x-48cos4x+18cos2x-1
the number pattern

emerges. The regularity of this pattern, that each number is the difference of twice of its upper-left neighbor and the one located two entries above, reflects the recursive application of the trigonometric formula

 cos(n+1)x = 2cos x cos nx -cos(n-1)x.

This curiosity naturally leads us to considering the generation of number pattern using the spreadsheet, an activity that develops our mental power of abstraction. If the trigonometric formula is regarded as a concrete experience from the secondary school mathematics, then the ability of realizing the recursive nature by detaching the powers of cosx to focus solely on the coefficients, requires mathematical maturity beyond that of brute-force computation. Computations with the spreadsheet keep our interests in this number pattern and confirm the validity of our conjecture raised while observing the regularity. After playing around with this number pattern for a while, it becomes appropriate to introduce the Chebyshev polynomials (of the first kind), a sequence of functions sharing the identical number pattern. Again the regularity suggests the following recursive relationship bonding the polynomials:

 T0(x)
 =
 1,T1(x) = x
 Tn+1(x)
 =
 2xTn(x) - Tn-1(x),   n = 1,2,¼
This is only the beginning of a whole series of experiments with the spreadsheet. Looking at the three-term recursive relation as
 2xTn(x) = Tn+1(x)+Tn-1(x),
we find an interesting rule of generating the Chebyshev expansion of the polynomial 2n-1xn: if
 2n-1xn = c0T0(x)+c1T1(x)+¼+cnTn(x),
and if
 2nxn+1 = d0T0(x)+d1T1(x)+¼+dn+1Tn+1(x),
then
 dk = ck-1+ck+1for k = 0,2,3,¼,n+1,
and
 d1 = 2c0+c2;
here we assume c-1 = cn+1 = cn+2 = 0. Since an empty cell in a spreadsheet always evaluates to 0, the coefficient dk  (k¹ 1) can be specified by the simple rule of taking the sum of values of its upper left-hand neighbor and the upper right hand neighbor. As for d1, just modify the formula to take twice the value of the upper left hand neighbor instead.
The interpretation of this table is now clear: the last row indicates
 26x7 = T7(x)+7T5(x)+21T3(x)+35T1(x),
for example.

While under the spreadsheet environment, we might as well take advantage of the built-in graphics capability by plotting the graph of the first few polynomials. To achieve this goal, the numerical values forming the sequence

 Tn(xk),  k = 1,2,¼,100,
where the sequence {xk} comprises an evenly spaced nodes in the interval  [-1,1] , are to be generated. Again the recursive relationship allows a quick way of doing it, thanks to the COPYING command.
This way, the spreadsheet reveals interesting graphic patterns as well as interesting number patterns.

Further Investigations

There is a whole series of exercises patterned after the experiment with Chebyshev polynomials describes above:
1. The shifted Chebyshev polynomials are given by
2.  Tk*(x) = Tk(2x-1)
1. Find the recursive relation satisfied by Tk*.
2. Compute the coefficients of Tk* for k = 1, 2, ¼, 10.
3. Express 22n-1xn as a linear combination of Tk* for n = 1,2,¼,12.
3. Chebyshev Polynomials of the Second Kind: If cos q = x then
4.  sin(n+1)q = sin q Un(x)
for some polynomial Un.
1. Find the recurrence relations satisfied by Un(x).
2. Compute the coefficients of Un(x) for n = 1,2,...,12.
3. Express 2nxn as a linear combination of Uk(x) for n = 0,1,2,¼,12.
5. The shifted Chebyshev polynomials of the second kind are given by
6.  Uk*(x) = Uk(2x-1)
1. Find the recursive relation satisfied by Uk*.
2. Compute the coefficients of Uk* for k = 1, 2, ¼,10.
3. Express 22n-1xn as a linear combination of Uk* for n = 1,2,¼,10.
7. Chebyshev Polynomials Cn(x): Let the polynomials Cn(x) be given recursively by
8.  C0(x) = 2 , C1(x) = x , Cn+1(x) = xCn(x)-Cn-1(x) for n ³ 1.
1. List the coefficients of Cn(x) for n = 0,1,2,¼,12.
2. Express xn as a linear combination of Ck(x) for n = 0,1,2,¼,12.
9. Chebyshev Polynomials Sn(x): Let the polynomials Sn(x) be given recursively by
10.  S0(x) = 1 , S1(x) = x , Sn+1(x) = xSn(x)-Sn-1(x) for n ³ 1.
1. List the coefficients of Sn(x) for n = 0,1,2,¼,12.
2. Express xn as a linear combination of Sk(x) for n = 0,1,2,¼,12.
11. xn + 1/ xn as a function of x+1/x: Let yn denote the rational function xn+ 1/ xn .
1. Find a recurrence relation satisfied by yn for n³ 0.
2. Express yn as a function of y1( = x+1/x) for n = 1,2,...,20.
12. Trace of An: Let A be a two-by-two matrix. Suppose that A has eigenvalues a and b. It is known that the determinant m of A is m = ab and the trace t of A is t = a+b. Since an and bn are the eigenvalues of An, it follows that the determinant of An is anbn while the trace Tn of An is Tn = an+bn.
1. Find the recurrence relations satisfied by Tn for n ³ 0.
2. Express Tn as a function of t and m for n = 0,1,2,...,20.
13. Determinant of a tridiagonal matrix: Let Mn be an n×n tridiagonal matrix of the form
14. Mn é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
 x
 1
 1
 x
 1
 1
 x
 1
 ¼
 1
 x
 1
 1
 x
 1
 1
 x
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û
.
1. Find a recursive formula satisfied by the determinant dn(x) of Mn.
2. Compute dn(x) for n = 2,3,...,12.
15. Powers of a two-by-two matrix:
1. Find the pattern of

2. for n = 1,2,3,¼.
3. Use the spreadsheet to compute its expansion.

Divided Differences

Let f be an arbitrary function. For any two distinct points x0,x1 set
 [ x0x1] = f(x0)-f(x1) x0-x1 .
For any three distinct points x0,x1,x2 set
 [ x0x1x2] = [ x0x1] -[x1x2] x0-x2 .
In general if x0,x1,¼,xn are distinct, set
 [ x0x1x2¼xn] = [x0x1x2¼xn-1] -[ x1x2x3¼xn] x0-xn .
It is a matter of algebraic substitution to verify that for a polynomial f(x) of degree n the identity
 f(x)
 =
 f(x0)+[ x0x1] (x-x0)+[x0x1x2] (x-x0)(x-x1)
 +¼+[ x0x1x2¼xn](x-x0)(x-x1)¼(x-xn-1)
holds for any distinct x0,x1,x2,¼,xn. The spreadsheet can then be applied to solve the polynomial interpolation problem based on this formula of divided differences. The basic procedure is to set up a table to compute the divided differences according to this scheme:

x0 x1 xx3

x1 x2 x3

x2 x3

x3

f(x0)   [x0 x1] [x0x1x2] [x0x1x2x3]

f(x1)   [x1 x2] [x1x2x3]

f(x2)   [x2 x3]

f(x3)

Having completed this scheme, we are ready to compute the numerical value of the interpolating polynomial according to Horner's algorithm. Working under the spreadsheet, this requires three steps:

1. listing of the sequence of numbers in the domain of f ,
2. entering only one formula,
3. copying the formula.
Taking advantage of the interactive graphics supplied by the spreadsheet, the shape of the function can be visualized at once.

When blessed with a ``fast'' CPU, any change in the initial data causes an instant redisplay of the graph of the resulting interpolating polynomial. In this experiment an activity combining the understanding of algebra, algorithm, and graph plotting takes place when the spreadsheet is appropriately manipulated.

Sums of Powers of Zeros

Let z1,z2,¼,zn be the zeros of a monic polynomial p(t) = tn+a1tn-1+¼+an. Newton discovered a way to compute the sums
 sk = z1k+z2k+¼+znk,   k = 1,2,3,¼
without finding z1,z2,¼,zn explicitly by using the formula
 sk = -(sk-1a1+sk-2a2+¼+s1ak-1+kak),
with the convention that ak = 0 for k > n. It takes nothing more than the understanding of ``relative address'' and ``absolute address'' of a cell in the spreadsheet to implement this recursive relationship. In this example the polynomial x4+x3+2x2-2 is considered:

Leverrier's Method

Let A be an n-by-n matrix. Let the characteristic polynomial of A be given by
 tn+a1tn-1+a2tn-2+¼+an.
If the eigenvalues of A are z1,z2,¼,zn, then z1k,z2k,¼,znk are the eigenvalues of Ak and so
 sk = z1k+z2k+¼+znk = tr(Ak).
It follows from Newton's formula for sums of powers of zeros that
 a1
 =
 -s1
 an
 =
 - 1 n (sn+a1sn-1+a2sn-1+¼+an-1s1),  n ³ 2.
The coefficients of the characteristic polynomial then be computed from a table as in this diagram:
 s1
 s2
 s3
 s4
 a1
 a1s1
 a1s2
 a1s3
 a1s4
 a2
 a2s1
 a2s2
 a2s3
 a2s4
 a3
 a3s1
 a3s2
 a3s3
 a3s4
 a4
 a4s1
 a4s2
 a4s3
 a4s4
The actual implementation of this algorithm based on the spreadsheet requires the skill of evaluating the powers of a matrix, and the ability to set up a multiplication table.
This computation shows that the characteristic polynomial of the four-by-four matrix
é
ê
ê
ê
ê
ê
ë
 2
 4
 5
 -1
 3
 1
 3
 2
 -1
 -3
 4
 3
 0
 1
 1
 -2
ù
ú
ú
ú
ú
ú
û
is
 t4-5t3-3t2+126t+98.

Solution of a Tridiagonal System

Consider the tridiagonal system of equations:
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
 b1
 c1
 a2
 b2
 c2
 a3
 b3
 c3
 ¼
 an
 bn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
 x1
 x2
 x3
 ·
 ·
 xn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë
 d1
 d2
 d3
 ·
 ·
 dn
ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û
The solution

can be computed as follows:
 b1 = b1     g1 = d1
 b2 = b2- a2c1 b1 g2 = d2- a2g1 b1
 b3 = b3- a3c2 b2 g3 = d3- a3g2 b2
 ...
 bk = bk- akck-1 bk-1 gk = dk- akgk-1 bk-1
 ...
 bn = bn- ancn-1 bn-1 gn = dn- angn-1 bn-1
The solution  xk is then computed ''backward''.
 xn = gn bn
 xn-1 = gn-1-cn-1xn bn-1
 xn-2 = gn-2-cn-2xn-1 bn-2
 ...
 xk = gk-ckxk+1 bk
 ...
 x1 = g1-c1x2 b1 .

In the spreadsheet solution, the given data are placed in columns A,B,C,D, the intermediate calculations of b and g are placed in column E and F. The solution x is placed in column G. Note that the formulas are entered in the order: E1, F1, E2, F2, E3..F5, G5, G1..G4.

Bisection Method

Question: What's the most practical application of the Intermediate Value Theorem?

The process of bisection method can be implemented on a spreadsheet when the built-in @IF function is used iteratively.

Stochastic Matrix

A square matrix M with nonnegative entries and with row sum equal to 1 is known as a stochastic matrix. For such a matrix its n-th power Mn converges to a matrix having the same value for entries on the same column. The phenomenon of the convergence of Mn can be illustrated with the spreadsheet. After entering the matrix and a column vector with all entries zeros except one, the crucial formula, reflecting matrix multiplication, is then entered and then copied.

The power method of finding the leading eigenvalue of a matrix can be carried out in a similar fashion.

Bezier Cubic Curve

For any t and any positive integer n, we have the identity
This formula expresses the fact that the function constant 1 is a sum of the n+1 polynomial functions
Trying to find a polynomial function approximating a continuous function uniformly, S.N. Bernstein came up with the idea of constructing the Bernstein polynomial
for a given continuous function f defined on [ 0,1] . This formula expresses the fact that at each point t the value p(t) is a suitable convex combination of the fixed values f(0),f(1/n),¼,f(1). It took almost 70 years before this formula was re-examined and made suitable for the design of automobiles. Bezier modified Bernstein's formula and considered the expression given by
for any given vectors Q0,Q1,¼,Qn (taken from the same vector space). This curve P(t) is now called the Bernstein-Bezier polynomial curve. We see at once that the Bernstein-Bezier polynomial curve enjoys the following properties:
1. P(0) = Q0, P(1) = Qn.
2. For 0 £ t £ 1, the vector P(t) lies in the convex hull of the vectors Q0,Q1,¼,Qn.
3. P¢(0) = n(Q1-Q0),P¢(1) = n(Qn-Qn-1)
Statements 1 and 3 say that this curve passes through the points Q0 and Qn with the tangents in the directions of Q0Q1 and Qn-1Qn. In case n = 3 the resulting curve takes the form
 P(t) = (1-t)3Q0+3(1-t)2tQ1+3(1-t)t2Q2+t3Q3
and is called the Bezier cubic curve. This may be expressed in the following matrix notation:
P(t) =  é
ë
 1
 t
 t2
 t3
ù
û
é
ê
ê
ê
ê
ê
ë
 1
 0
 0
 0
 -3
 3
 0
 0
 3
 -6
 3
 0
 1
 3
 -3
 1
ù
ú
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ê
ë
 Q0
 Q1
 Q2
 Q3
ù
ú
ú
ú
ú
ú
û

and statements 1 and 3 above become

 P(0) = Q0 , P(1) = Q3 , P¢(0) = 3(Q1-Q0) , P¢(1) = 3(Q3-Q2).
It is thus an easy matter to construct the Bezier cubic curve once the initial position P(0), the final position P(1), the initial velocity P¢(0) and the final velocity P¢(1) are given.

The basic task of a computer-aided design (CAD) program is to draw some smooth curve that meets the requirement of passing through two given points and be tangent to specified given vectors there. The Bezier cubic curve discussed here furnishes a neat solution to the problem. This application of linear algebra concept (convex combinations) to CAD can aptly be demonstrated with a spreadsheet program. We arrange the spreadsheet so the initial position P(0), the final position P(1), the final velocity P¢(1) are ``fixed'' and the initial vector P¢(0) is allowed to ``vary'' along the ``unit circles'', so each time the RECALC key is pressed a cubic curve corresponding to the solution of a new initial velocity is displayed.

Conclusion

The best kind of software lures users to explore possibilities that were NOT designed for in the first place. Spreadsheet belongs to this category.

References

1. M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover, 1967.
2. E. W. Cheney, Numerical Mathematics and Computing, 3rd Ed., Brooks/Cole, 1994.
3. G. Farin, Curves and Surfaces for Computer Aided Geometric Design, A Practical Guide, Academic Press, 1988.
4. U. Grenander, Mathematical Experiments on the Computer, Academic Press, 1982.
5. D. M. Wiberg, Schaum's Outline of Theory and Problems of State Space and Linear Systems, McGraw-Hill, 1971.
6. H.S. Wilf, Mathematics for the Physical Sciences, Dover, 1962.