Uses of the Spreadsheet in
Numerical Analysis
Jenchung 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 123 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 antiintellectual 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 nonresidues, 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 4x = 8cos^{4}x8cos^{2}x+1 

cos 5x = 16cos^{5}x20cos^{3}x+5cosx 

cos 6x = 32cos^{6}x48cos^{4}x+18cos^{2}x1 

the number pattern
emerges. The regularity of this pattern, that each number is the difference
of twice of its upperleft 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(n1)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 bruteforce
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:






2xT_{n}(x)  T_{n}_{1}(x),
n = 1,2,¼ 




This is only the beginning of a whole series of experiments with the spreadsheet.
Looking at the threeterm recursive relation as
2xT_{n}(x) = T_{n}_{+1}(x)+T_{n}_{1}(x), 

we find an interesting rule of generating the Chebyshev expansion
of the polynomial 2^{n1}x^{n}: if
2^{n1}x^{n} = c_{0}T_{0}(x)+c_{1}T_{1}(x)+¼+c_{n}T_{n}(x), 

and if
2^{n}x^{n}^{+1}
= d_{0}T_{0}(x)+d_{1}T_{1}(x)+¼+d_{n}_{+1}T_{n}_{+1}(x), 

then
d_{k} = c_{k}_{1}+c_{k}_{+1}for
k = 0,2,3,¼,n+1, 

and
here we assume c_{1} = c_{n}_{+1}
= c_{n}_{+2} = 0. Since an empty cell in a spreadsheet
always evaluates to 0, the coefficient d_{k} (k¹
1) can be specified by the simple rule of taking the sum of values of its
upper lefthand neighbor and the upper right hand neighbor. As for d_{1},
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
2^{6}x^{7} = T_{7}(x)+7T_{5}(x)+21T_{3}(x)+35T_{1}(x), 

for example.
While under the spreadsheet environment,
we might as well take advantage of the builtin graphics capability by
plotting the graph of the first few polynomials. To achieve this goal,
the numerical values forming the sequence
T_{n}(x_{k}),
k = 1,2,¼,100, 

where the sequence {x_{k}} 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:

The shifted Chebyshev polynomials are given by
T_{k}^{*}(x)
= T_{k}(2x1) 


Find the recursive relation satisfied by T_{k}^{*}.

Compute the coefficients of T_{k}^{*}
for k = 1, 2, ¼, 10.

Express 2^{2n1}x^{n} as a linear combination
of T_{k}^{*} for n
= 1,2,¼,12.

Chebyshev Polynomials of the Second Kind: If cos q
= x then
sin(n+1)q = sin
q
U_{n}(x) 

for some polynomial U_{n}.

Find the recurrence relations satisfied by U_{n}(x).

Compute the coefficients of U_{n}(x) for n
= 1,2,...,12.

Express 2^{n}x^{n} as a linear combination of U_{k}(x)
for n = 0,1,2,¼,12.

The shifted Chebyshev polynomials of the second kind are given by
U_{k}^{*}(x)
= U_{k}(2x1) 


Find the recursive relation satisfied by U_{k}^{*}.

Compute the coefficients of U_{k}^{*}
for k = 1, 2, ¼,10.

Express 2^{2n1}x^{n} as a linear combination
of U_{k}^{*} for n
= 1,2,¼,10.

Chebyshev Polynomials C_{n}(x): Let the polynomials
C_{n}(x)
be given recursively by
C_{0}(x) = 2 , C_{1}(x)
= x , C_{n}_{+1}(x) = xC_{n}(x)C_{n}_{1}(x)
for
n ³ 1. 


List the coefficients of C_{n}(x) for n =
0,1,2,¼,12.

Express x^{n} as a linear combination of C_{k}(x)
for n = 0,1,2,¼,12.

Chebyshev Polynomials S_{n}(x): Let the polynomials
S_{n}(x)
be given recursively by
S_{0}(x) = 1 , S_{1}(x)
= x , S_{n}_{+1}(x) = xS_{n}(x)S_{n}_{1}(x)
for
n ³ 1. 


List the coefficients of S_{n}(x) for n =
0,1,2,¼,12.

Express x^{n} as a linear combination of S_{k}(x)
for n = 0,1,2,¼,12.

