Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 6
Текст из файла (страница 6)
If jobswithin a large program are not broken into subroutines it can be much harder to isolate theblock of coding that deals with a particular function and remove it without affecting thewhole works.For another reason, if one be needed, it may well happen that even though the job thatis done by the subroutine occurs only once in the current program, it may recur in otherprograms as yet undreamed of. If one is in the habit of writing small independent modulesand stringing them together to make large programs, then it doesn’t take long before onehas a library of useful subroutines, each one tested, working and documented, that willgreatly simplify the task of writing future programs.Finally, the practice of subdividing the large jobs into the smaller jobs of which they arecomposed is an extremely valuable analytical skill, one that is useful not only in programming, but in all sorts of organizational activities where smaller efforts are to be pooled inorder to produce a larger effect.
It is therefore a quality of mind that provides much of itsown justification.In this book, the major programs that are the objects of study have been broken up intosubroutines in the expectation that the reader will be able to start writing and checking outthese modules even before the main ideas of the current subject have been fully explained.This was done in part because some of these programs are quite complex, and it would beunreasonable to expect the whole program to be written in a short time.
It was also doneto give examples of the process of subdivision that we have been talking about.For instance, the general linear algebra program for solving systems of linear simultaneous equations in Chapter 3, has been divided into six modules, and they are described insection 3.3. The reader might wish to look ahead at those routines and to verify that eventhough their relationship to the whole job of solving equations is by no means clear now,nonetheless, because of the fact that they are independent and self-contained, they can beprogrammed and checked out at any time without waiting for the full explanation.One more ingredient that is needed for the production of useful software is:2.3 Systems and equations of higher order293. StyleWe don’t mean style in the sense of “class,” although this is as welcome in programmingas it is elsewhere.
There have evolved a number of elements of good programming style,and these will mainly be discussed as they arise. But two of them (one trivial and one quitedeep) are:(a) Indentation: The instructions that lie within the range of a loop are indented in theprogram listing further to the right than the instructions that announce that the loopis about to begin, or that it has just terminated.(b) Top-down structuring: When we visualize the overall logical structure of a complicated program we see a grand loop, within which there are several other loops andbranchings, within which .
. . etc. According to the principles of top-down design thelooping and branching structure of the program should be visible at once in the listing. That is, we should see an announcement of the opening of the grand loop, thenindented under that perhaps a two-way branch (if-then-else), where, under the “then”one sees all that will happen if the condition is met, and under the “else” one seeswhat happens if it is not met.When we say that we see all that will happen, we mean that there are not any “go-to”instructions that would take our eye out of the flow of the if-then-else loop to someother page.
It all happens right there on the same page, under “then” and under“else”.These few words can scarcely convey the ideas of structuring, which we leave to thenumerous examples in the sequel.2.3Systems and equations of higher orderWe have already remarked that the methods of numerical integration for a single firstorder differential equation carry over with very little change to systems of simultaneousdifferential equations of first order.
In this section we’ll discuss exactly how this is done,and furthermore, how the same idea can be applied to equations of higher order than thefirst. Euler’s method will be used as the example, but the same transformations will applyto all of the methods that we will study.In Euler’s method for one equation, the approximate value of the unknown function atthe next point xn+1 = xn + h is calculated fromyn+1 = yn + hf (xn , yn ).(2.3.1)Now suppose that we are trying to solve not just a single equation, but a system of Nsimultaneous equations, sayyi0 (x) = fi (x, y1 , y2 , . .
. yN )i = 1, . . . , N.(2.3.2)30The Numerical Solution of Differential EquationsEquation (2.3.2) represents N equations, in each of which just one derivative appears,and whose right-hand side may depend on x, and on all of the unknown functions, but noton their derivatives. The “fi ” indicates that, of course, each equation can have a differentright-hand side.Now introduce the vector y(x) of unknown functionsy(x) = [y1 (x), y2 (x), y3 (x), . . . , yN (x)](2.3.3)and the vector f = f (x, y(x)) of right-hand sidesf = [f1 (x, y), f2 (x, y), . . .
, fN (x, y)].(2.3.4)In terms of these vectors, equation (2.3.2) can be rewritten asy 0 (x) = f (x, y(x)).(2.3.5)We observe that equation (2.3.5) looks just like our standard form (2.1.1) for a singleequation in a single unknown function, except for the bold face type, i.e., except for thefact that y and f now represent vector quantities.To apply a numerical method such as that of Euler, then, all we need to do is to takethe statement of the method for a single differential equation in a single unknown function,and replace y(x) and f (x, y(x)) by vector quantities as above.
We will then have obtainedthe generalization of the numerical method to systems.To be specific, Euler’s method for a single equation isyn+1 = yn + hf (x, yn )(2.3.6)so Euler’s method for a system of differential equations will bey n+1 = y n + hf (xn , y n ).(2.3.7)This means that if we know the entire vector y of unknown functions at the point x = xn ,then we can find the entire vector of unknown functions at the next point xn+1 = xn + hby means of (2.3.7).In detail, if y i (xn ) denotes the computed approximate value of the unknown function yiat the point xn , then what we must calculate arey i (xn+1 ) = y i (xn ) + hfi (xn , y 1 (xn ), y 2 (xn ), .
. . , y N (xn ))(2.3.8)for each i = 1, 2, . . . , N .As an example, take the pair of differential equations(y10 = x + y1 + y2y20 = y1 y2 + 1together with the initial values y1 (0) = 0, y2 (0) = 1.(2.3.9)2.3 Systems and equations of higher order31Now the vector of unknown functions is y = [y1 , y2 ], and the vector of right-hand sidesis f = [x + y1 + y2 , y1 y2 + 1]. Initially, the vector of unknowns is y = [0, 1]. Let’s choose astep size h = 0.05.
Then we calculate"and"y 1 (0.10)y 2 (0.10)y1 (0.05)y2 (0.05)#"=0.051.05#"=01##+ 0.05"+ 0.05"0+0+10∗1+10.05 + 0.05 + 1.050.05 ∗ 1.05 + 1#"=#"=0.051.05#0.10751.102625(2.3.10)#(2.3.11)and so forth. At each step we compute the vector of approximate values of the two unknownfunctions from the corresponding vector at the immediately preceding step.Let’s consider the preparation of a computer program that will carry out the solution,by Euler’s method, of a system of N simultaneous equations of the form (2.3.2), which wewill rewrite just slightly, in the formyi0 = fi (x, y)i = 1, .
. . , N.(2.3.12)Note that on the left is just one of the unknown functions, and on the right there mayappear all N of them in each equation.Evidently we will need an array Y of length N to hold the values of the N unknownfunctions at the current point x. Suppose we have computed the array Y at a point x, andwe want to get the new array Y at the point x + h. Exactly what do we do?Some care is necessary in answering this question because there is a bit of a snare inthe underbrush.
The new values of the N unknown functions are calculated from (2.3.8) or(2.3.12) in a certain order. For instance, we might calculate y 1 (x + h), then y 2 (x + h), theny 3 (x + h),. . . , then yN (x + h).The question is this: when we compute y 1 (x+h), where shall we put it? If we put it intoY[1], the first position of the Y array in storage, then the previous contents of Y[1] arelost, i.e., the value of y1 (x) is lost. But we aren’t finished with y 1 (x) yet; it’s still needed tocompute y 2 (x+h), y 3 (x+h), etc.
This is because the new value y i (x+h) depends (or mightdepend), according to (2.3.8), on the old values of all of the unknown functions, includingthose whose new values have already been computed before we begin the computation ofy i (x + h).If the point still is murky, go back to (2.3.11) and notice how, in the calculation ofy 2 (0.10) we needed to know y 1 (0.05) even though y1 (0.10) had already been computed.Hence if we had put y 1 (0.10) into an array to replace the old value y 1 (0.05) we would nothave been able to obtain y 2 (0.10).The conclusion is that we need at least two arrays, say YIN and YOUT, each of length N .The array YIN holds the unknown functions evaluated at x, and YOUT will hold their valuesat x + h. Initially YIN holds the given data at x0 . Then we compute all of the unknowns atx0 + h, and store them in YOUT as we find them.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
