Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 7
Текст из файла (страница 7)
When all have been done, we print themif desired, move all entries of YOUT back to YIN, increase x by h and repeat.32The Numerical Solution of Differential EquationsThe principal building block in this structure would be a subroutine that would advancethe solution exactly one step. The main program would initialize the arrays, call thissubroutine, increase x, move date from the output array YOUT back to the input array YIN,print, etc. The single-step subroutine is shown below.
We will use it later on to help usget a solution started when we use methods that need information at more than one pointbefore they can start the integration.Eulerstep:=proc(xin,yin,h,n) local i,yout;# This program numerically integrates the system# y’=f(x,y) one step forward by Euler’s method using step size# h. Enter with values at xin in yin.
Exit with values at xin+h# in yout. Supply f as a function subprogram.yout:=[seq(evalf(yin[i]+h*f(xin,yin,i)),i=1..n)];return(yout);end;A few remarks about the program are in order. One structured data type in Mapleis a list of things, in this case, a list of floating point numbers. The seq command (for“sequence”) creates such a list, in this case a list of length n since i goes from 1 to n in theseq command.
The brackets [ and ] convert the list into a vector. The evalf commandensures that the results of the computation of the components of yout are floating pointnumbers.Our next remark about the program concerns the function subprogram f, which calculates the right hand sides of the differential equation. This subprogram, of course, must besupplied by the user. Here is a sample of such a program, namely the one that describesthe system (2.3.9).
In that case we have f1 (x, y) = x + y1 + y2 and f2 (x, y) = y1 y2 + 1. Thistranslates into the following:f:=proc(x,y,i);# Calculates the right-hand sides of the system of differential# equations.if i=1 then return(x+y[1]+y[2]) else return(y[1]*y[2]+1) fi;end;Our last comment about the program to solve systems is that it is perfectly possible touse it in such a way that we would not have to move the contents of the vector YOUT backto the vector YIN at each step.
In other words, we could save N move operations, where Nis the number of equations. Such savings might be significant in an extended calculation.To achieve this saving, we write two blocks of programming in the main program. Oneblock takes the contents of YIN as input, advances the solution one step by Euler’s methodand prints, leaving the new vector in YOUT.
Then, without moving anything, another blockof programming takes the contents of YOUT as input, advances the solution one step, leavesthe new vector in YIN, and prints. The two blocks call Euler alternately as the integrationproceeds to the right. The reader might enjoy writing this program, and thinking abouthow to generalize the idea to the situation where the new value of y is computed from two2.3 Systems and equations of higher order33previously computed values, rather than from just one (then three blocks of programmingwould be needed).Now we’ve discussed the numerical solution of a single differential equation of first order,and of a system of simultaneous differential equations of first order, and there remains thetreatment of equations of higher order than the first.
Fortunately, this case is very easilyreduced to the varieties that we have already studied.For example, suppose we want to solve a single equation of the second order, sayy 00 + xy 0 + (x + 1) cos y = 2.(2.3.13)The strategy is to transform the single second-order equation into a pair of simultaneousfirst order equations that can then be handled as before. To do this, choose two unknownfunctions u and v.
The function u is to be the unknown function y in (2.3.13), and thefunction v is to be the derivative of u. Then u and v satisfy two simultaneous first-orderdifferential equations:(u0 = v(2.3.14)v 0 = −xv − (x + 1) cos u + 2and these are exactly of the form (2.3.5) that we have already discussed!The same trick works on a general differential equation of N th ordery (N ) + G(x, y, y 0 , y 00 , . .
. , y (N −1) ) = 0.(2.3.15)We introduce N unknown functions u0 , u1 , . . . , uN −1 , and let them be the solutions of thesystem of N simultaneous first order equationsu00 = u1u01 = u2...u0N −2u0N −1(2.3.16)= uN −1= −G(x, u0 , u1 , . . . , uN −2 ).The system can now be dealt with as before.Exercises 2.31. Write each of the following as a system of simultaneous first-order initial-value problems in the standard form (2.3.2):(a) y 00 + x2 y = 0; y(0) = 1; y 0 (0) = 0(b) u0 + xv = 2; v 0 + euv = 0; u(0) = 0; v(0) = 0(c) u0 + xv 0 = 0; v 0 + x2 u = 1; u(1) = 1; v(1) = 0(d) y iv + 3xy 000 + x2 y 00 + 2y 0 + y = 0; y(0) = y 0 (0) = y 00 (0) = y 000 (0) = 1(e) x000 (t) + t3 x(t) + y(t) = 0; y 00 (t) + x(t)2 = t3 ; x(0) = 2; x0 (0) = 1; x00 (0) = 0;y(0) = 1; y 0 (0) = 034The Numerical Solution of Differential Equations2.
For each of the parts of problem 1, write the function subprogram that will computethe right-hand sides, as required by the Eulerstep subroutine.3. For each of the parts of problem 1, assemble nad run on the computer the Eulerprogram, together with the relevant function subprogram of problem 2, to print outthe solutions for fifty steps of integration, each of size h = 0.03. Begin with x = x0 ,the point at which the initial data was given.4. Reprogram the Eulerstep subroutine, as discussed in the text, to avoid the movementof YOUT back to YIN.5.
Modify your program as necessary (in Maple, take advantage of the plot command)to produce graphical output (graph all of the unknown functions on the same axes).Test your program by running it with Euler as it solves y 00 + y = 0 with y(0) = 0,y 0 (0) = 1, and h = π/60 for 150 steps.6. Write a program that will compute successive values yp , yp+1 , . . .
from a differenceequation of order p. Do this by storing the y’s as they are computed in a circular list,so that it is never necessary to move back the last p computed values before findingthe next one. Write your program so that it will work with vectors, so you can solvesystems of difference equations as well as single ones.2.4How to document a programOne 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 serve asguides to the preparation of readable software.In this section we discuss further the all-important question of program documentation,already touched upon in section 2.2.
Some very nontrivial skills are called for in the creationof good user-oriented program descriptions. One of these is the ability to enter the head ofanother person, the user, and to relate to the program that you have just written throughthe user’s eyes.It’s hard to enter someone else’s head. One of the skills that make one person a betterteacher than another person is of the same kind: the ability to see the subject matter thatis being taught through the eyes of another person, the student. If one can do that, oreven make a good try at it, then obviously one will be able much better to deal with thequestions that are really concerning the audience. Relatively few actually do this to anygreat extent not, I think, because it’s an ability that one either has or doesn’t have, butbecause few efforts are made to train this skill and to develop it.We’ll try to make our little contribution here.2.4 How to document a program35(A) What does it do?The first task should be to describe the precise purpose of the program.
Put yourself inthe place of a potential user who is looking at a particular numerical instance of the problemthat needs solving. That user is now thumbing through a book full of program descriptionsin the library of a computer center 10,000 miles away in order to find a program that willdo the problem in question. Your description must answer that question.Let’s now assume that your program has been written in the form of a subroutine orprocedure, rather than as a main program. Then the list of global, or communicatingvariables is plainly in view, in the opening statement of the subroutine.As we noted in section 2.2, you should state the purpose of your program using theglobal variables of the subroutine in the same sentence.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
