University of Nottingham School of Mathematical Sciences
MATH4063 Scientific Computing and C
Submission Date: Monday, 4th November 2019 Assessed Coursework 1
The following questions are to be used for the coursework assessment in the module MATH4063. Credit will be given for relevant working associated with solving the problems, as well as the answer. A single zip or tar file containing your solution should be submitted through the module webpage on Moodle as per the instructions on the Coursework Submission form.
The style and efficiency of your programs is important. A barelypassing solution will include programs which include some of the correct elements of a solution. A borderline distinction solution will include working code that neatly and efficiently implements the relevant algorithms, and that shows evidence of testing.
An indication is given of the approximate weighting of each question by means of a figure enclosed by square brackets, eg 12.
All calculations should be done in double precision.
1. The Horner method is the standard method for evaluating polynomials. Suppose that a poly nomial is given in the form
pnxan an1xa0xn. 1 The Horner method is based on the observation that pn can be rewritten in the form
pnxan xan1 xan2 xa1 xa0. 2 Then, we can evaluate pn at a point z recursively, working our way from the innermost paren
thesis in the expression above:
bz0a0 3 bzi ai bzi1z, i1,2,…,n, 4
as we have bznpnz. Apart from providing a stable and efficient evaluation method1, Horner method also provides a convenient way to compute the roots of the polynomial.
Here we are going to combine it with Newton method to find a real solution to pnx0: Fix z0 and iterate, for k0, zk1zkpnzkpnzk,
until a given tolerance is reached. To this end, define the associated polynomial qz xbzbz xbzxn1.
n1 n1 n2 0
Then, from the uniqueness property of polynomial division, it can be seen that
p xbz xzqz x. 5 n n n1
Differentiating 5 we get pnxqnz1xxzqnz1x and hence pnzqnz1z. In other words, Horner Method also gives an efficient way to compute derivatives of polynomials! Using this formula we see that the Newton iteration reduces to:
Fixz anditeratefork0, z z p z qz z , 6 0 k1 k n k n1 k
1It is indeed more efficient than the obvious method of evaluation based on 1, using 2n flops instead of 3n flops.
until a given tolerance is reached. This is the NewtonHorner method: at each iteration, all is required is to use pnzkbzn and to evaluate the associated polynomial qnz1 at the point zk, which we can do once again using Horner method!
After a root z has been found, then 0pnzbzn and so
p xxzqz x, thatis qz xp xxz. 7
n n1 n1 n
So, we can proceed to find a subsequent root by working on the quotient polynomial qnz1, which is again just the associated polynomial, and so on. Note that here we are assuming that all roots are real. Every time one root has been found, the new polynomial has one degree less: this is called deflation procedure.
a In a file newtonhorner.cpp, write a C function called horner which, given the n1 coefficients defining a polynomial, evaluates the n1 Horner coefficients at a point x using the formulas 34. Write a test program called q1a.cpp to evaluate the polynomial pxx32×23 at x0.5 by calling the function horner and output the computed value. 10
b Add to your file newtonhorner.cpp, a new C function newton which applies the NewtonHorner method 6 to find a root z of a polynomial starting from some initial guess z0. The function should iterate 6 while zk1zkTOL and the number of iterations is less then some nmax. Write a test program called q1b.cpp to perform the following tasks. 1 Fix a polynomial degree, the polynomial coefficients, and the initial guess for the root. The coefficients should be stored into a dynamically allocated array. 2 Calculates a root of the given polynomial by calling the function newton once. 3 Outputs the computed root. 4 Destroy storage. Set TOL10e8, nmax5, z04. Test your program with the polynomial pxx2 1 and check that your program outputs the root z1. Now rerun your program with pxx5 x4 9×3 x2 20×12 What number do you get at output? Now set nmax10, rerun your program and report again the output. 15
c Write a test program called q1c.cpp that calculates all the roots of a polynomial by following the deflation procedure described above. On each step, your program should call the function newton from part b to compute a root. Set TOL10e8, nmax20. For the first iteration use z04 while for subsequent iterations you may use the previous root as initial guess. Test your program with the polynomial pxx21 and check that your program correctly outputs the two roots z1. Now change the polynomial topxx5 x4 9×3 x2 20x12andreportthefiverootsinoutput. 10
2. The Chebyshev polynomials are defined for x1, 1 by
Pnxcoskarccosx, n0,1,2,… 8
Although written in this form they may not seem like polynomials, they are. This can be shown by proving that they can be equally be generated using the following Chebyshev recursion formula :
P0 x
P1 x Pn1x
1,
x,
2 x PnxPn1x, n1.
They possess the property that they are orthogonal to each other on the interval 1, 1 with respect to the weight function wx1x212; more precisely,
1 , PnxPmxwx dx2,
if nm0, if nm0, otherwise.
1
0,
2
As all sets of orthogonal polynomials, they provide an effective way to represent functions as power series expansions. It is therefore useful to have an efficient algorithm to evaluate the partial sums
n
pnx aiPix. 9
i0
This can be done using the Clenshaw algorithm which is used to calculate sums of functions that are defined recursively. A special case of the Clenshaw algorithm is the Horner method of Question 1! When applied to the Chebyshev sum 9, the Clenshaw algorithm reads as follows: giventhen1coefficientsof 9,computerecursivelycoefficientsbxn2,bxn1,…,bx0 as:
bxn1 bxn20
bxiai2xbxi1xbxi2x forin,n1,…,1.
The desired sum is then obtained as pnxa0xbx1bx2.
a In a file chebyshev.cpp write a C function called chebeval which calculates the Chebyshev polynomial Pnx, n0, using the Chebyshev recursion formula given above. Write a test program q2a.cpp to compute the following values: P20.5, P30.5 and P9 0.7. 10
b Add to your file chebyshev.cpp a C function called chebsum which calculates 9 using the Chebyshev recursion formula given above. Your function should input an integer n, an array of n1 coefficients ai, i1,…,n, and the value x at which the sum is to be computed and return the value Snx. Write a test program called q2b.cpp to computesuchsumforn3,a0 a1 a2 a3 1andx0.5andcheckthatthe result corresponds to that obtained with pencil and paper. Then let n10, and aii, i0,1,…,10, and output the sum for x1,0.5,0,0.5,1. 10
3. The Chebyshev polynomials have a number of great properties. One of them is that they yield the socalled Chebyshev points which are of great importance, for instance, for the interpo lation problem: given a set of interpolation nodes x0x1xn and corresponding values fifxi, find the unique polynomial pnx of order n such that pnxifi for all i0, 1, . . . , n. A solution can be found in the form of the Lagrange interpolant . Define
lixn xxj i0, . . . , n. 10 j0 xixj
j i
Then, clearly, lixjij where ij denotes as usual the Kronecker delta, namely ij1 if
ij and ij0 if ij. Then the Lagrange form of the interpolating polynomial
n
pnx fj lj x, 11
j0
is clearly interpolatory as pnxinj1 fjljxifilixifi for all i.
In what follows we are going to assume that n1 and the interpolation problem is posed over the interval 1, 1 and that x01 and xn1. Interpolation on a generic interval a, b can be easily implemented by an affine mapping the problem. Some families of interpolation
3
nodes are particularly important. One of them is simply given by uniformly spaced nodes. Namely,
xj 1jh for 0jn, with h2n. 12 Another one is given by the Chebyshev nodes
xj cosj for 0jn. 13 n
These nodes are named after Chebyshev because the internal nodes all nodes apart from 1 are the extremal points of the Chebyshev polynomial Tn. These nodes are not uniform: they are more clustered towards the endpoints 1. This property makes them somehow more suitable to the interpolation problem than uniform spaced nodes, as we shall see. They are important because they constitute a class of GaussLobatto nodes2.
Given a set of nodes such as the uniform nodes 12 or the Chebyshev nodes 13, the inter polating polynomial can be evaluated as some x1, 1 by using 11 and 10.
a In a file lagrange.cpp write a C function lagrange which evaluates the Lagrange form of the interpolating polynomial of degree n at a given point x using 11 and 10. The function should input the number of nodes, the interpolation nodes xi and values fi, i0, 1, . . . , n, and the point x where the interpolating polynomial pn is to be evaluated and return the value pnx. Write a test program called q3.cpp designed to use the function lagrange with the uniformly spaced nodes 12 and with fjfxj with fxx4x. You may find it convenient to write a C function performing the evaluation of fx. Use it to perform the following tasks.
Set n2 and calculate the interpolating polynomial pnx at x1, x0, and x1 and verify that pn is interpolating by checking that the error fxpnx0 in all cases. Then, set x16 and output the error fxpnx.
For n4, the uniqueness of the interpolating polynomial implies p4xfx for all x. Check that this is the case by setting n4 and evaluating pnx at different locations. Report the result obtained when x16. 15
b Modify program q3.cpp to include, as a new set of nodes, the set of Chebyshev nodes given by 13. Use once again your function lagrange to evaluate the corresponding interpolating polynomial. Perform the following tasks.
Repeat the tasks in Question 3a to verify that the result is interpolating the given function. Report, as before, for n2 and x16, the error fxpnx.
Modify your program to use, instead, the function fx1116×2. Let again x16 and evaluate the interpolating polynomial obtained with both the uniform and Chebyshev nodes with n2 and n10. Report the resulting errors fxpnx in all cases. Consider now x.9, repeat the experiment and report, once again, the errors fxpnx. You should observe that as expected the interpolation error reduces as the number of nodes grows when the Chebyshev points are used. Instead, the error grows when the uniform nodes are used! This is a manifestation of the Runge phenomenon: equispaced nodes are not necessarily a good choice for global interpolation. Instead, the Chebyshev points can be proven to be the best choice. To visualise the phenomenon, write a test program q3b.cpp which, for n10 and for both the uniform nodes and Chebyshev nodes interpolant, evaluates at 100 uniformly distributed points zj1jk for 0jn with k2100 the value of the
2Gauss nodes and GaussLobatto nodes are defined as set of nodes that, together with a set of weights, define quadrature rules which are exact for all polynomials of degree 2n1, which is the highest possible order of exactness achievable.
4
interpolant. Output the resulting values without including them to your solution and use the output to plot the values of the two interpolants as well as the corresponding values of the function fzj using MATLAB or any other plotting tool of your liking. Save your plot as as the pdf file runge.pdf and include it within your solutions. 15
An alternative exists for the Chebyshev nodes which exploits the orthogonality property of the Chebyshev polynomials mentioned in Question 2. To this end, we look directly for the unique linear combination of the first n1 Chebyshev polynomials, that is in the form of the partial sum 9. The Chebyshev interpolation problem reads: find the coefficients ai of the polynomial
n
pnxaiPix such that pnxjfj. 14
i0
where xj are the Chebyshev nodes 13. This is the Chebyshev form of the interpolating
polynomial. The problem here is to evaluate the coefficients ai. In the special case in which xj cosj for 0jn, 15
n
note the change of sign with respect to 13: the nodes, in this case, are given as 1x0xn1 then there is a very efficient formula to do so. First of all, enforcing the interpolation equations and using 8 we immediately see that the interpolation condition can be written as: find the coefficients ai, i0,…,n such that
ai2n fjcosij i0,…,n, 16 ndi j0 dj n
where d0dn2 and dj1 otherwise. This is a type of discrete transform in that it yields the discrete coefficients in terms of the nodal values. It is known as the Chebyshev discrete transform CDT3.
c Use the above alternative way to evaluate the interpolating polynomial at the Chebyshev points given by 16. Within the file chebyshev.cpp, write a new function chebcoef that computes the coefficients of the Chebyshev form of the interpolating polynomial 14 given by 16. Modify the program q3.cpp to include, as a third set of nodes, the set of Chebyshev nodes xj, j0,1,…,n given by 15 and perform the following tasks. Evaluate the function fx1116×2 at the new set of Chebyshev nodes to obtain the interpolation values fjf xj . Call the function chebcoef to compute the corresponding Chebyshev coefficients ai. Report the coefficients in the case n2. Evaluate the resulting Chebyshev sum 14 at a point x using your function chebsum from Question 2. Report the resulting errors fxpnx when n2 and x.9. Compare your result with your findings from Question 3b. 15
October 2019
3The CDT can be efficiently computed using the Fast Fourier Transform FFT algorithm, although we are not going to take advantage of this here.
n ij
aicos n , j0,…,n. polynomials it can be shown that, indeed,
fj
This problem has a very elegant solution. Using the orthogonality property of the Chebyshev
i0
5
Reviews
There are no reviews yet.