x^{n }+ 1/ x^{n} as a function of x+^{1}/_{x}:
Let y_{n} denote the rational function x^{n}+
1/
x^{n} .

Find a recurrence relation satisfied by y_{n} for n³
0.

Express y_{n} as a function of y_{1}( = x+^{1}/_{x})
for n = 1,2,...,20.

Trace of A^{n}: Let A be a twobytwo 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 a^{n}
and b^{n} are the eigenvalues of A^{n}, it
follows that the determinant of A^{n} is a^{n}b^{n}
while the trace T_{n} of A^{n} is T_{n}
= a^{n}+b^{n}.

Find the recurrence relations satisfied by T_{n} for n
³
0.

Express T_{n} as a function of t and m for
n
= 0,1,2,...,20.

Determinant of a tridiagonal matrix: Let M_{n} be
an n×n tridiagonal matrix of the form
M_{n} = 
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û 
. 


Find a recursive formula satisfied by the determinant d_{n}(x)
of M_{n}.

Compute d_{n}(x) for n = 2,3,...,12.

Powers of a twobytwo matrix:

Find the pattern of
for n = 1,2,3,¼.

Use the spreadsheet to compute its expansion.
Divided Differences
Let f be an arbitrary function.
For any two distinct points x_{0},x_{1} set
[ x_{0}x_{1}] = 
f(x_{0})f(x_{1})
x_{0}x_{1} 
. 

For any three distinct points x_{0},x_{1},x_{2}
set
[ x_{0}x_{1}x_{2}]
= 
[ x_{0}x_{1}] [x_{1}x_{2}]
x_{0}x_{2} 
. 

In general if x_{0},x_{1},¼,x_{n}
are distinct, set
[ x_{0}x_{1}x_{2}¼x_{n}]
= 
[x_{0}x_{1}x_{2}¼x_{n}_{1}]
[ x_{1}x_{2}x_{3}¼x_{n}]
x_{0}x_{n} 
. 

It is a matter of algebraic substitution to verify that for a polynomial
f(x)
of degree n the identity



f(x_{0})+[ x_{0}x_{1}]
(xx_{0})+[x_{0}x_{1}x_{2}]
(xx_{0})(xx_{1}) 



+¼+[ x_{0}x_{1}x_{2}¼x_{n}](xx_{0})(xx_{1})¼(xx_{n}_{1}) 




holds for any distinct x_{0},x_{1},x_{2},¼,x_{n}.
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:
x_{0} x_{1} x_{2
}x_{3}
x_{1} x_{2} x_{3}
x_{2} x_{3}
x_{3}
f(x_{0}) [x_{0}
x_{1}]
[x_{0}x_{1}x_{2}]
[x_{0}x_{1}x_{2}x_{3}]
f(x_{1}) [x_{1}
x_{2}]
[x_{1}x_{2}x_{3}]
f(x_{2}) [x_{2}
x_{3}]
f(x_{3})
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:

listing of the sequence of numbers in the domain of f ,

entering only one formula,

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 z_{1},z_{2},¼,z_{n}
be the zeros of a monic polynomial p(t) = t^{n}+a_{1}t^{n}^{1}+¼+a_{n}.
Newton discovered a way to compute the sums
s_{k} = z_{1}^{k}+z_{2}^{k}+¼+z_{n}^{k},
k = 1,2,3,¼ 

without finding z_{1},z_{2},¼,z_{n}
explicitly by using the formula
s_{k} = (s_{k}_{1}a_{1}+s_{k}_{2}a_{2}+¼+s_{1}a_{k}_{1}+ka_{k}), 

with the convention that a_{k} = 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 x^{4}+x^{3}+2x^{2}2
is considered:
Leverrier's Method
Let A be an nbyn matrix.
Let the characteristic polynomial of A be given by
t^{n}+a_{1}t^{n}^{1}+a_{2}t^{n}^{2}+¼+a_{n}. 

If the eigenvalues of A are z_{1},z_{2},¼,z_{n},
then z_{1}^{k},z_{2}^{k},¼,z_{n}^{k}
are the eigenvalues of A^{k} and so
s_{k} = z_{1}^{k}+z_{2}^{k}+¼+z_{n}^{k}
= tr(A^{k}). 

