Deturck, Wilf - Lecture Notes on Numerical Analysis

Lectures on Numerical AnalysisDennis Deturck and Herbert S. WilfDepartment of MathematicsUniversity of PennsylvaniaPhiladelphia, PA 19104-6395Copyright 2002, Dennis Deturck and Herbert WilfApril 30, 20022Contents1 Differential and Difference Equations51.1Introduction .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.2Linear equations with constant coefficients . . . . . . . . . . . . . . . . . . .81.3Difference equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .111.4Computing with difference equations .

. . . . . . . . . . . . . . . . . . . . .141.5Stability theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .161.6Stability theory of difference equations . . . . . . . . . . . . . . . . . . . . .192 The Numerical Solution of Differential Equations232.1Euler’s method . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .232.2Software notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .262.3Systems and equations of higher order . . . . . . . . . . . . . . . . . . . . .292.4How to document a program .

. . . . . . . . . . . . . . . . . . . . . . . . .342.5The midpoint and trapezoidal rules . . . . . . . . . . . . . . . . . . . . . . .382.6Comparison of the methods . . . . . . . . . . . . . . . . . . . . . . . . . . .432.7Predictor-corrector methods . . . .

. . . . . . . . . . . . . . . . . . . . . . .482.8Truncation error and step size . . . . . . . . . . . . . . . . . . . . . . . . . .502.9Controlling the step size . . . . . . . . . . . . . . . . . . . . . . . . . . . . .542.10 Case study: Rocket to the moon . . . . . . . . . . . . .

. . . . . . . . . . .602.11 Maple programs for the trapezoidal rule . . . . . . . . . . . . . . . . . . . .652.11.1 Example: Computing the cosine function . . . . . . . . . . . . . . .672.11.2 Example: The moon rocket in one dimension . . . . . . . . . . . . .682.12 The big leagues . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .692.13 Lagrange and Adams formulas . . . . . . . . . . . . . . . . . . . . . . . . .744CONTENTS3 Numerical linear algebra813.1Vector spaces and linear mappings . . . . . . . . . . . . . . . . . . . . . . .813.2Linear systems . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .863.3Building blocks for the linear equation solver . . . . . . . . . . . . . . . . .923.4How big is zero? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .973.5Operation count . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . .1023.6To unscramble the eggs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .1053.7Eigenvalues and eigenvectors of matrices . . . . . . . . . . . . . . . . . . . .1083.8The orthogonal matrices of Jacobi . .

. . . . . . . . . . . . . . . . . . . . .1123.9Convergence of the Jacobi method . . . . . . . . . . . . . . . . . . . . . . .1153.10 Corbató’s idea and the implementation of the Jacobi algorithm . . . . . . .1183.11 Getting it together . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . .1223.12 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .124Chapter 1Differential and DifferenceEquations1.1IntroductionIn this chapter we are going to study differential equations, with particular emphasis on howto solve them with computers. We assume that the reader has previously met differentialequations, so we’re going to review the most basic facts about them rather quickly.A differential equation is an equation in an unknown function, say y(x), where theequation contains various derivatives of y and various known functions of x.

The problemis to “find” the unknown function. The order of a differential equation is the order of thehighest derivative that appears in it.Here’s an easy equation of first order:y 0 (x) = 0.(1.1.1)The unknown function is y(x) = constant, so we have solved the given equation (1.1.1).The next one is a little harder:y 0 (x) = 2y(x).(1.1.2)A solution will, now doubt, arrive after a bit of thought, namely y(x) = e2x . But, if y(x)is a solution of (1.1.2), then so is 10y(x), or 49.6y(x), or in fact cy(x) for any constant c.Hence y = ce2x is a solution of (1.1.2). Are there any other solutions? No there aren’t,because if y is any function that satisfies (1.1.2) then(ye−2x )0 = e−2x (y 0 − 2y) = 0,(1.1.3)and so ye−2x must be a constant, C.In general, we can expect that if a differential equation is of the first order, then themost general solution will involve one arbitrary constant C.

This is not always the case,6Differential and Difference Equationssince we can write down differential equations that have no solutions at all. We would have,for instance, a fairly hard time (why?) finding a real function y(x) for which(y 0 )2 = −y 2 − 2.(1.1.4)There are certain special kinds of differential equations that can always be solved, andit’s often important to be able to recognize them. Among there are the “first-order linear”equationsy 0 (x) + a(x)y(x) = 0,(1.1.5)where a(x) is a given function of x.Before we describe the solution of these equations, let’s discuss the word linear.

To saythat an equation is linear is to say that if we have any two solutions y1 (x) and y2 (x) of theequation, then c1 y1 (x) + c2 y2 (x) is also a solution of the equation, where c1 and c2 are anytwo constants (in other words, the set of solutions forms a vector space).Equation (1.1.1) is linear, in fact, y1 (x) = 7 and y2 (x) = 23 are both solutions, and sois 7c1 + 23c2 . Less trivially, the equationy 00 (x) + y(x) = 0(1.1.6)is linear. The linearity of (1.1.6) can be checked right from the equation itself, withoutknowing what the solutions are (do it!). For an example, though, we might note thaty = sin x is a solution of (1.1.6), that y = cos x is another solution of (1.1.6), and finally,by linearity, that the function y = c1 sin x + c2 cos x is a solution, whatever the constants c1and c2 .

