Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 3
Текст из файла (страница 3)
Use it tocalculate the largest root of the equationx8 = x7 + x6 + x5 + · · · + 1.1.4(1.3.15)Computing with difference equationsThis is, after all, a book about computing, so let’s begin with computing from differenceequations since they will give us a chance to discuss some important questions that concernthe design of computer programs. For a sample difference equation we’ll useyn+3 = yn+2 + 5yn+1 + 3yn(1.4.1)together with the starting values y0 = y1 = y2 = 1.
The reader might want, just for practice,to find an explicit formula for this sequence by the methods of the previous section.Suppose we want to compute a large number of these y’s in order to verify some propertythat they have, for instance to check thatyn+1lim=3(1.4.2)n→∞ ynwhich must be true since 3 is the root of largest absolute value of the characteristic equation.As a first approach, we might declare y to be a linear array of some size large enough toaccommodate the expected length of the calculation. Then the rest is easy.
For each n, wewould calculate the next yn+1 from (1.4.1), we would divide it by its predecessor yn to geta new ratio. If the new ratio agrees sufficiently well with the previous ratio we announcethat the computation has terminated and print the new ratio as our answer. Otherwise, wemove the new ratio to the location of the old ratio, increase n and try again.If we were to write this out as formal procedure (algorithm) it might look like:1.4 Computing with difference equations15y0 := 1; y1 := 1; y2 := 1; n := 2;newrat := −10; oldrat := 1;while |newrat − oldrat| ≥ 0.000001 dooldrat := newrat; n := n + 1;yn := yn−1 + 5 ∗ yn−2 + 3 ∗ yn−3 ;newrat := yn /yn−1endwhileprint newrat; Halt.We’ll use the symbol ‘:=’ to mean that we are to compute the quantity on the right, ifnecessary, and then store it in the place named on the left.
It can be read as ‘is replacedby’ or ‘is assigned.’ Also, the block that begins with ‘while’ and ends with ‘endwhile’represents a group of instructions that are to be executed repeatedly until the conditionthat follows ‘while’ becomes false, at which point the line following ‘endwhile’ is executed.The procedure just described is fast, but it uses lots of storage. If, for instance, such aprogram needed to calculate 79 y’s before convergence occurred, then it would have used79 locations of array storage. In fact, the problem above doesn’t need that many locationsbecause convergence happens a lot sooner.
Suppose you wanted to find out how muchsooner, given only a programmable hand calculator with ten or twenty memory locations.Then you might appreciate a calculation procedure that needs just four locations to holdall necessary y’s.That’s fairly easy to accomplish, though. At any given moment in the program, whatwe need to find the next y are just the previous three y’s. So why not save only those three?We’ll use the previous three to calculate the next one, and stow it for a moment in a fourthlocation. Then we’ll compute the new ratio and compare it with the old. If they’re notclose enough, we move each one of the three newest y’s back one step into the places wherewe store the latest three y’s and repeat the process. Formally, it might be:y := 1; ym1 := 1; ym2 := 1;newrat := −10; oldrat := 1;while |newrat − oldrat| ≥ 0.000001 doym3 := ym2; ym2 := ym1; ym1 := y;oldrat := newrat;y := ym1 + 5 ∗ ym2 + 3 ∗ ym3;newrat := y/ym1 endwhile;print newrat; Halt.The calculation can now be done in exactly six memory locations (y, ym1, ym2, ym3,oldrat, newrat) no matter how many y’s have to be calculated, so you can undertake iton your hand calculator with complete confidence.
The price that we pay for the memorysaving is that we must move the data around a bit more.One should not think that such programming methods are only for hand calculators. Aswe progress through the numerical solution of differential equations we will see situationsin which each of the quantities that appears in the difference equation will itself be an16Differential and Difference Equationsarray (!), and that very large numbers, perhaps thousands, of these arrays will need to becomputed. Even large computers might quake at the thought of using the first methodabove, rather than the second, for doing the calculation. Fortunately, it will almost neverbe necessary to save in memory all of the computed values simultaneously.
Normally, theywill be computed, and then printed or plotted, and never needed except in the calculationof their immediate successors.Exercises 1.41. The Fibonacci numbers F0 , F1 , F2 , . . . are defined by the recurrence formula Fn+2 =Fn+1 + Fn for n = 0, 1, 2, . . . together with the starting values F0 = 0, F1 = 1.(a) Write out the first ten Fibonacci numbers.(b) Derive an explicit formula for the nth Fibonacci number Fn .(c) Evaluate your formula for n = 0, 1, 2, 3, 4.(d) Prove directly from your formula that the Fibonacci numbers are integers (This isperfectly obvious from their definition, but is not so obvious from the formula!).(e) Evaluatelimn→∞Fn+1Fn(1.4.3)(f) Write a computer program that will compute Fibonacci numbers and print outthe limit in part (e) above, correct to six decimal places.(g) Write a computer program that will compute the first√40 members of the modifiedFibonacci sequence in which F0 = 1 and F1 = (1 − 5)/2.
Do these computednumbers seem to be approaching zero? Explain carefully what you see and whyit happens.(h) Modify the program of part (h) to run in higher (or double) precision arithmetic.Does it change any of your answers?2. Find the most general solution of each of the following difference equations:(a) yn+1 − 2yn + yn−1 = 0(b) yn+1 = 2yn(c) yn+2 + yn = 0(d) yn+2 + 3yn+1 + 3yn + yn−1 = 01.5Stability theoryIn the study of natural phenomena it is most often true that a small change in conditionswill produce just a small change in the state of the system being studied. If, for example,a very slight increase in atmospheric pollution could produce dramatically large changesin populations of flora and fauna, or if tiny variations in the period of the earth’s rotation1.5 Stability theory17produced huge changes in climatic conditions, the world would be a very different place tolive in, or to try to live in.
In brief, we may say that most aspects of nature are stable.When physical scientists attempt to understand some facet of nature, they often willmake a mathematical model. This model will usually not faithfully reproduce all of thestructure of the original phenomenon, but one hopes that the important features of thesystem will be preserved in the model, so that predictions will be possible. One of the mostimportant features to preserve is that of stability.For instance, the example of atmospheric pollution and its effect on living things referredto above is important and very complex. Therefore considerable effort has gone into theconstruction of mathematical models that will allow computer studies of the effects ofatmospheric changes.
One of the first tests to which such a model should be subjected isthat of stability: does it faithfully reproduce the observed fact that small changes producesmall changes? What is true in nature need not be true in a man-made model that is asimplification or idealization of the real world.Now suppose that we have gotten ourselves over this hurdle, and we have constructeda model that is indeed stable. The next step might be to go to the computer and docalculations from the model, and to use these calculations for predicting the effects ofvarious proposed actions.
Unfortunately, yet another layer of approximation is usuallyintroduced at this stage, because the model, even though it is a simplification of the realworld, might still be too complicated to solve exactly on a computer.For instance, may models use differential equations. Models of the weather, of the motion of fluids, of the movement of astronomical objects, of spacecraft, of population growth,of predator-prey relationships, of electric circuit transients, and so forth, all involve differential equations. Digital computers solve differential equations by approximating them bydifference equations, and then solving the difference equations.
Even though the differentialequation that represents our model is indeed stable, it may be that the difference equationthat we use on the computer is no longer stable, and that small changes in initial data onthe computer, or small roundoff errors, will produce not small but very large changes in thecomputed solution.An important job of the numerical analyst is to make sure that this does not happen,and we will find that this theme of stability recurs throughout our study of computerapproximations.As an example of instability in differential equations, suppose that some model of asystem led us to the equationy 00 − y 0 − 2y = 0(1.5.1)together with the initial datay(0) = 1;y 0 (0) = −1.(1.5.2)We are thinking of the independent variable t as the time, so we will be interested inthe solution as t becomes large and positive.The general solution of (1.5.1) is y(t) = c1 e−t + c2 e2t .