### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# Numerical Analysis MATH 128A

GPA 3.93

### View Full Document

## 19

## 0

## Popular in Course

## Popular in Mathematics (M)

This 16 page Class Notes was uploaded by Kavon Feest on Thursday October 22, 2015. The Class Notes belongs to MATH 128A at University of California - Berkeley taught by Staff in Fall. Since its upload, it has received 19 views. For similar materials see /class/226600/math-128a-university-of-california-berkeley in Mathematics (M) at University of California - Berkeley.

## Reviews for Numerical Analysis

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/22/15

Math 128 B Iterate April 7 2004 1043 am Conjugate Gradients Positive De nite A A39 To solve AZ b for z Ailb choose x71 x0 arbitrarily and for n 0 1 2 3 in turn compute Xn1 xn Snb 7AxIQ anxn 7 x1171 with scalars Sn and um chosen to minimize xn17z39Axn17z n139Axn1 7 2b39xn z39Az Simplify typography by dropping subscripts from rn b 7 Axn dn xn 7 x1171 Sn and um to nd that rd39Ard PAY Md 5 7 fr rd39r u d Ar d Ad u d39r Solve this equation for B and u At least one solution always exists if A is positive de nite In the absence of roundoff successive residuals rn b 7 Axn turn out to be orthogonal in fact rm39rn 0 dm39Adn for all m gt n 2 0 Consequently rm o and xm 2 at least as soon as m equals the dimension of A But the point of the iteration is not to iterate that many times when the dimension is huge Instead take advantage of the tendency of residuals rn and increments d to dwindle as n increases and stop iterating when they both become small enough Overrelaxation Positive De nite A A39 7L V 7 L39 in which V DiagA V39 is Positive De nite too and 7L SubdiagA Given A and b we seek z Ailb Starting from an arbitrary initial guess x0 for n 0 1 2 3 in turn ordinary GaussSeidel iteration solves Vxn1 7x1 b anH 7 Vxn L39xn for Xn1 V 7 L 1 b L39xn Then Xn1 7z Exn7z where E V 7 L 1L39 can be shown to have eigenvalues all with magnitudes less than 1 though not necessarily much less unless HLH is rather smaller than the smallest eigenvalue of V To accelerate convergence consider using an OverUnderrelaxation parameter 5 con ned to 71 lt 5 lt 1 for reasons to be explained later To solve AZ b for z Ailb choose x0 arbitrarily and for n 0 1 2 3 in turn solve Vxn17 n 15b anH 7Vxn L39xIQ for Xn1 V 7 15L 1Vxn 15b 7 Vxn L39xIQ Then xn 7z Exn7z where E V 7 15L 143V 15L39 Note that because L lies strictly below the diagonal product of all E s eigenvalues detE det l Therefore at least one eigenvalue of E has magnitude at least as big as 151 This is why we keep 71 lt 5 lt 1 E and its eigenvalues depend upon 5 as well as A but the dependence is obscure except in special cases An important special case arises when A 1B f In this case every eigenvalue 0c of A has the form 0L 1 i B where B is either asingular value of B or if B is not square 0 And in this case every eigenvalue 8 of E can be shown easily to satisfy 85 iB15E e 715152 4 15B22 7 5 2 The largest ofthe quot J of 39 39 e is 39 39 39 J when 5 HBH1 V1711B1122 and then every eigenvalue 8 has the same magnitude 5 Prof W Kahan Page 11 Filename Ellipsoi Ellipsoidal ErrorBounds for Trajectory Calculations November 14 1999 1 1 12 pm Abbreviated Lecture Notes on Ellipsoidal Error Bounds for Trajectory Calculations Prof W Kahan Mathematics Dept and Elect Eng amp Computer Science Dept University of California at Berkeley For Presentation Oct 12 1993 Abstract A practical way is outlined to bound error accrued during numerical calculations of traj ectories The algorithm accommodates uncertainties in the governing differential equation as well as error due to the numerical process The accrued error is constrained to lie within an ellipsoid that can be proved to grow as time 17 increases to 00 bigger than the worst possible accrual by a factor no worse than 1 NE for some constant B rather than exponentially bigger until nonlinearity in the differential equation forces a singularity to manifest itself Introduction A trajectory is the graph of y l7 the solution of an Autonomous Initial Value Problem dyd c y39 c fy c for E 2 0 y0 yo AIVP Here yo is a given vector and f a given vectorvalued function of its vector argument Perturbing yo to yo Syo and f to f of perturbs the trajectory to the solution y 5y of a perturbed AIVP y SYYW f 50W 5301 for T 2 0 y 5300 yo 5yo Given only some kind of bound for Syo and 5f how can we infer a bound upon 5y Actually in practice we know y 5y but not y this will be overlooked now to simplify the exposition To the same end we shall not discuss how the truncation and rounding errors incurred by the numerical process that solves the AIVP can be incorporated into 5f along with errors due to idealizations that model a complex physical situation by a simpli ed expression f Were all perturbations 5 in nitesimal the accrued error 5y would satisfy the Aaljoint or Variational Initial Value Problem associated with the given AIVP 5y3917 J175y c of 5y17 Syo VIVP Here Jacobian matrix J l7 f y BfyBy at y y l7 this J l7 is the matrix of rst partial derivatives of f evaluated at y l7 the presumedtobeknown trajectory Perturbation 5f 5fy l7 however whereas J l7 can be computed from f and y l7 all we know about of is an upper bound Likewise for Syo Presumably computable perhaps functions of y l7 these bounds upon of and Syo are the data from which we wish to infer a bound upon 5y l7 Computing this bound for accrued error along with y l7 is the problem addressed in these notes In practice the perturbations 5 are not in nitesimal The effect in practice of their finiteness is to increase the bound upon of by an amount roughly proportional to the product of the square of the computed bound upon 5y l7 and a bound upon the second derivative of f This increase turns the linear VIVP into something nonlinear like a Riccati equation whose solution may Prof W Kahan Page 15 Filename Ellipsoi Ellipsoidal ErrorBounds for Trajectory Calculations November 14 1999 l l 12 pm become in nite at a nite time 17 even in some cases when 5y l7 is known to stay bounded for all nite 17 Despite its importance this nonlinear contribution is not discussed in these notes Instead to simplify the exposition of a subject that is already complicated enough 5y is assumed to stay so small that its secondorder contributions may safely be neglected We assume J l7 and of to be continuous for all 17 in practice of may be only piecewise continuous This is a technicality we could dispatch by converting differential equations into equivalent Volterra integral equations For simplicity s sake we won t do that either Despite all our simpli cations the computation of a bound upon 5y remains so challenging that every previously published scheme I know about is prone to producing bounds too big by a factor that grows like an exponential function of 17 when J 17 behaves in a way that the scheme dislikes The algorithm described here is far less pessimistic its bound cannot grow too big by a factor bigger than 1 SWE for some constant B that depends upon moduli of continuity of J 17 and of the bound for of The proof of this claim is too long to t here Space and time barer suf ce for an outline of my algorithm My algorithm adjoins to the given AIVP another differential equation whose solution intended to be computed simultaneously with y l7 is a symmetric matrix A E that describes an ellipsoid centered at y l7 and surely big enough to enclose 5y l7 but not too much bigger The adjoined differential equation s dimension is about half the square of y s dimension so A E may well cost enormously more to compute than y l7 Since previously published schemes typically cost twice as much as mine I need not apologize for it The Reachable Set Let us now write z in place of 5y and u in place of of in the VIVP above it becomes 2Y1 J E z c um for E 2 0 z0 zo VIVP The matrix J l7 is assumed computable but for u l7 and zo only bounds are available We construe these bounds as constraints that restrict u l7 and zo to certain small regions say u l7 e Uh and zo e A for given centrally symmetric convex bodies TONI and A characterized by parameters to be discussed later For instance these bodies could be spheres characterized by their radii which are then upper bounds for the lengths of u l7 and zo Parallelepipeds have been used too characterized by the matrices that map a unit hypercube onto them But we use ellipsoids for reasons to be discussed later The Reachable Set 17 consists of all values that z l7 can take compatible with the given constraints imposed by TDY E and A upon u and zo in the VIVP Reachable Set is a term coined by Control theorists From the given hypotheses about lift and A it follows that T must be a centrally symmetric convex body too But generally the shape of 51 is not so simple as the shapes of 17 and A Regardless of whether the latter are ellipsoids or parallelepipeds 17 need not be any of those The best we can expect to do computationally is to approximate 51 by one of those simpler gures Thus our task is to compute whatever parameters characterize a simpler centrally symmetric convex body ellipsoid or parallelepiped that circumscribes reachable set T as tightly as is possible at a tolerable cost Prof W Kahan Page 25 Filename Ellipsoi Ellipsoidal ErrorBounds for Trajectory Calculations November 14 1999 1 1 12 pm Ellipsoidal Bounds Any open ellipsoid A centered at the origin 0 is characterized by an appropriate symmetric positive de nite matrix A as follows A consists of all vectors x that satisfy xTAilx lt 1 The eigenvectors of A point in the directions of principal semiaxes of A and the eigenvalues of A are the squared lengths of those semiaxes Flattened ellipsoids those with some semiaxes of zero length are represented by singular positive semide nite matrices A for which we must understand the expression xTAilx to require that x be con ned to the range of A A theorem of Fritz John 1948 asserts that any centrally symmetric convex body 13 in an N dimensional space can be circumscribed by an ellipsoid Y tightly enough that the boundary of lies inside Y but outside YWN A short proof of his theorem is posted in lecture notes at lthttp www cs berkeley eduwkahanMathHllONORMlite pdfgt His chosen candidate Y is the circumscribing ellipsoid of minimum volume its characterizing matrix Y minimizes detY Thus little can be lost in spaces of modest dimension if the bounding bodies 101 and A containing respectively u l7 and zo above are taken to be ellipsoids otherwise they can be circumscribed by ellipsoids at the cost of worsening our bound Aft upon the reachable set 51 by at worst a factor WN As error bounds go this degree of exacerbated pessimism will go unremarked Thus we may assume that along with the Jacobian matrix J 17 of partial derivatives we are supplied with symmetric positive semide nite matrices U E and A0 whose ellipsoids 1611 and A we are told contain u l7 and zo Often U E and A0 will be diagonal Our task is to compute a symmetric positive de nite matrix A l7 whose ellipsoid AFC circumscribes the reachable set T as tightly as possible at a tolerable cost The Auxiliary Differential Equation First let s summarize the situation as it stands now The given AIVP y39 l7 fy l7 for 17 2 0 and y0 yo AIVP is being solved numerically for y but unknown perturbations 5y0 zo and 5f u l7 induce in y an accrued perturbation 5y z that satis es 2Y1 J Ez c u l7 for 17 2 0 and z0 zo VIVP Here J l7 f y l7 is known and so are symmetric positive de nite matrices U E and A0 A0 for which uTUilu lt1 for 17 20 and zOTAoilzoltl at 17 0 By constraining u and zo these inequalities compel the accrual z to lie in a Reachable Set T about which we wish to circumscribe an ellipsoid AFC by computing its symmetric positive de nite matrix A E such that zTAilz lt 1 for every z in for every 17 2 0 Here all matrices except A0 are functions of 17 In nitely many matrices A ful ll these requirements we seek one of the smaller ones Prof W Kahan Page 35 Filename Ellipsoi Ellipsoidal ErrorBounds for Trajectory Calculations November 14 1999 1 1 12 pm Here is how to construct one First let 0L be any strictly positive piecewise continuous scalar function of 17 later we shall exhibit a good choice for 0L Then compute A E as the solution of the following Auxiliary Di erential Equation A39 J39AA39JTOL39UAOL for 120 and A0A0 ADE The solution A l7 of this ADE is intended to be computed numerically and simultaneously with the solution y l7 of the AIVP Theorem zTAilz lt 1 for every 2 in the reachable set 5 Proof of the Theorem The theorem s inequality starts true at 17 0 so if increasing 1 ever makes it false it must do so for the rst time at some 17 gt 0 at which both zTAilz 1 and zTAile 2 0 To show that the last inequality is incompatible with the previous equality we need only expand zTAile zTA lzy Jz uTA 1z 7 zTA laA AJT ocU NooA lz zTA laz u uTAilz zTAilu 7 OLZTATIUATIZ 7 ZTATIZOL OLA IZ 7 U luTU0LA lz 7 U lu0L 7 1 7 uTU lu0L lt 0 as claimed so the theorem is true for all 17 gt 0 How to Choose 0L Having just proved that ellipsoid A encloses reachable set 6 we are dismayed to observe from the ADE that A may well enclose a great deal more than if 0c is either too big or too small Can 0L be so chosen that A is not much bigger than In fact there is a way albeit impractical so to choose 0L that A and 6 will share a common support plane tangent for all 17 2 0 provided the normal to that plane is fixed in advance Even so A may still exceed the size of enormously in directions parallel to that plane Apparently choosing 0L well is a subtle problem All the more surprising then is the existence of a computationally simple choice that always turns out to be adequate Choose 0L traceA traceU For this choice and some others I have proved that the diameter of A cannot exceed that of by a factor bigger than 1 SWE where B is a constant that depends upon various attributes of J and U My proofis still so long thatI am too embarrassed to publish it Computational Experience My earliest experiments with ellipsoidal bounds are still the most satisfying The AIVP was the equations of motion of one of the moons of Jupiter and the ADE included crude bounds for the gravitational in uences of the rest of the solar system plus the contributions of all numerical errors The computation ran for hundreds of orbits during which the bounding ellipsoid A became ever more needleshaped growing roughly like 1 E2 Computation was halted only because I wished not to hog the University of Toronto39s IBM 7094 in 1968 The results con rmed that in a situation where cof nshaped or more general parallepiped shaped bounds became in nite after a few orbits ellipsoidal bounds did not Prof W Kahan Page 45 Filename Ellipsoi Ellipsoidal ErrorBounds for Trajectory Calculations November 14 1999 1 1 12 pm Acknowledgements I am indebted to Gerry Gabel who programmed the differential equation solver that I most preferred on the 7094 I have yet to see a better program it was a predictorcorrector with automatically chosen continually varying order and stepsize The late Prof Tom Hull endured my many attempts to explain why ellipsoidal bounds could succeed where all shapes promoted previously were doomed to exponentially excessive growth and I am grateful for his patience And I am grateful to the late Prof Fred C Schweppe of MIT for his encouragement and for including some of the foregoing material in his book Uncertain Dynamic Systems 1973 PrenticeHall NJ where it has languished through no fault of his These notes were extracted from voluminous classroom handouts for a rarely offered graduate course Math 273 at the Univ of Calif at Berkeley at the instigation of Prof Jerrold Marsden for presentation at his workshop on Integration Algorithms for Classical Mechanics 1417 Oct 1993 held at the Fields Inst for Research in Math Sciences then at the Univ of Waterloo Ontario Canada Prof W Kahan Page 55 Re Rotn March 13 2006 743 am Re ections Rotations and QR Factorization QR Factorization gures in Least Squares problems and Singular Value Decompositions among other things numerical These notes explain some re ections and rotations that do it and oiTer MATLAB implementations in its notation x39 complex conjugate transpose of x Householder Re ections A HouseholderRe ection is W I 7 ww39 W39 W 1 for any column w satisfying w39w 2 If yWx then y39y x39x and y39xx39y is real even if x y and w are complex W is a re ection because Ww 7w and Wp p whenever w39p 0 Numerical analysts name this re ection after Alston S Householder because he introduced it to them in the mid 1950s as part of an improved way to solve Least Squares problems Given columns x and e 1 0 0 039 so that e39x x1 we seek column w so that w39w2 and W I7ww39 re ects x to WxeB for some scalar 5 Its lBl HxH Vx39x and x39Wxx39eBx139B must be real so BinHx1lx1l if x1 0 Construction Set N HXH 15 iNxllxll d x 7 eB and w dd39d2 unless d o in which case set w o But all bets are off if UNDERFLOW degrades x39x Proof Let p x eB so that p39d 0 p39w Then Wd 7d and Wp p whereupon 2Wx Wpd p7d 2e 3915 as desired How is the sign i in B chosen The simplest way maximizes 22 d39d2 NN 7 ilx1l by setting 5 7Nx1lX1l as we ll see Ofcourse any i sign works when x1 0 Detailed Construction Let V x 7 exl so that e39V V1 0 and let u V39V gt 0 so that N Vx39x W M lx1l2 Next set 9 xllxll signx1 except that we reset 91 if x1 0 Next we choose 5 icN Numerical stability requires two cases to be distinguished If B7cN set d x7eBxecN by copying x to d and then resetting d1 x1 7159lx1l N If BcN set d x7eBx7ecN by copying x to d and then resetting d1 x1 B 39Hlxll N Next 2 ld1l2 u2 ldllN and w dQ Return w B hshldrwx Prof W Kahan Math 128B Page 16 Re Rotn March 13 2006 743 am QR Factorization Given an mbyn matrix F with no fewer rows than columns so m 2 n we wish to factorize F QR with Q39Q I and R uppertriangular by using Householder re ections thus Wn3939W239W139F IS in which each re ection WJ Wj39 WJT1 is constructed to annihilate all subdiagonal elements in column j of Fj1 Wj1W2W1F Then Q W1W2WnH Each WJ I 7 wjwj39 has wJ wj 2 or 0 and no nonzero element in wJ above row j Each FJ WjFJq has the same rst jil rows as Fj1 and no nonzero subdiagonal elements in its rst j columns Each QJ WjWJH Wn has ones on the diagonal and zeros elsewhere in its rst jil rows and columns so QJ WJ Qj1 is obtained by altering only rows and columns of QJH with indices no less than j Detailed Construction in MATLAB Start with F0 F For j l 2 n in turn get w 5 hshldrwFJjmj as above and store E in place of FJQ39mj then if j lt n overwrite FJjmjln iwjwj39Fjjmjln onto FJjmjln to get Fj1jmjln Next R DiagBl 52 511 triuFnln ln 1 Finally set Gn FH and for j n nil l in turn extract Wj from Gj1jmj overwrite column OJ1 l omj onto Gj1 j and then onto Gjljm jn overwrite Gjjmjn Gjljmjn 7vVJvVJ39Gj1jmjn Then Q G1 Return Q R hshldrqrF Numerical experiments indicate that MATLAB uses the same method to get Q R qrF 0 QR Factorization by Givens Rotations A Givens Rotation is Q CY 5 so chosen that a 2vector v is rotated to Qv is c y wherein lrl2 v39v so c2 s39s 1 when by convention we choose c 2 0 Here v39 is the complex conjugate transpose of v and s39 is the complex conjugate of s The rotation is named after Wallace Givens who introduced this rotation to numerical analysts in the 1950s while he was working at Argonne National Labs near Chicago The rotation is encoded in one complex number t yx from which are derived c ll t39t s ct and r xc In the special case that t oo presumably because x 0 we set c 0 s l and r y In any event note that Q 1 Q39 Return c s t r givenstx y Prof W Kahan Math 128B Page 26 ReflRotn March 13 2006 743 am Bottom Up QR Factorization Given an mbyn matrix F with no fewer rows than columns so m 2 n we wish to factorize F QR with Q39Q I and R uppertriangular by using Givens rotations thus For 1 Si Smil and l S j S n let Qij be the Givens rotation that acts upon an mbyn matrix Z Z to overwrite Qij ZN fiij onto ZN We shall premultiply F by asequence of 0 z i1 i1 rotations Qij in this order from right to left for j 1 up to n in turn for i mil down to j in turn premultiply by Qij Since each Qij affects only rows i and il of columns j to n of the product we may store ti in place of the product s zero element in position il j since it will not gure in subsequent premultiplications After the last premultiplication we nd R in the product s rst n rows and columns after ignoring the subdiagonal elements that hold tils Then these are used to construct Q as a product of inverse rotations Qij39 premultiplying m in this reverse order from right to 0 left for j n down to l in turn for i j up to mil in turn premultiply by Q3 Each premultiplication by Qij39 affects only rows i and il of columns j to n of the product after tij was extracted from location il j and replaced by 0 Return Q R gvnsupqrF This is not the only way to use Givens rotations for QR factorizations Another is Top Down QR Factorization Given an mbyn matrix F with no fewer rows than columns so m 2 n we wish to factorize F QR with Q39Q I and R uppertriangular by using Givens rotations thus For 1 S j S n and jl Si S m let Qij be the Givens rotation that acts upon an mbyn matrix Z to overwrite QijZjvjj Ziyj r161 onto We shall premultiply F byasequence of rotations Qij in this order from right to left H for j 1 up to n in turn for i jl up to m in turn premultiply by Qij Since each Qij affects only rows i and j of columns j to n of the product we may store tij in place of the product s zero element in position i j since it will not gure in subsequent Prof W Kahan Math 128B Page 36 Re Rotn March 13 2006 743 am premultiplications After the last premultiplication we nd R in the product s rst n rows and columns after ignoring the subdiagonal elements that hold tijs Then these are used to construct Q as a product of inverse rotations Qij39 premultiplying in this reverse order from right to left for j n down to l in turn for i m down to jl in turn premultiply by Q3 Each premultiplication by Qij39 affects only rows i and j of columns j to n of the product after tij was extracted from location i j and replaced by 0 Return Q R gvnsdnqrF MATLAB appears to use Householder re ections to get its Q R qrF 0 There is no reason to expect any two of the three dilTerent Q R factorizations to agree though they must be related in the absence of roundolT R1 R2 1 Q139Q2 must be a diagonal unitary matrix MATLAB Programs function F R hshldrqrF Q R hshldrqrF uses Householder Reflections to factorize F QR so that R is upper triangular and do do Q has orthonormal columns Q Q I This works only if F has no more columns than rows and if underflow does not degrade F F Uses hshldrwm m n sizeF if m lt n error F has more columns than rows in hshldrqrF end Z zerosl n w zerosm l for j lzn w Zj hshlerFjm j Fjm j W i Fjm jln ww Fjm jln 1 1 lzn lzn l i w Fjm j F j zerosml Fj j l Fjm jzn Fjm jzn WW Fjm jzn i 96 39 Prof W Kahan Math 128B Page 46 Re Rotn March 13 2006 743 am function w 2 hshldrwX w 2 hshldrwX yields w with w w 2 or O so W I ww W WA l reflects the given column X to WX 2 O O O with Z normX But all bets are off if UNDERFLOW degrades X X w X m lengthw X1 wl a1 absXl if m lt 2 w O Z Xl return end if a1 s Xlal else s l end vv w2 m w2m aX sgrtalal vv Z saX a1 a1 aX wl sal dd2 alaX if dd2 w wsgrtdd2 end so w w 2 unless w o function c s t r givenstX y c s t r givenstX y satisfies c gt 0 cA2 sA2 l r cX sy t Xy s c So c s X r and c lsgrtl t t s c y 0 s ct r Xc if X O t conjyX u sgrtl t t r uX c lu s ct else t inf r y c O s 1 end 95 function F R gvnsupgrF Q R gvnsupgrF uses Givens Rotations to factorize F QR so that R is upper triangular and Q has orthonormal columns Q Q I This works only if F has no more columns than rows and if underflow does not degrade F F Uses givenstm bottom up m n sizeF if m lt n error F has more columns than rows in gvnsupgrF end for j lzn for i mlzlzj c s Fil j Fi j givenstFi j Fil j if j lt 11 Fiil jln c s s cFiil jln end end end i mlzlzj j lzn R triuFln ln for j n 11 Flzj j zerosjl i FjIj 1 i for i jm l t Fil j Fil j O c lsgrtl t t if c0 s ct else s l end Fiil jn c s s cFiil jn end end i jm l j nz lzl Prof W Kahan Math 128B Page 56 Re Rotn March 13 2006 743 am function F R gvnsdnqrF Q R gvnsdnqrF uses Givens Rotations to factorize F QR so that R is upper triangular and Q has orthonormal columns Q Q I This works only if F has no more columns than rows and if underflow does not degrade F F Uses givenstm top down m n sizeF if m lt n error F has more columns than rows in gvnsdnqrF end for j lzn for i jlm 0 S Fi j Fj j givenSHFU j Fi j i if j lt 11 Fji jln c s s cFji jln end end end i jlm j R triuFln ln for j nzl 1 Flzj j zerosjl i FjIj 1 i for i m l jl t Fi j Fi j O c lsqrtl t t if c s ct else s l end Fji J 1391 C s 5 cFji jzn i end end i m ljl j nz lzl Prof W Kahan Math 128B Page 66 Filename SingNewt February 5 2007 750 am Notes on Nonlinear Newton Iterations at a Typical Singularity Given anonlinear nvectorvalued function fX of an nvector X a solution z of the equation fz o is often sought by means of Newton 3 Iteration Xk1 N fXk wherein Newton 3 Iterating Function NfX X if39X7139fX Here f39X BfXBX is the nbyn Jacobian Matrix of rst partial derivatives it is usually presumed adequately diiTerentiable and invertible at all X in some open neighborhood of the desired solution z and then N fX 7 z f39X 1 fquotXX7z X7X z OOH2 as x gt z Those presumptions imply at least Quadratic Convergence This happens almost always The iteration s behavior is not so easy to characterize in some singular situations Among these the most common has detf39X 0 but 3 detf39XBX i 0 at X z Thus detf39X 0 along some curve r1 2 or hypersurface r1 2 3 that passes through z but does not intersect itself there not even tangentially Call this locus It divides some open neighborhood of z into two open regions in which detf39X takes opposite signs N fX is either ambiguous or in nite on because f 39X is not invertible there These notes eXplore the behavior of Newton s iteration when it converges to z but no iterate falls into Convergence turns out to be Linear rather than quadratic and iterates approaching z usually tend to avoid as we shall see Our conclusions are summarized on p 3 and illustrated by eXamples on p 4 The Derivative row h x Bdetf39xBx Jacobi 3 Formula for the derivative of a determinant says d detB TraceAdjBdB wherein the Trace is the sum of all diagonal elements and AdjB is the ClassicalAaljoint orAaljugate AdjB detBB71 when detB 0 and is otherwise de ned by the continuity of what turns out to be a polynomial function of the elements of B de ned by BAdj B E AdjBB E detBl in general There are other equivalent de nitions of AdjB in terms of determinants or the Characteristic Polynomial of B but all we need from them is this fact Adj B i 0 just when the nbyn matriX B has RankB 2 n71 and then RankAdjB n 7 n71n 7 RankB Jacobi s formula is derived at ltwww csberkeley eduwkahanMathHllOjacobi pdfgt Jacobi s formula says d det f 39X TraceAdj f 39X f quotX dX wherein the second derivative fquot is a bilinear operator that maps nvectors y and z each linearly to an nvector f quotX yz that in general varies nonlinearly with X and if continuously fquotXyz fquotXzy In particular f quotX 39dX is a matriX whose every element depends linearly upon column dX so there must be some row nvector h X B det f 39X BX satisfying d det f 39X h XdX for all dX and if h z i 0 then the equation of the plane tangent to at z is h Z39XZ 0 We assume henceforth that h z i o and detf39z 0 whence Rank f 39z n71 and hence RankAdjf39z l so Adj f 39z wv for two nonzero vectors satisfying f39zw o and v f39z o Consequently h zdX TraceAdjf39zfquotzdX Tracewv fquotzdX Tracev fquotzwdX and therefore h z v f quot zw B det f 39XBX at X z We shall need this formula later Prof W Kahan Math 128B Page 14 Filename SingNeWt February 5 2007 750 am A Taylor Series expansion of f will be presumed valid for x in some open neighborhood of z fx o f39zx7z fquotzx7zx7z t13f quotzx7zx7zx7z f39X 1quot2 fquotZ39XZ t f39quotZ39XZ39XZ t detf39Z 0 Adjf39Z W39V detf39x h x7z wherein h v fquotzw i o Needed next is a Taylor series expansion for Newton s iterating function NfX x 7 f39X 139fX z 0 fquotx7z f 39x7zx7z 71fquotx7z f 39x7zx7z x7z in which all derivatives of fx are evaluated at x z This is too messy It must be simpli ed A Change of Coordinates By a process akin to Gaussian Elimination with Pivotal Exchanges of both rows and columns we can obtain a diagonal matrix L 1 f 39zU 1 M I 3 using suitably permuted lower and 0 uppertriangular matrices L and U These gure in a simplifying change of variables from x to u Ux7z Let gu LT1 fz U71 u so that Newton s iterating function for g is Ngu u 7 g39u 1 gu UNfz U71 u 7 z In other words the iteration Xk1 N fxk starting from x0 is mimicked by the iteration uk1 Nguk starting from uo Uxo 7 z in so far as every uk Uxk 7 z Iterates xk gt 2 just as fast if at all as uk gt o Thus no generality is lost by assuming 2 o fo f39o M I 3 Adjf39o 0 vv 039 039 in which w v o l and the foregoing Taylor series for Nfx is simpli ed to L 1 1 1 2 Nfx7 iMfquot x if xx fquot x 5f39quotxx L 1 1 1 71 1 7 x7 Mfquotx f39quotxx M76f39quotxx x in which all derivatives of fx are evaluated at x o Further progress requires a partition of fquot x XE in which Xx is an n7lbyn7l matrix whose every element is a linear X39 39 quotX function of x Similarly for the column Cx the row x R and the scalar h x In fact h x v fquotxv v fquotvx v R h vx so h v R h v Bdetf39xBx at x o as expected and h i 0 has been assumed Later much of our analysis will be affected by whether the last element of h namely h v is zero Let LW We f x 1 M fquotx f39mxx 4r 11xX OX f39mxx 1 r39 u X39R h39X Here a process akin to Gaussian elimination provides estimates gx 7 detI Xx 0x2detf39x 7 1 h x Ox2 column cx 7I7XxCx 0x2 row r x 7x R I7Xx 0x2 and matrix Wx I7Xx Mxcxr x 0x2 as x gt o Prof W Kahan Math 128B Page 24 Filename SingNeWt February 5 2007 750 am These estimates produce an estimate of NfX 1 W 1 39X MX39OX3 7H r39 What happens neXt depends upon whether X approaches z o too nearly tangentially to the surface on which detf39X 0 If so hquotX 0X2 and then pX lOX2 can be arbitrarily big so big that the computation of f39X 1fX and thus N fX malfunctions because of roundolT if not division by zero To avoid that malfunction we must keep lh gtgt 0X2 though h 39X OX so that terms like MX39OX2 stay no bigger than 0X as X gt o Under these circumstances WX I 0X and then N fX gainr l39Xv 002 OX39v 0x2 Thus if N fX is not 002 it is likely to take the form N fX Bv 032 for some tiny scalar B OX Whether these circumstances persist and avoid malfunctions depends upon whether the last element h v of h is zero If as is most likely h v i 0 then iterates of the form X Bv 032 for suf ciently tiny nonzero scalars 5 turn into NfX Bv 032 rs X because pX lBh v 0B2 is not too tiny and then convergence is linear with rate log2 and iterates approach z along a direction V that is Transverse not tangential to the surface In the unlikely case that h v 0 the iteration s behavior is di icult to predict because although iterates X tend often to come close to the form X Bv 032 it is too nearly tangential to Summary If fz o and detf39z 0 and Newton s iteration is started close enough to z but not too much closer to the locus on which detf39X 0 the iteration s convergence to z is linear with rate log2 provided a technical condition B det f 39X BX at X z Adj f39z i o is satis ed by the second derivative f quotz as is usually the case This technical condition can be described in terms of nonzero nullvectors v i o v f39z and w i o f39zw of singular matriX f39z it is that v f quotzww 0 And then iterates Xk1 Xk 7 f39Xk 1fXk approach the desired zero z very nearly like z BZkw for some small nonzero scalar constant B Example p For column 2vector arguments X let column 2vector pX a BX M1 2 in which x39C2x a 11512B1 2C1 jC2z Now p 172 o and p39 32 1342 is 0 2 3 792 2 nonsingular so Newton s iteration converges quadratically to this zero of p as eXpected But Prof W Kahan Math 128B Page 34 Filename SingNewt February 5 2007 750 am MD 0 has asingular p39 L74 g now v 1 0 w jz and v fquotzww 5 0 so Newton s iterates Xk tend to this double zero z of p linearly like z BZkw Try it Newton s iterating function NpX X 7 p39X 1 pX may tend to in nity as X tends to the parabola on which detp39X 0 The equation 7detp39 ET l772 SE 7 51 0 of the n parabola is solved by its parametrization 5 07 7 817 7 5172 and T 111 5172 7 317 7 4 NpX becomes indeterminate at two points on this parabola One is the double zero z 74 and the second is the point X 44 As X approaches a third point 31388 on the parabola the direction of NpX7X 7p39X 1 pX approaches tangency with the parabola elsewhere than near these three points the direction of NpX7X throws Newton s iterate NpX violently across the parabola as X approaches it Parabola on thch det p39xy D I I 1 I I Example q For column 2vector arguments X let column 2vector qX a BX x xv ch112 in which Czx aj B t2 C1B C2 t 2 Now q o and q39 45 O has rank 0 so the analysis displayed above cannot eXplain the linear convergence of Newton s iteration to this zero z Worse qX 0 all along the line whose equation is 1 1 39X 72 and thereon eXcept at this zero z q39X is a nonzero scalar multiple of 1 11 1 so v and w are nonzero scalar multiples of 1 1 whence v f quotz ww 0 Consequently the analysis displayed above does not eXplain why Newton s iterates Xk converge to no zero of q on that line other than the zero z t6 and converges to it like Z BZkw I Try it Prof W Kahan Math 128B Page 44

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

#### "Selling my MCAT study guides and notes has been a great source of side revenue while I'm in school. Some months I'm making over $500! Plus, it makes me happy knowing that I'm helping future med students with their MCAT."

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "It's a great way for students to improve their educational experience and it seemed like a product that everybody wants, so all the people participating are winning."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.