Heath - Scientific Computing (523150), страница 42
Текст из файла (страница 42)
. . , λn ) U −1 ,where λ1 , . . . , λn are the eigenvalues of A andU is a matrix whose columns are corresponding eigenvectors. Then the matrix exponentialis given by1494.14 Write a routine for solving an arbitrary,possibly rank-deficient, linear least squaresproblem Ax ≈ b using the singular value decomposition. You may call a library routineto compute the SVD, then use its output tocompute the least squares solution (see Section 4.5.2). The input to your routine shouldinclude the matrix A, right-hand-side vectorb, and a tolerance for determining the numerical rank of A. Test your routine on some ofthe linear least squares problems in Chapter 3.4.15 (a) Write a routine that uses a one-sidedplane rotation to symmetrize an arbitrary 2×2matrix. That is, given a 2×2 matrix A, choosec and s so thatc−sCompare your results with those for a libraryroutine for computing the matrix exponential,if you have access to one.
Which of your tworoutines is more accurate and robust? Try toexplain why. See [179] for several additionalmethods for computing the matrix exponential.a11a21a12a22b= 11b12b12b22is symmetric.(b) Write a routine that uses a two-sided planerotation to annihilate the off-diagonal entriesof an arbitrary 2 × 2 symmetric matrix. Thatis, given a symmetric 2 × 2 matrix B, choosec and s so thatc −ssc=exp(A) = U diag(eλ1 , . . . , eλn ) U −1 .Write a second program to evaluate exp(A)using this method.Test both methods using each of the followingtest matrices:2 −1A=,−12113 −114B=.152 −153scb11b12d110b12b220d22c s−s cis diagonal.(c) Combine the two routines developed inparts a and b to obtain a routine for computingthe singular value decomposition A = U ΣV Tof an arbitrary 2 × 2 matrix A. Note thatU will be a product to two plane rotations,whereas V will be a single plane rotation.
Testyour routine on a few randomly chosen 2 × 2matrices and compare the results with thosefor a library SVD routine.By systematically solving successive 2 × 2 subproblems, the module you have just developedcan be used to compute the SVD of an arbitrary m × n matrix in a manner analogous tothe Jacobi method for the symmetric eigenvalue problem.150CHAPTER 4. EIGENVALUES AND SINGULAR VALUES4.16 We will revisit Computer Problem 3.5concerning the elliptical orbit of a planet, represented in a Cartesian (x, y) coordinate system by the equationay 2 + bxy + cx + dy + e = x2 .The orbital parameters a, b, c, d, e can be determined by a linear least squares fit to thefollowing observations of the planet’s position:xyxy1.020.390.560.150.950.320.440.130.870.270.300.120.770.220.160.130.670.180.010.15(a) Use a library routine to compute the singular value decomposition of the 10 × 5 leastsquares matrix.(b) Use the singular value decomposition tocompute the solution to the least squares problem.
With the singular values in order of decreasing magnitude, compute the solutions using the first k singular values, k = 1, . . . , 5.For each of the five solutions obtained, printthe values for the orbital parameters and alsoplot the resulting orbits along with the givendata points in the (x, y) plane.(c) Perturb the input data slightly by addingto each coordinate of each data point a random number uniformly distributed on the interval [−0.005, 0.005] (see Section 13.5).
Compute the singular value decomposition of thenew least squares matrix, and solve the leastsquares problem with the perturbed data as inpart b. Compare the new values for the parameters with those previously computed foreach value of k. What effect does this difference have on the plot of the orbits? Can youexplain this behavior? Which solution wouldyou regard as better: one that fits the datamore closely, or one that is less sensitive tosmall perturbations in the data? Why?4.17 Write a routine for computing the pseudoinverse of an arbitrary m × n matrix. Youmay call a library routine to compute the singular value decomposition, then use its output to compute the pseudoinverse (see Section 4.5.2). Consider the use of a tolerancefor declaring relatively small singular values tobe zero. Test your routine on both singularand nonsingular matrices. In the latter case,of course, your results should agree with thoseof standard matrix inversion.
What happenswhen the matrix is nonsingular, but severelyill-conditioned (e.g., a Hilbert matrix)?4.18 Consider the problem of fitting themodel functionf (t, x) = xt(i.e., a straight line through the origin, withslope x to be determined) to the following datapoints:ty−2−1−133−2(a) Perform a standard linear least squares fitof such a line to y as a function of t, minimizingthe vertical distances between the data pointsand the line (this procedure is appropriate if yis subject to error and t is known exactly).(b) Perform a standard linear least squares fitof such a line to t as a function of y, minimizing the horizontal distances between the datapoints and the line (this procedure is appropriate if t is subject to error and y is knownexactly).(c) Perform a total least squares fit of the lineto the data, minimizing the orthogonal distances between the data points and the line(this procedure is appropriate if both variablesare subject to error).
Such a fit can be done using the singular value decomposition (see Section 4.5.2).(d ) What is the resulting slope x of the line ineach case? Plot the data points and all threelines on a single graph.Chapter 5Nonlinear Equations5.1Nonlinear EquationsWe will now consider methods for solving nonlinear equations. Given a nonlinear functionf , we seek a value x for whichf (x) = 0.Such a solution value for x is called a root of the equation, and a zero of the function f .Though technically they have distinct meanings, these two terms are informally used moreor less interchangeably, with the obvious meaning.
Thus, this problem is often referred toas root finding or zero finding.In discussing numerical methods for solving nonlinear equations, we will distinguish twocases:f : R → R (scalar),andf : Rn → Rn(vector).The latter is referred to as a system of nonlinear equations, in which we seek a vector xsuch that all the component functions of f (x) are zero simultaneously.Example 5.1 Nonlinear Equations. An example of a nonlinear equation in one dimension isf (x) = x2 − 4 sin(x) = 0,for which one approximate solution is x = 1.9. An example of a system of nonlinearequations in two dimensions is 2 x1 − x2 + 0.250f (x) ==,2−x1 + x2 + 0.250for which the solution vector is x = [ 0.50.5 ]T .151152CHAPTER 5. NONLINEAR EQUATIONS5.1.1Solutions of Nonlinear EquationsA system of linear equations always has a unique solution unless the matrix of the systemis singular. The existence and uniqueness of solutions for nonlinear equations are oftenmuch more complicated and difficult to determine, and a much wider variety of behavior ispossible.
Curved lines can intersect, or fail to intersect, in many more ways than straightlines can. For example, unlike straight lines, two curved lines can be tangent without beingcoincident. Whereas for systems of linear equations the number of solutions must be eitherzero, one, or infinitely many, nonlinear equations can have any number of solutions.Example 5.2 Solutions of Nonlinear Equations.
For example:•••••ex + 1 = 0 has no solution.e−x − x = 0 has one solution.x2 − 4 sin(x) = 0 has two solutions.x3 + 6x2 + 11x − 6 = 0 has three solutions.sin(x) = 0 has infinitely many solutions.In addition, a nonlinear equation may have a multiple root, where both the function andits derivative are zero, i.e., f (x) = 0 and f 0 (x) = 0. In one dimension, this property meansthat the curve has a horizontal tangent on the x axis. If f (x) = 0 and f 0 (x) 6= 0, then x issaid to be a simple root.Example 5.3 Multiple Root. Examples of equations having a multiple root includex2 − 2x + 1 = 0and x3 − 3x2 + 3x − 1 = 0,which are illustrated in Fig. 5.1.......................................................................................................................................................................................................................................................................................................................................................Figure 5.1: Nonlinear equations having a multiple root.What do we mean by an approximate solution x̂ to a nonlinear system,kf (x̂)k ≈ 0 or kx̂ − x∗ k ≈ 0,where x∗ is the “true” solution to f (x) = 0? The first corresponds to having a small residual, whereas the second measures closeness to the (usually unknown) true solution.
As with5.1. NONLINEAR EQUATIONS153linear systems, these two criteria for a solution are not necessarily “small” simultaneously.This feature is illustrated for one dimension in Fig. 5.2, where the two functions have aboutthe same uncertainty in their values (e.g., due to rounding error or measurement error) butvery different uncertainties in the locations of their roots (compare with Fig. 2.2). Thus,we see that the same concept of sensitivity or conditioning applies to nonlinear equations:it is the relative change in the solution due to a given relative change in the input data.
Forexample, a multiple root is ill-conditioned, since by definition the curve has a horizontaltangent at such a root, and is therefore locally approximately parallel to the x axis.... ... .... ... .... ... .... ... .... ..... ..... ... ..... ..... ..... ...
......... ..... .... ......... ... .......... .. ......... ....... ....................................... ...................................................................................................................................................................................................................................................................................................................................................................
....... ....... ....... ....... ..... .... .......... ....... ......... .. ............ ..... ..... ... ...... ..... ..... ... ..... ..... .... .. .... .... ...well-conditionedill-conditionedFigure 5.2: Conditioning of roots of nonlinear equations.The conditioning of the root-finding problem for a given function is the opposite of thatfor evaluating the function: if the function value is insensitive to the value of the argument,then the root will be sensitive, whereas if the function value is sensitive to the argument,then the root will be insensitive.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
