Derivation of Numerical Methods Using Computer Algebra.
The use of computer algebra systems in a course on scientific computation is demonstrated. Various examples, such as the derivation of Newton's iteration formula, the secant method, Newton--Cotes and Gaussian integration formulas, as well as Runge--Kutta formulas, are presented. For the derivations, the computer algebra system Maple is used.
Introduction
At ETH Z-urich we have redesigned our former courses on numerical analysis. We do not only run
numerical programs but also introduce the students to computer algebra and make heavy use of
computer algebra systems both in the lectures and the assignments.
Computer algebra may be used to generate numerical algorithms, to compute discretization errors,
to simplify proofs, etc., but also to run examples and to generate plots.
We claim, that it is easier for students to follow a derivation which is carried out with the help of a
computer algebra system than by hand. Computer algebra systems take over the hard hand work
such as e.g. solving systems of equations. Students do not need to be concerned with all the details
(and all the small glitches) of a manual derivation and can understand and keep the overview over
the general steps of the derivation. A computer supported derivation is also more convincing than
a presentation of the bare results without any reasoning.
Moreover, using computer algebra systems rather complex numerical formulas can be derived, far
more complex than what can be done in class by hand. E.g. all useful Newton-Cotes rules can be
computed without problems, in contrast to hand derivations, which usually end with Simpson's
rule.
We will prove these statements with some examples taken from our introductory courses in scientific
computing. We use Maple V Release 4, but the examples could also be reproduced e.g. with
Mathematica, MuPad or any other computer algebra system.
One of the first formulas that students learn is Newton's iteration to solve a nonlinear equation
Given an approximation x k for the root s, a better approximation x can be
obtained using the iteration function
This iteration converges quadratically to a single root s of f(x). This can be proven by computing
the first derivative of F (x) at which is zero.
dF := f(x) D (2) (f)(x)
(2)
D (2) (f)(s)
The arbitrary precision floating point arithmetic, which is provided by most computer algebra
systems, can be used to demonstrate what quadratic convergence means "in real life". As example,
we compute the square root of 9 using Newton's iteration to solve equation x starting
with As expected, the number of correct digits doubles with each iteration.
Digits := 80:
to 8 do xk := F(xk); lprint(xk); od:3.40000000000000000000000000000000000000000000000000000000000000000000000000000003.00009155413138017853055619134813458457312886243991760128175783932249942778667893.00000000000000000032526065174565133021986825552331682604822127820501315653354913.0000000000000000000000000000000000000000000000000000000000000000000000000000518
The convergence result is only valid for single roots, as the first derivative of f appears in the
denominator of (2), i.e. the result is only valid if f 0 The behavior of Newton's iteration for
an equation with multiple roots is the next topic we want to discuss using Maple. Let us assume
that f(x) has a zero of multiplicity n at We therefore define f(x) to be
where Again we inspect the first derivative of F (x). If F 0 then the iteration
converges only linearly.
x\Gammas
Taking the limit of the above expression for
We have just proven, that Newton's iteration converges linearly with the factor (n \Gamma 1)=n if f(x)
has a zero of multiplicity n. Thus e.g. convergence is linear with the factor 1=2 for double roots.
Newton's iteration also has a nice geometrical interpretation. Starting with the approximation x k ,
the next value x k+1 of the iteration is the intersection of the tangent at f(x) in [x k ; f(x k )] with
the x-axes. This property can also be proven with Maple. We set up an equation
for the tangent line. p(x) must interpolate [x k ; f(x k )] and must have the same derivative as f(x)
at these two conditions the parameters a and b of the tangent p(x) are defined.
We have claimed that the intersection of the tangent p(x) with the x-axes is the next Newton
approximation. If the equation solved again the iteration function (1) is obtained. This
proves that the geometrical interpretation is correct.
3 Secant Method
The secant method (3) approximates the derivative which appears in Newton's formula by a finite
difference.
As with Newton's method, we want to analyse the convergence with Maple. Using equation (3)
we obtain for the error e s the recurrence
s:
The right hand side can be expanded into a multivariate Taylor series at e
We assume that s is a single root, (f 0 must hold. We set
compute the first term of the Taylor series expansion:
(D (2)
If we divide this leading coefficient by e 0 and e 1 , we see that the limit of the quotient e 2 =e 0 =e 1
is a constant different from zero. We assume that the convergence coefficient is p and substitute
This equation is valid for all errors e 0 . Since the right hand side is constant, the left hand side must
also be independent of e 0 . This is only the case, if the exponent of the power of e0 is zero. This
condition is an equation for p which solved gives the well known convergence factor
5)=2
for the secant method.
Having considered the Newton and the secant method, we now want to derive and analyze a new
iteration method (a combination of the Newton and the secant method) to compute a root of a
nonlinear equation. The new method shall use the function values and the first derivatives at two
points. These four data define a degree three (Hermite) interpolation polynomial. A zero of this
polynomial can be taken as next approximation of the root. Unfortunately, the explicit expression
for this zero is rather complex, therefore we propose to use inverse interpolation for the given data.
We then only need to evaluate the resulting polynomial at to obtain the new approximation
for the root. In Maple, this is done with the following commands.
p(-f(x1))=x1, D(p)(-f(x1))=-1/D(f)(x1)-,
The resulting expression is still not very simple However, if the evaluation of f and f 0 is very
expensive, it may still pay off since the convergence rate is 2.73 as we will see. For the convergence
analysis we expand
s:
into a multivariate Taylor series at e as we have done for the secant method.
(D (3) )(f)(s) D(f)(s) (D (2) (D (2) (D (4)
As before with the secant and Newton method we consider only single roots, i.e., we assume that
condition holds, then the above equation tells us that in the limit
const
Let us again introduce the convergence coefficient p and make the following substitutions:
const
const
This equation must hold for all errors e 0 . Since K, p and const are all constant, the exponent of
must be zero.
[2:732050808; \Gamma:732050808]
Thus, the convergence factor is
3 and we have super-quadratic convergence.
With the help of Maple we want to demonstrate the above convergence rate for an example. We
use our algorithm to compute the zero of the function starting with x
1. For every iteration we print the number of correct digits (first column) and its ratio to the
number of correct digits in the previous step (second column). This ratio should converge to the
convergence rate 2:73. We see that this is the case.
Digits := 500:
for i to 6 do x2 := evalf(F(x0,x1));
else lprint(evalf(d2,20), evalf(d2/d1,20))
-7.8214175487946676893 2.8597847216407742880
-23.118600850923949542 2.9558070140989569391
-63.885801375143936026 2.7633939349141321183
-481.80650538330786806 2.7373103717659224742
5 Newton-Cotes Rules
Simpson's rule is a well known quadrature rule for approximating a definite integral
a
using three equidistant function values at the endpoints and in the middle of the integration
interval. The formula is obtained by interpolating the three function values and by computing the
integral of the interpolating polynomial of degree two. The following Maple statements may be
used to derive Simpson's rule. We first define a polynomial of degree two. We then state the
interpolation conditions and solve the resulting linear system for the coefficients of the polynomial.
Finally we integrate the polynomial and simplify the result.
What is the error of this integration rule? Discretization errors can be computed very simply by
appropriate series expansions. Let then we obtain
En :=
a
90 (D (4) )(f)(a) h 5
This shows that the error is proportional to h 5 . For the composite Simpson rule (n intervals of
length 2h, b \Gamma a = 2nh) the error therefore is b\Gammaa
in (a; b).
Instead of computing the interpolation polynomial, we can use the Maple function interp to
interpolate a polynomial through the points (0; y 0 ); (h; y 1 ) and (2h; y 2 ). By integration we obtain
the same result as above.
Notice that we have normalized the rule by dividing through the interval length.
More generally, we obtain Newton-Cotes Rules by interpolating values with given
equidistant nodes and by integrating the degree n interpolation polynomial. The following procedure
generates such a n + 1-point normalized Newton-Cotes Rule.
int(interp([seq(i*h, i=0.n)], [seq(y[i], i=0.n)], z), z=0.n*h))/(n*h):
With this procedure we can e.g. construct the Trapezoidal rule (n = 1), the Milne Rule
the Weddle Rule
For we obtain the following equidistant 9-point rule which is used by the Matlab function
quad8.
In [7] we can find one sided formulas that can also be generated with Maple in the same way. E.g.
given four equidistant function values (three intervals), find an approximation for the integral over
the 3rd interval:
6 Approximating Derivatives
When replacing derivatives by finite differences one uses relations like e.g.
They are obtained by computing derivatives of the corresponding interpolation polynomial. Relation
(4) is obtained by the Maple statement
Again, we can determine the discretization error with the help of a series expansion.
(D (2) )(y)(x)
Similarly, with the statements
an approximation for y 0 (x) is obtained. The discretization error of this approximation is also of
3 (D (3) )(y)(x)
4 (D (4) )(y)(x) h 3
7 Gauss-Quadrature
The idea of Gauss-Quadrature is to find nodes x i and weights w i so that the quadrature rule
is exact for polynomials of degree as high as possible. For we have to determine the six
unknowns We demand exact values for the integrals of the monomes x j
eqns
eqns := fw1 x1
We can solve this system with Maple:
5g
However, this brute force approach will not work for all values of n. For larger n the system of
nonlinear equations become too complicated for Maple. One has to add some more sophisticated
theory to compute the rules.
It is our goal to find nodes and weights to get an exact rule for polynomials of degree up to 2n \Gamma 1.
We can argue as follows: consider the decomposition of P 2n\Gamma1 by dividing by some polynomial
Qn (x) of degree n:
Applying rule (7) on both sides and subtracting yields the following expression for the error:
error :=
Now it is easy to see that we can make the error to zero by the following choices. First take Qn (x)
as the orthogonal polynomial on the interval [\Gamma1; 1] to the scalar product
By this choice and by the definition of orthogonal polynomial the first term in the error vanishes:
is a Legendre Polynomialavailable in Maple as orthopoly[P](n,x).
Second, choose as the nodes the (real) zeros of Qn . Then the second term in the error will also
Finally compute the weights according to Newton-Cotes by integrating the interpolation polynomial
for which is of course again Rn\Gamma1 by the uniqueness of the interpolation polynomial.
Thus Z 1
and the two last error term cancel.
So, we can compute a Gauss quadrature rule e.g. for with the following Maple statements:
\Gamma:1252334085; :1252334085; :3678314990; :5873179543; :7699026742; :9041172564;
We note that numerical errors occur (the weights should be symmetric) because we are computing
the rules here in a well-known unstable way. However, Maple offers us more precision by increasing
the value of Digits. With two runs of the above statements with different precision we are able
to obtain the rules correct to the amount of decimal digits we want.
The Gauss-Lobatto quadrature rule on [\Gamma1; 1] using the end points and two intermediate points
can be computed with a computer algebra system as follows ([1]). Considering the symmetry of
the formula, we make the ansatz
A
and require it to be exact for
If we like to compute a Kronrod extension by adding three more points (by symmetry one of them
will be 0), we could make the ansatz
A
and require exactness for
Thus our rule becomes
For a more elaborated treatment on generating Gauss Quadrature formulas symbolically we refer
to [8].
8 Generation of Explicit Runge-Kutta Formulas
In this section we show how a computer algebra system can be used to derive explicit Runge-Kutta
formulas. Such formulas are used to solve systems of differential equations of the first order. The
solution of the initial value problem
can be approximated by a Taylor series around x k , which is obtained from (8) by repeated differentiation
and replacing y 0 (x) by f
appears.
y (i)
@
@x
f
f
f y
The idea of the Runge-Kutta methods is to approximate the Taylor series (9) up to order m by
using only values of f
\Delta and no derivatives of it.
The general form of a s-stage explicit Runge-Kutta method is
a
s
where s is the number of "stages" and a i;j , b i and c i are real coefficients
To derive the coefficients of such a method, the series expansion of (9) and (10) are equated. This
leads to a set of nonlinear equations for the parameters a i;j , b i and c i which has to be solved.
For this derivation we have to compute the Taylor series expansions of (9) and (10). Maple knows
how to expand a function with two parameters that both depend on x, but we have to inform
Maple, that y 0 (x) has to be replaced by f
whenever it appears. We do this by overwriting
the derivative of the operator y.
In this result, D 1 (f)(x; y(x)) stands for the derivative of f with respect to the first argument, i.e.
\Delta . In order to make the result more readable, we define some alias substitutions.
We are now ready to derive the parameters of a Runge-Kutta formula for which is of order
3:
The variable TaylorPhi corresponds to \Phi in equation (9).
For the Runge-Kutta scheme we get the following Taylor series. Note that we keep the parameters
a i;j , b i and c i in symbolic form. RungeKuttaPhi corresponds to \Phi in equation (10).
The difference d of the two polynomials TaylorPhi and RungeKuttaPhi should be zero. We
consider d as a polynomial in the unknowns h, F, Fx, Fy, Fxx, etc. and set the coefficients of that
polynomial to zero. This gives us a nonlinear system of equations which has to be solved.
eqns := -coeffs(d, [h,F,Fx,Fy,Fxx,Fxy,Fyy])-;
eqns
\Gammab3
\Gammab3
For we found two (parameterized) solutions which can be represented by the following
coefficient schemes. Note that in the first solution, the unknowns c 3 and b 3 are free parameters,
i.e. they can take on any value. This is indicated by the entries in the
solution set. For the second solution, a 3;2 is a free parameter (see also Figure 1). From the second
solution we get the method of Heun of third order if we set a
Figure

1: Three Stage Runge-Kutta Methods of Order 3.
For further information on this topic we refer to [4].
9 Methods of Ritz and Galerkin
This section shows an example where the computer algebra systems takes over the hard hand work
and the students can understand the general steps of the method. Consider the boundary value
problem
For
2 x) the solution is
Equation (11) is obtained if we want to minimize the functional
In the method of Ritz the solution u is approximated by an ansatz (a linear combination of known
functions OE
Introducing y(x) in (12) we obtain a quadratic form in the unknown coefficients c j . Minimizing
this quadratic form gives us values for c j and an approximation y(x).
The principle is easy to describe but to compute a concrete example is rather tedious, even if we
choose only
and OE 2
Both ansatz-functions do satisfy the
boundary conditions.
With Maple we can compute and plot the approximation with the following statements:
x
Not every differential equation minimizes a functional. In the method of Galerkin one tries to
solve (11) also by the ansatz (13). The goal is to choose the coefficients c j such that the residual
small in some sense. For this we have to choose another set of functions
(x)g. The coefficients c j are now computed in such a way that the residual is orthogonal
to the space spanned by the functions /
We obtain this way again a system of linear equations for the coefficients c j .
x
Conclusions
In this paper we have given some examples on how we use computer algebra systems in our scientific
computing courses. Many numerical methods and classical proofs can be developed with only a
few statements in a computer algebra system.
However, it might be difficult sometimes to find the right ones, this method therefore is in general
still restricted to the reproduction of classical results. The use of a computer algebra system also
requires much experience, as it is not always easy to find an elegant way to the result one is
expecting.
We also made the experience that it may be particularly complicated to convince a computer
algebra system to perform a specific task. As an example take the convergence analysis of section 3.
We came across the expression (a were interested in the exponent if this expression is
written as a b . How is b obtained? Right, taking the logarithm to the base a. The result however
does not simplify, even not after using the simplify command.
a p a
a
What is the problem? Maple does not know as much as we do. Maple cannot simplify this
expression as it assumes a and p to be complex numbers. Obviously, a and p are real and positive
in our context, but Maple has to be informed about this fact using the assume facility. This
can be done directly in the simplify command for all the indeterminants which appear in the
expression to be simplified or for each unknown with the assume command. The symbol ~ signals
that assumptions have been made on a variable.
a is assumed to be real and positive
is assumed to be real
Look also at the discussion of the discretization error in chapter 6. We have computed a series
expansion of an expression comparable to the following one:
@x
Which leading term do you expect if this expression is expanded into a series? Maple gives you
the following answer:
@x
This result may be surprising for a student. The leading coefficient is indeed zero, but Maple does
not recognize this zero automatically. In general, it is a particularly difficult problem to recognize
zeros, but in this example the above result can be simplified using a special option to the command
simplify.
Computer algebra systems still have to progress. They are not yet a replacement for paper and
pencil.
We also do not want to conceal that computer algebra systems still have bugs and produce erroneous
results or results which are only valid under some assumptions. It is also very important to
demonstrate this fact to students. Students should learn, that results cannot blindly be trusted.
Whenever a numerical method is derived, the result has to be compared with ones expectations.
Nevertheless, a computer algebra system is a very powerful tool to be used in teaching numerical
methods. Many examples in this article have proven it. Further examples on the use of computer
algebra systems can be found e.g. in [3, 2, 5].



--R

Adaptive Quadrature - Revisited
The Billiard Problem

Symbolic Computation of Explicit Runge-Kutta Formulas
The Maple Technical Newsletter
Scientific Computing

chapter
--TR
finite elements;numerical methods;computer algebra;maple;rayleigh-ritz;galerkin method;runge-kutta method;nonlinear equations;quadrature formulas