Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 5
Текст из файла (страница 5)
At each of thesepoints xn we will compute yn , our approximation to y(xn ).Hence, suppose that the spacing h between consecutive points has been chosen. Wepropose to start at the point x0 where the initial data are given, and move to the right,obtaining y1 from y0 , then y2 from y1 and so forth until sufficiently many values have beenfound.Next we need to derive a method by which each value of y can be obtained from itsimmediate predecessor. Consider the Taylor series expansion of the unknown function y(x)24The Numerical Solution of Differential Equationsabout the point xny 00 (X),(2.1.2)2where we have halted the expansion after the first power of h and in the remainder term,the point X lies between xn and xn + h.y(xn + h) = y(xn ) + hy 0 (xn ) + h2Now equation (2.1.2) is exact, but of course it cannot be used for computation becausethe point X is unknown.
On the other hand, if we simply “forget” the error term, we’llhave only an approximate relation instead of an exact one, with the consolation that wewill be able to compute from it. The approximate relation isy(xn + h) ≈ y(xn ) + hy 0 (xn ).(2.1.3)Next define yn+1 to be the approximate value of y(xn+1 ) that we obtain by using theright side of (2.1.3) instead of (2.1.2). Then we getyn+1 = yn + hyn0 .(2.1.4)Now we have a computable formula for the approximate values of the unknown function,because the quantity yn0 can be found from the differential equation (2.1.1) by writingyn0 = f (xx , yn ),(2.1.5)and if we do so then (2.1.4) takes the formyn+1 = yn + hf (xn , yn ).(2.1.6)This is Euler’s method, in a very explicit form, so that the computational procedure isclear. Equation (2.1.6) is in fact a recurrence relation, or difference equation, whereby eachvalue of yn is computed from its immediate predecessor.Let’s use Euler’s method to obtain a numerical solution of the differential equationy 0 = 0.5y(2.1.7)together with the starting value y(0) = 1.
The exact solution of this initial-value problemis obviously y(x) = ex/2 .Concerning the approximate solution by Euler’s method, we have, by comparing (2.1.7)with (2.1.1), f (x, y) = 0.5y, soyn2h= 1+yn .2yn+1 = yn + h(2.1.8)Therefore, in this example, each yn will be obtained from its predecessor by multiplicationby 1 + h2 . To be quite specific, let’s take h to be 0.05. Then we show below, for each valueof x = 0, 0.05, 0.10, 0.15, 0.20, . . . the approximate value of y computed from (2.1.8) andthe exact value y(xn ) = exn /2 :2.1 Euler’s method25x0.000.050.100.150.200.25...Euler(x)1.000001.025001.050631.076891.103811.13141...Exact(x)1.000001.025321.051271.077881.105171.13315...1.002.003.00...1.638622.685064.39979...1.648722.718284.48169...5.0010.0011.81372139.5638912.18249148.41316table 1Considering the extreme simplicity of the approximation, it seems that we have donepretty well by this equation.
Let’s continue with this simple example by asking for a formulafor the numbers that are called Euler(x) in the above table. In other words, exactly whatfunction of x is Euler(x)?To answer this, we note first that each computed value yn+1 is obtained according to(2.1.8) by multiplying its predecessor yn by 1 + h2 . Since y0 = 1, it is clear tha we willcompute yn = (1 + h2 )n .
Now we want to express this in terms of x rather than n. Sincexn = nh, we have n = x/h, and since h = 0.05 we have n = 20x. Hence the computedapproximation to y at a particular point x will be 1.02520x , or equivalentlyEuler(x) = (1.638616 . . .)x .(2.1.9)The approximate values can now easily be compared with the true solution, sincex1Exact(x) = e 2 = e 2x= (1.648721 . . .)x .(2.1.10)Therefore both the exact solution of this differential equation and its computed solutionhave the form (const.)x . The correct value of “const.” is e1/2 , and the value that is, ineffect, used by Euler’s method is (1 + h2 )1/h . For a fixed value of x, we see that if we useEuler’s method with smaller and smaller values of h (neglecting the increase in roundofferror that is sure to result), the values Euler(x) will converge to Exact(x), becausehlim 1 +h→02Exercises 2.11/h1= e2 .(2.1.11)26The Numerical Solution of Differential Equations1.
Verify the limit (2.1.11).2. Use a calculator or a computer to integrate each of the following differential equationsforward ten steps, using a spacing h = 0.05 with Euler’s method. Also tabulate theexact solution at each value of x that occurs.(a) y 0 (x) = xy(x); y(0) = 1(b) y 0 (x) = xy(x) + 2; y(0) = 1y(x)(c) y 0 (x) =; y(0) = 11+x(d) y 0 (x) = −2xy(x)2 ; y(0) = 102.2Software notesOne of the main themes of our study will be the preparation of programs that not onlywork, but also are easily readable and useable by other people. The act of communicationthat must take place before a program can be used by persons other than its author is adifficult one to carry out, and we will return several times to the principles that have evolvedas guides to the preparation of readable software.
Here are some of these guidelines.1. DocumentationThe documentation of a program is the set of written instructions to a user that informthe user about the purpose and operation of the program. At the moment that the jobof writing and testing a program has been completed it is only natural to feel an urge toget the whole thing over with and get on to the next job.
Besides, one might think, it’sperfectly obvious how to use this program. Some programs may be obscure, but not thisone.It is amazing how rapidly our knowledge of our very own program fades. If we comeback to a program after a lapse of a few months’ time, it often happens that we will have noidea what the program did or how to use it, at least not without making a large investmentof time.For that reason it is important that when our program has been written and tested itshould be documented immediately, while our memory of it is still green. Furthermore, thebest place for documentation is in the program itself, in “comment” statements. That wayone can be sure that when the comments are needed they will be available.The first mission of program documentation is to describe the purpose of the program.State clearly the problem that the program solves, or the exact operation that it performson its input in order to get its output.Already in this first mission, a good bit of technical skill can be brought to bear thatwill be very helpful to the use, by intertwining the description of the program purpose withthe names of the communicating variables in the program.Let’s see what that means by considering an example.
Suppose we have written asubroutine that searches through a specified row of a matrix to find the element of largest2.2 Software notes27absolute value, and outputs a column in which it was found. Such a routine, in Maple forinstance, might look like this:search:=proc(A,i)local j, winner, jwin;winner:=-1;for j from 1 to coldim(A) doif (abs(A[i,j])>winner) thenwinner:=abs(A[i,j]) ; jwin:=j fiod;return(jwin);end;Now let’s try our hand at documenting this program:“The purpose of this program is to search a given row of a matrix to find anelement of largest absolute value and return the column in which it was found.”That is pretty good documentation, perhaps better than many programs get. But wecan make it a lot more useful by doing the intertwining that we referred to above.
Therewe said that the description should be related to the communicating variables. Thosevariables are the ones that the user can see. They are the input and output variables ofthe subroutine. In most important computer languages, the communicating variables areannounced in the first line of the coding of a procedure or subroutine. Maple follows thisconvention at least for the input variables, although the output variable is usually specifiedin the “return” statement.In the first line of the little subroutine above we see the list (A, i) of its input variables(or “arguments”). These are the ones that the user has to understand, as opposed to theother “local” variables that live inside the subroutine but don’t communicate with theoutside world (like j, winner, jwin, which are listed on the second line of the program).The best way to help the user to understand these variables is to relate them directlyto the description of the purpose of the program.“The purpose of this program is to search row I of a given matrix A to find anentry of largest absolute value, and returns the column jwin where that entrylives.”We’ll come back to the subject of documentation, but now let’s mention another ingredient of ease of use of programs, and that is:2.
ModularityIt is important to divide a long program into a number of smaller modules, each with aclearly stated set of inputs and outputs, and each with its own documentation. That meansthat we should get into the habit of writing lots of subroutines or procedures, because28The Numerical Solution of Differential Equationsthe subroutine or procedure mode of expression forces one to be quite explicit about therelationship of the block of coding to the rest of the world.When we are writing a large program we would all write a subroutine if we found thata certain sequence of steps was being called for repeatedly. Beyond this, however, there arenumerous inducements for breaking off subroutines even if the block of coding occurs justonce in the main program.For one thing it’s easier to check out the program.
The testing procedure would consistof first testing each of the subroutines separately on test problems designed just for them.Once the subroutines work, it would remain only to test their relationships to the callingprogram.For another reason, we might discover a better, faster, more elegant, or what-have-youmethod of performing the task that one of these subroutines does. Then we would be ableto yank out the former subroutine and plug in the new one, while being careful only to makesure that the new subroutine relates to the same inputs and outputs as the old one.