Now let’s consider an instance of the first order linear equation (1.1.5):y 0 (x) + xy(x) = 0.(1.1.7)So we’re looking for a function whose derivative is −x times the function. Evidently y =22e−x /2 will do, and the general solution is y(x) = ce−x /2 .If instead of (1.1.7) we hady 0 (x) + x2 y(x) = 0,then we would have found the general solution ce−xAs a last example, take3 /3.y 0 (x) − (cos x) y(x) = 0.(1.1.8)The right medicine is now y(x) = esin x . In the next paragraph we’ll give the general ruleof which the above are three examples. The reader might like to put down the book at thispoint and try to formulate the rule for solving (1.1.5) before going on to read about it.Ready? What we need is to choose some antiderivative A(x) of a(x), and then thesolution is y(x) = ce−A(x) .Since that was so easy, next let’s put a more interesting right hand side into (1.1.5), byconsidering the equationy 0 (x) + a(x)y(x) = b(x)(1.1.9)1.1 Introduction7where now b(x) is also a given function of x (Is (1.1.9) a linear equation? Are you sure?).To solve (1.1.9), once again choose some antiderivative A(x) of a(x), and then note thatwe can rewrite (1.1.9) in the equivalent forme−A(x)d A(x)ey(x) = b(x).dxNow if we multiply through by eA(x) we see thatd A(x)ey(x) = b(x)eA(x)dx(1.1.10)so , if we integrate both sides,ZeA(x) y(x) =xb(t)eA(t) dt + const.

,(1.1.11)where on the right side, we mean any antiderivative of the function under the integral sign.ConsequentlyZ xy(x) = e−A(x)b(t)eA(t) dt + const. .(1.1.12)As an example, consider the equationy0 +y= x + 1.x(1.1.13)We find that A(x) = log x, then from (1.1.12) we getZx1y(x) =(t + 1)t dt + Cxx2 x C=+ + .32x(1.1.14)We may be doing a disservice to the reader by beginning with this discussion of certaintypes of differential equations that can be solved analytically, because it would be erroneousto assume that most, or even many, such equations can be dealt with by these techniques.Indeed, the reason for the importance of the numerical methods that are the main subjectof this chapter is precisely that most equations that arise in “real” problems are quiteintractable by analytical means, so the computer is the only hope.Despite the above disclaimer, in the next section we will study yet another importantfamily of differential equations that can be handled analytically, namely linear equationswith constant coefficients.Exercises 1.11.

Find the general solution of each of the following equations:(a) y 0 = 2 cos x8Differential and Difference Equations2y=0x(c) y 0 + xy = 31(d) y 0 + y = x + 5x(e) 2yy 0 = x + 1(b) y 0 +2. Show that the equation (1.1.4) has no real solutions.3. Go to your computer or terminal and familiarize yourself with the equipment, theoperating system, and the specific software you will be using.

Then write a programthat will calculate and print the sum of the squares of the integers 1, 2, . . . , 100. Runthis program.4. For each part of problem 1, find the solution for which y(1) = 1.1.2Linear equations with constant coefficientsOne particularly pleasant, and important, type of linear differential equation is the varietywith constant coefficients, such asy 00 + 3y 0 + 2y = 0 .(1.2.1)It turns out that what we have to do to solve these equations is to try a solution of a certainform, and we will then find that all of the solutions indeed are of that form.Let’s see if the function y(x) = eαx is a solution of (1.2.1).

If we substitute in (1.2.1),and then cancel the common factor eαx , we are left with the quadratic equationα2 + 3α + 2 = 0whose solutions are α = −2 and α = −1. Hence for those two values of α our trial functiony(x) = eαx is indeed a solution of (1.2.1). In other words, e−2x is a solution, e−x is asolution, and since the equation is linear,y(x) = c1 e−2x + c2 e−x(1.2.2)is also a solution, where c1 and c2 are arbitrary constants. Finally, (1.2.2) must be the mostgeneral solution since it has the “right” number of arbitrary constants, namely two.Trying a solution in the form of an exponential is always the correct first step in solvinglinear equations with constant coefficients. Various complications can develop, however, asillustrated by the equationy 00 + 4y 0 + 4y = 0 .(1.2.3)Again, let’s see if there is a solution of the form y = eαx .

This time, substitution into(1.2.3) and cancellation of the factor eαx leads to the quadratic equationα2 + 4α + 4 = 0,(1.2.4)1.2 Linear equations with constant coefficients9whose two roots are identical, both being −2. Hence e−2x is a solution, and of course so isc1 e−2x , but we don’t yet have the general solution because there is, so far, only one arbitraryconstant. The difficulty, of course, is caused by the fact that the roots of (1.2.4) are notdistinct.In this case, it turns out that xe−2x is another solution of the differential equation (1.2.3)(verify this), so the general solution is (c1 + c2 x)e−2x .Suppose that we begin with an equation of third order, and that all three roots turnout to be the same.

For instance, to solve the equationy 000 + 3y 00 + 3y 0 + y = 0(1.2.5)we would try y = eαx , and we would then be facing the cubic equationα3 + 3α2 + 3α + 1 = 0 ,(1.2.6)whose “three” roots are all equal to −1. Now, not only is e−x a solution, but so are xe−xand x2 e−x .To see why this procedure works in general, suppose we have a linear differential equationwith constant coeficcients, sayy (n) + a1 y (n−1) + a2 y (n−2) + · · · + an y = 0(1.2.7)If we try to find a solution of the usual exponential form y = eαx , then after substitution into(1.2.7) and cancellation of the common factor eαx , we would find the polynomial equationαn + a1 αn−1 + a2 αn−2 + · · · + an = 0 .(1.2.8)The polynomial on the left side is called the characteristic polynomial of the givendifferential equation.

Suppose now that a certain number α = α∗ is a root of (1.2.8) ofmultiplicity p. To say that α∗ is a root of multiplicity p of the equation is to say that(α − α∗ )p is a factor of the characteristic polynomial. Now look at the left side of the givendifferential equation (1.2.7).