Numerical Analysis

613 Pages·2013·15.901 MB·English

by  Sauer, Timothy

N
u
m
e
r
i
c
a
l

A
n
a
l
y
s
i
s

S
a
u
e
r

S
e
c
o
n
d
Numerical Analysis

E
d
Timothy Sauer
i
t
i
Second Edition
o
n
ISBN 978-1-29202-358-8
9 781292 023588
Numerical Analysis
Timothy Sauer
Second Edition
ISBN 10: 1-292-02358-9
ISBN 13: 978-1-292-02358-8
Pearson Education Limited
Edinburgh Gate
Harlow
Essex CM20 2JE
England and Associated Companies throughout the world
Visit us on the World Wide Web at: www.pearsoned.co.uk
© Pearson Education Limited 2014
All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted
in any form or by any means, electronic, mechanical, photocopying, recording or otherwise, without either the
prior written permission of the publisher or a licence permitting restricted copying in the United Kingdom
issued by the Copyright Licensing Agency Ltd, Saffron House, 6–10 Kirby Street, London EC1N 8TS.
All trademarks used herein are the property of their respective owners. The use of any trademark
in this text does not vest in the author or publisher any trademark ownership rights in such
trademarks, nor does the use of such trademarks imply any affi liation with or endorsement of this
book by such owners.
ISBN 10: 1-292-02358-9
ISBN 13: 978-1-292-02358-8
British Library Cataloguing-in-Publication Data
A catalogue record for this book is available from the British Library
Printed in the United States of America
11223344455273848473693618831841751541
P E A R S O N C U S T O M L I B R AR Y
Table of Contents
Chapter 0. Fundamentals
Timothy Sauer 1
Chapter 1. Solving Equations
Timothy Sauer 24
Chapter 2. Systems of Equations
Timothy Sauer 71
Chapter 3. Interpolation
Timothy Sauer 138
Chapter 4. Least Squares
Timothy Sauer 188
Chapter 5. Numerical Differentiation and Integration
Timothy Sauer 243
Chapter 6. Ordinary Differential Equations
Timothy Sauer 281
Chapter 7. Boundary Value Problems
Timothy Sauer 348
Chapter 8. Partial Differential Equations
Timothy Sauer 374
Chapter 9. Random Numbers and Applications
Timothy Sauer 431
Chapter 10. Trigonometric Interpolation and the FFT
Timothy Sauer 467
Chapter 11. Compression
Timothy Sauer 495
Chapter 12. Eigenvalues and Singular Values
Timothy Sauer 531
Answers to Selected Exercises
Timothy Sauer 565
I
569023
Bibliography
Timothy Sauer 592
Index 603
II
C H A P T E R
0
Fundamentals
This introductory chapter provides basic building Afterdiscussingefficientmethodsforevaluating
blocksnecessaryfortheconstructionandunderstand- polynomials,westudythebinarynumbersystem,the
ing of the algorithms of the book.They include fun- representationoffloatingpointnumbersandthecom-
damentalideasofintroductorycalculusandfunction mon protocols used for rounding.The effects of the
evaluation,thedetailsofmachinearithmeticasitiscar- smallroundingerrorsoncomputationsaremagnified
riedoutonmoderncomputers,anddiscussionofthe in ill-conditioned problems.The battle to limit these
lossofsignificantdigitsresultingfrompoorly-designed perniciouseffectsisarecurringthemethroughoutthe
calculations. restofthechapters.
Thegoalofthisbookistopresentanddiscussmethodsofsolvingmathematicalprob-
lemswithcomputers.Themostfundamentaloperationsofarithmeticareadditionand
multiplication. These are also the operations needed to evaluate a polynomial P(x) at a
particularvaluex.Itisnocoincidencethatpolynomialsarethebasicbuildingblocksfor
manycomputationaltechniqueswewillconstruct.
Because of this, it is important to know how to evaluate a polynomial. The reader
probably already knows how and may consider spending time on such an easy problem
slightlyridiculous!Butthemorebasicanoperationis,themorewestandtogainbydoingit
right.Thereforewewillthinkabouthowtoimplementpolynomialevaluationasefficiently
aspossible.
0.1
EVALUATINGAPOLYNOMIAL
Whatisthebestwaytoevaluate
P(x) 2x4 3x3 3x2 5x 1,
= + − + −
say, at x 1/2?Assume that the coefficients of the polynomial and the number 1/2 are
=
storedinmemory,andtrytominimizethenumberofadditionsandmultiplicationsrequired
From Numerical Analysis, Second Edition. Timothy Sauer. Copyright © 2012 by Pearson Education, Inc.
All rights reserved.
1
2 | CHAPTER0 Fundamentals
to get P(1/2). To simplify matters, we will not count time spent storing and fetching
numberstoandfrommemory.
METHOD1 Thefirstandmoststraightforwardapproachis
1 1 1 1 1 1 1 1 1 1 1 5
P 2 3 3 5 1 . (0.1)
2 = ∗ 2 ∗ 2 ∗ 2 ∗ 2 + ∗ 2 ∗ 2 ∗ 2 − ∗ 2 ∗ 2 + ∗ 2 − = 4
(cid:2) (cid:3)
Thenumberofmultiplicationsrequiredis10,togetherwith4additions.Twooftheadditions
areactuallysubtractions,butbecausesubtractioncanbeviewedasaddinganegativestored
number,wewillnotworryaboutthedifference.
There surely is a better way than (0.1). Effort is being duplicated—operations can
be saved by eliminating the repeated multiplication by the input 1/2.Abetter strategy is
tofirstcompute(1/2)4,storingpartialproductsaswego.Thatleadstothefollowingmethod:
METHOD2 Findthepowersoftheinputnumberx 1/2first,andstorethemforfutureuse:
=
1 1 1 2
2 ∗ 2 = 2
(cid:2) (cid:3)
1 2 1 1 3
2 ∗ 2 = 2
(cid:2) (cid:3) (cid:2) (cid:3)
1 3 1 1 4
.
2 ∗ 2 = 2
(cid:2) (cid:3) (cid:2) (cid:3)
Nowwecanadduptheterms:
1 1 4 1 3 1 2 1 5
P 2 3 3 5 1 .
2 = ∗ 2 + ∗ 2 − ∗ 2 + ∗ 2 − = 4
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
Therearenow3multiplicationsof1/2,alongwith4othermultiplications.Countingup,
wehavereducedto7multiplications,withthesame4additions.Isthereductionfrom14
to11operationsasignificantimprovement?Ifthereisonlyoneevaluationtobedone,then
probablynot.WhetherMethod1orMethod2isused,theanswerwillbeavailablebefore
youcanliftyourfingersfromthecomputerkeyboard.However,supposethepolynomial
needs to be evaluated at different inputs x several times per second.Then the difference
maybecrucialtogettingtheinformationwhenitisneeded.
Isthisthebestwecandoforadegree4polynomial?Itmaybehardtoimaginethat
we can eliminate three more operations, but we can. The best elementary method is the
followingone:
METHOD3 (NestedMultiplication)Rewritethepolynomialsothatitcanbeevaluatedfromtheinside
out:
P(x) 1 x(5 3x 3x2 2x3)
=− + − + +
1 x(5 x( 3 3x 2x2))
=− + + − + +
1 x(5 x( 3 x(3 2x)))
=− + + − + +
1 x (5 x ( 3 x (3 x 2))). (0.2)
=− + ∗ + ∗ − + ∗ + ∗
Herethepolynomialiswrittenbackwards,andpowersofx arefactoredoutoftherestof
thepolynomial.Onceyoucanseetowriteitthisway—nocomputationisrequiredtodo
therewriting—thecoefficientsareunchanged.Nowevaluatefromtheinsideout:
2
0.1 EvaluatingaPolynomial | 3
1
multiply 2, add 3 4
2 ∗ + →
1
multiply 4, add 3 1
2 ∗ − →−
1 9
multiply 1, add 5
2 ∗− + → 2
1 9 5
multiply , add 1 . (0.3)
2 ∗ 2 − → 4
Thismethod,callednestedmultiplicationorHorner’smethod,evaluatesthepolynomial
in 4 multiplications and 4 additions.Ageneral degree d polynomial can be evaluated in
d multiplications and d additions. Nested multiplication is closely related to synthetic
divisionofpolynomialarithmetic.
Theexampleofpolynomialevaluationischaracteristicoftheentiretopicofcomputa-
tionalmethodsforscientificcomputing.First,computersareveryfastatdoingverysimple
things.Second,itisimportanttodoevensimpletasksasefficientlyaspossible,sincethey
maybeexecutedmanytimes.Third,thebestwaymaynotbetheobviousway.Overthe
last half-century, the fields of numerical analysis and scientific computing, hand in hand
withcomputerhardwaretechnology,havedevelopedefficientsolutiontechniquestoattack
commonproblems.
While the standard form for a polynomial c c x c x2 c x3 c x4 can be
1 2 3 4 5
+ + + +
writteninnestedformas
c x(c x(c x(c x(c )))), (0.4)
1 2 3 4 5
+ + + +
someapplicationsrequireamoregeneralform.Inparticular,interpolationcalculationsin
Chapter3willrequiretheform
c (x r )(c (x r )(c (x r )(c (x r )(c )))), (0.5)
1 1 2 2 3 3 4 4 5
+ − + − + − + −
wherewecallr ,r ,r ,andr thebasepoints.Notethatsettingr r r r 0in
1 2 3 4 1 2 3 4
= = = =
(0.5)recoverstheoriginalnestedform(0.4).
The following Matlab code implements the general form of nested multiplication
(comparewith(0.3)):
%Program 0.1 Nested multiplication
%Evaluates polynomial from nested form using Horner’s Method
%Input: degree d of polynomial,
% array of d+1 coefficients c (constant term first),
% x-coordinate x at which to evaluate, and
% array of d base points b, if needed
%Output: value y of polynomial at x
function y=nest(d,c,x,b)
if nargin<4, b=zeros(d,1); end
y=c(d+1);
for i=d:-1:1
y = y.*(x-b(i))+c(i);
end
RunningthisMatlabfunctionisamatterofsubstitutingtheinputdata,whichconsist
of the degree, coefficients, evaluation points, and base points. For example, polynomial
(0.2)canbeevaluatedatx 1/2bytheMatlabcommand
=
3
4 | CHAPTER0 Fundamentals
>> nest(4,[-1 5 -3 3 2],1/2,[0 0 0 0])
ans =
1.2500
aswefoundearlierbyhand.Thefilenest.m,astherestoftheMatlabcodeshownin
this book, must be accessible from the Matlab path (or in the current directory) when
executingthecommand.
Ifthenestcommandistobeusedwithallbasepoints0asin(0.2),theabbreviated
form
>> nest(4,[-1 5 -3 3 2],1/2)
may be used with the same result. This is due to the nargin statement in nest.m.
If the number of input arguments is less than 4, the base points are automatically set to
zero.
BecauseofMatlab’sseamlesstreatmentofvectornotation,thenestcommandcan
evaluateanarrayofx valuesatonce.Thefollowingcodeisillustrative:
>> nest(4,[-1 5 -3 3 2],[-2 -1 0 1 2])
ans =
-15 -10 -1 6 53
Finally,thedegree3interpolatingpolynomial
1 1 1
P(x) 1 x (x 2) (x 3)
= + 2 + − 2 + − −2
(cid:2) (cid:2) (cid:2) (cid:3)(cid:3)(cid:3)
fromChapter3hasbasepointsr 0,r 2,r 3.Itcanbeevaluatedatx 1by
1 2 3
= = = =
>> nest(3,[1 1/2 1/2 -1/2],1,[0 2 3])
ans =
0
(cid:2)EXAMPLE0.1 FindanefficientmethodforevaluatingthepolynomialP(x) 4x5 7x8 3x11 2x14.
= + − +
Some rewriting of the polynomial may help reduce the computational effort
required for evaluation. The idea is to factor x5 from each term and write as a polyno-
mialinthequantityx3:
P(x) x5(4 7x3 3x6 2x9)
= + − +
x5 (4 x3 (7 x3 ( 3 x3 (2)))).
= ∗ + ∗ + ∗ − + ∗
For each input x, we need to calculate x x x2, x x2 x3, and x2 x3 x5 first.
Thesethreemultiplications,combinedwith∗the=multiplic∗ation=ofx5,andthe∗thre=emultipli-
cationsandthreeadditionsfromthedegree3polynomialinthequantityx3 givethetotal
operationcountof7multipliesand3addsperevaluation. (cid:3)
4
0.2 BinaryNumbers | 5
0.1Exercises
1. Rewritethefollowingpolynomialsinnestedform.Evaluatewithandwithoutnestedformat
x 1/3.
=
(a) P(x) 6x4 x3 5x2 x 1
= + + + +
(b) P(x) 3x4 4x3 5x2 5x 1
=− + + − +
(c) P(x) 2x4 x3 x2 1
= + − +
2. Rewritethefollowingpolynomialsinnestedformandevaluateatx 1/2:
=−
(a) P(x) 6x3 2x2 3x 7
= − − +
(b) P(x) 8x5 x4 3x3 x2 3x 1
= − − + − +
(c) P(x) 4x6 2x4 2x 4
= − − +
3. EvaluateP(x) x6 4x4 2x2 1atx 1/2byconsideringP(x)asapolynomialinx2
= − + + =
andusingnestedmultiplication.
4. EvaluatethenestedpolynomialwithbasepointsP(x) 1 x(1/2 (x 2)(1/2 (x 3)
= + + − + −
( 1/2)))at(a)x 5and(b)x 1.
− = =−
5. EvaluatethenestedpolynomialwithbasepointsP(x) 4 x(4 (x 1)(1 (x 2)
= + + − + −
(3 (x 3)(2))))at(a)x 1/2and(b)x 1/2.
+ − = =−
6. Explainhowtoevaluatethepolynomialforagiveninputx,usingasfewoperationsas
possible.Howmanymultiplicationsandhowmanyadditionsarerequired?
(a)P(x) a a x5 a x10 a x15
0 5 10 15
= + + +
(b)P(x) a x7 a x12 a x17 a x22 a x27.
7 12 17 22 27
= + + + +
7. Howmanyadditionsandmultiplicationsarerequiredtoevaluateadegreenpolynomialwith
basepoints,usingthegeneralnestedmultiplicationalgorithm?
0.1ComputerProblems
1. UsethefunctionnesttoevaluateP(x) 1 x x50atx 1.00001.(Usethe
= + +···+ =
Matlabonescommandtosavetyping.)Findtheerrorofthecomputationbycomparingwith
theequivalentexpressionQ(x) (x51 1)/(x 1).
= − −
2. Usenest.mtoevaluateP(x) 1 x x2 x3 x98 x99atx 1.00001.Finda
= − + − +···+ − =
simpler,equivalentexpression,anduseittoestimatetheerrorofthenestedmultiplication.
0.2
BINARYNUMBERS
In preparation for the detailed study of computer arithmetic in the next section, we need
tounderstandthebinarynumbersystem.Decimalnumbersareconvertedfrombase10to
base2inordertostorenumbersonacomputerandtosimplifycomputeroperationslike
additionandmultiplication.Togiveoutputindecimalnotation,theprocessisreversed.In
thissection,wediscusswaystoconvertbetweendecimalandbinarynumbers.
Binarynumbersareexpressedas
...b b b .b b ...,
2 1 0 1 2
− −
5