It follows from Newton's formula for sums of powers of zeros that






 
1
n 
(s_{n}+a_{1}s_{n}_{1}+a_{2}s_{n}_{1}+¼+a_{n}_{1}s_{1}),
n ³ 2. 




The coefficients of the characteristic polynomial then be computed from
a table as in this diagram:
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 fourbyfour
matrix

é
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
û 


is
t^{4}5t^{3}3t^{2}+126t+98. 

Solution of a Tridiagonal System
Consider the tridiagonal system
of equations:

é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û 

é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û 
= 
é
ê
ê
ê
ê
ê
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
ú
ú
ú
ú
ú
û 
. 

The solution
can be computed as follows:
b_{1} = b_{1
}g_{1} = d_{1} 

b_{2} = b_{2} 
a_{2}c_{1
}b_{1} 
g_{2}
= d_{2} 
a_{2}g_{1
}b_{1} 


b_{3} = b_{3} 
a_{3}c_{2
}b_{2} 
g_{3}
= d_{3} 
a_{3}g_{2
}b_{2} 


b_{k}
= b_{k} 
a_{k}c_{k}_{1
}b_{k1} 
g_{k}
= d_{k} 
a_{k}g_{k1
}b_{k1} 


b_{n}
= b_{n} 
a_{n}c_{n}_{1
}b_{n1} 
g_{n}
= d_{n} 
a_{n}g_{n1
}b_{n1} 


The solution x_{k} is then computed ''backward''.
x_{n}_{1} = 
g_{n1}c_{n}_{1}x_{n
}b_{n1} 


x_{n}_{2} = 
g_{n2}c_{n}_{2}x_{n}_{1
}b_{n2} 


x_{k} = 
g_{k}c_{k}x_{k}_{+1
}b_{k} 


x_{1} = 
g_{1}c_{1}x_{2
}b_{1} 
. 

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?
Answer: the bisection method.
The process of bisection method
can be implemented on a spreadsheet when the builtin @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 nth power M^{n} converges
to a matrix having the same value for entries on the same column. The phenomenon
of the convergence of M^{n} 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 reexamined and made suitable
for the design of automobiles. Bezier modified Bernstein's formula and
considered the expression given by
for any given vectors Q_{0},Q_{1},¼,Q_{n}
(taken from the same vector space). This curve P(t) is now
called the BernsteinBezier polynomial curve. We see at once that
the BernsteinBezier polynomial curve enjoys the following properties:

P(0) = Q_{0}, P(1) = Q_{n}.

For 0 £ t £
1, the vector P(t) lies in the convex hull of the vectors
Q_{0},Q_{1},¼,Q_{n}.

P^{¢}(0) = n(Q_{1}Q_{0}),P^{¢}(1)
= n(Q_{n}Q_{n}_{1})
Statements 1 and 3 say that this curve passes through the points Q_{0}
and Q_{n} with the tangents in the directions of Q_{0}Q_{1}
and Q_{n}_{1}Q_{n}. In case n
= 3 the resulting curve takes the form
P(t) = (1t)^{3}Q_{0}+3(1t)^{2}tQ_{1}+3(1t)t^{2}Q_{2}+t^{3}Q_{3} 

and is called the Bezier cubic curve. This may be expressed in the
following matrix notation:
P(t) = 
é
ë 



ù
û 

é
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
û 

é
ê
ê
ê
ê
ê
ë 



ù
ú
ú
ú
ú
ú
û 


and statements 1 and 3 above become
P(0) = Q_{0} , P(1)
= Q_{3} , P^{¢}(0)
= 3(Q_{1}Q_{0}) , P^{¢}(1)
= 3(Q_{3}Q_{2}). 

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 computeraided
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

M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover,
1967.

E. W. Cheney, Numerical Mathematics and Computing, 3rd Ed., Brooks/Cole,
1994.

G. Farin, Curves and Surfaces for Computer Aided Geometric Design, A Practical
Guide, Academic Press, 1988.

U. Grenander, Mathematical Experiments on the Computer, Academic Press,
1982.

D. M. Wiberg, Schaum's Outline of Theory and Problems of State Space and
Linear Systems, McGrawHill, 1971.

H.S. Wilf, Mathematics for the Physical Sciences, Dover, 1962.