Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 68
Текст из файла (страница 68)
. ,y ( N - 1 ) (x))(8.79)this reduction can always be accomplished as follows: With y1 = y, we sety´1 = y2y´2 = y3y´3 = y4= . . .(8.80)y´ N - 1 = y Ny´N = f(x, y1, y2, . . . yN)The system (8.80) is equivalent to (8.79). Not every system of equationswill be expressible in the simple form of (8.80). More generally, a system ofN first-order equations will have the formy´1 = f1(x, y1, y2, .
. . ,yN)y´ = f (x, y , y , . . . ,y ). .2 . . 2. . . 1. . 2. . . . . .N .y´N = fN(x, y1, y2, . . . ,yN)(8.81)All the numerical methods considered in this chapter can be adapted to thesystem (8.81). The system (8.81) can be expressed more compactly invector form,y´ = f(x, y)where y´, f, and y are vectors with N components.We illustrate the procedure for the Runge-Kutta method for two*8.12SYSTEMS OF DIFFERENTIAL EQUATIONS399equations, which we write in the formy´ = f(x, y, z)z´ = g(x, y, z)(8.82)The Runge-Kutta formulas corresponding to (8.37) will now be(8.83)whereExtension to a system of equations is obvious.
Note that all the incrementswith lower subscript must be computed before proceeding to those of nexthigher subscript.The Adams-Moulton formulas adapted to the pair of equations (8.82)proceed as follows:(8.84)400THE SOLUTION OF DIFFERENTIAL EQUATIONSIn Sec. 8.6 we described a subroutine named DVERK from the IMSLprograms and used it to solve a single differential equation. Here we willuse this subroutine to solve a system of first-order differential equations.
InDVERK, X will denote the independent variable while Y(K), K =1 . . . , N is used to denote the vector of dependent variables of length Nassuming that we have a system of N first-order equations in the form(8.81). YPRIME(K), K = 1, . . . , N is used to denote the vector of functions f1, . . .
, fN in the right-hand side of (8.81). The subroutine FCN isused to define YPRIME(K). Usage of DVERK for a system of equationsis otherwise identical to its usage for a single equation. The example belowillustrates this usage.Example 8.9 Express the following system of equations as a system of first-orderequations and solve it from x = 0 to x = 1 using the subroutine DVERK:(8.85)z (0) = z´(0) = 0y(0) = 1y´(0) = -2In this example x is the independent variable while z(x) and y(x) are the dependentvariables.
To express this as a first-order system we set z(x) = y1(x), y(x) = y2(x) andthen the first-order system together with the initial conditions becomesy´1 (x) = y 3 (x)y´2 (x) = y 4 (x)y 1 (0) = 0.0y 2 (0) = 1.0y´3 (X ) = y 1 2 (x) - y 2 (x) + e x2y´ 4 (x) = y 1 ( X ) - y 2 (x) - exy 3 (0) = 0.0(8.86)y 4 (0) = -2.0The FORTRAN program and partial results are given below. The values arecorrect to at least eight significant digits. It appears that 16 function evaluations peroutput step were required to achieve this accuracy. This implies that an internal step ofroughly h = 0.05 was used.CPROGRAM TO SOLVE EXAMPLE 6.9 USING D V E R K(IMSL).INTEGER IER,IND,K,N,NWREALC(24),TOL,W(5,9),X,XEND,Y(4)DATA N , X ,, TOL,IND,NW*/ 4 , 0., 0.,l.,0.,-2.,l.E-9, 1 , 5 /EXTERNAL FCN2DO 12 K=1,10XEND = FLOAT(K)/l0.CALL DVERK ( N, FCN2, X, Y, XEND, TOL, IND, C, NW, W, IER )PRINT 600, XEND,Y(l),Y(2),C(24)600FORMAT(3X,F3.1,3X,2(E16.8,3X),F4.0)12 CONTINUESTOPENDSUBROUTINE FCN2 ( N, X, Y, YPRIME )INTEGER NREAL X, Y(N), YPRIME(N)YPRIME(1) = Y(3)YPRIME(2) = Y(4)YPRIME(3) = Y(l)**2 - Y(2) + EXP(X)YPRIME(4) = Y(1) - Y(2)**2 - EXP(X)RETURNEND*8.13STIFF DIFFERENTIAL EQUATIONS401COMPUTER RESULTS FOR EXAMPLE 8.9X.l.2.3.4.5.6.7.8.91.0Y(l)5.12342280E4.19528369E1.44796017E3.50756908E699842327E1.23532042E2.0046026E3.05983760E446147292E6.28019076EY(2)-04030202020101010101FCN EVALS790476884E - 015.63595308E - 013.21283135E - 01644861308E - 02-2.07035152E - 01-494906488E - 01-8.02372169E - 01-1.13460479E + 00-1.49915828E + 00-190666076E + 00163248648096112128144160As this example illustrates, DVERK is a very simple subroutine to use,and it is extremely efficient when high accuracy is required.EXERCISES8.12-l Write the second-order equationy´´(x) = 2(e2x - y2)1/2y (0) = 0y´(0) = 1as a system of first-order equations and solve it from x = 0 to x = 1 using the classicalfourth-order Runge-Kutta method with fixed step sizes of h = 1/64 and h = 1/128.
Estimate theaccuracy of your results.8.12-2 Solve the following second-order equation from x = 1 to x = 2 using the AdamsMoulton formulas (8.84) with a fixed step size of h = 0.1:y´´(x) = 2y 3y (1) = 1y´(1) = -1You will need four starting values for y(x) and f(x,y) = 2y3. Generate these from the exactsolution y(x) = 1/x and then compare your results with the exact solution.8.12-3 Check with your computing center to see if they subscribe to the IMSL programs. Ifthey do, solve the equation in Exercise 8.12-1 using DVERK with the XEND values K/10.with K = l, .
. . ,10.*8.13 STIFF DIFFERENTIAL EQUATIONSApplications in a number of important areas, including chemical reactions,control systems, and electronic networks, lead to systems of differentialequations which are especially difficult to solve because different processesin the system behave with significantly different time scales. If, for example, the solution of a differential equation is given by y(x) = C1 e-x +C 2 e -l00x , the second component of the solution will decay much more402THE SOLUTION OF DIFFERENTIAL EQUATIONSrapidly than the first component as x increases. Most of the methods wehave described for solving differential equations exhibit extreme instabilitywhen applied to problems which have solutions of this type. Problems withsolution components containing widely different time scales are said to bestiff problems.Consider for instance the second order equation(8.87)The general solution of (8.87) isy(x) = Ae - x + Be - 1 0 0 0 xIf we impose the initial conditions y(0) = 1, y´(0) = -1, the exact solutionisy(x) = e-xWe now try to solve (8.87) with these initial conditions using the RK 4method.
The system rewritten as a first-order system (see Sec. 8.12) isy 1 (0) = 1(8.88)y 2 (0) = -1For steps h < 0.002, the Runge-Kutta method yields solutions whichapproximate e-x very nicely. However, h = 0.002 means that we must take500 integration steps per unit interval. Since the desired solution is y(x) =e-x, it would appear safe to take a much larger step h. However, if we takeh = 0.003, still quite small, the numerical solution essentially explodes tocc.
The explanation for this behavior is related to the stability requirements of the method being used. For the RK 4 method, the region ofstability is such that we must have (see Gear [30])1000h < 2.8or h < 0.0028. That is, the step h is for stability reasons restricted by themost rapidly changing component of the solution, namely e-1000y , for theproblem above. Adams-Moulton and other standard multistep methodswould similarly restrict the step h.Extensive research is still going on to find suitable methods for solvingstiff differential equations.
The most successful methods apparently areimplicit. The trapezoidal method (8.52), for example, has been used withsome success. For this method the region of stability is the entire negativehalf-plane, so that h is unrestricted by stability requirements (see Gear[30]). As applied to a system of two equations in two unknowns of theformy´1 = f 1 (x,y 1 ,y 2 )y´2 = f 2 (x,y 1 ,y 2 )*8.13STIFF DIFFERENTIAL EQUATIONS403the trapezoidal method becomes(8.89)Specializing these to the system (8.88), which is linear, leads toNormally, these equations are solved by iteration but because of thelinearity we can obtain an explicit system for the unknowns y1,n+1 andy2,n+l:(8.90)We now choose h = 0.1 so that (8.90) becomesy 1,n+l - 0.05y2,n+1 = y1n + .05y 2 n50y1,n+1 + 51.05y2,n+1 = -49.05y2n - 50y 1 nFor n = 0 we have y10 = 1, y20 = -1, and from (8.91) we obtain(8.91)y11 = 0.904762y 21 = -0.904762which is a reasonable approximation to the exact solution y(0.1) = e-0.1 =0.904837, considering the large step size being used.
After 10 steps withh = 0.1 we obtain y1(1.0)0.367573 which compares very favorably withthe exact result y(1.0) = e-1.0 = 0.367879.In using the trapezoidal method for stiff nonlinear problems, however,there is one essential modification which must be used. For the singleequation y´ = f(x,y) the trapezoidal method is implicit and defined by(8.92)With n fixed, this is an implicit equation which must be solved for yn+1 bysome iterative method. Normally, one uses fixed-point iteration defined bym = 0, l, . . .whereis an approximation to yn+1 obtained by some other methodsuch as Euler’s method.
This fixed-point iteration will converge as shown404THE SOLUTION OF DIFFERENTIAL EQUATIONSin Sec. 8.8 ifand sincefor stiff problems is very large,this requires very small step sizes for convergence. We can, however, solve(8.92) for yn+1 by Newton’s iteration method as follows. We set = y n + 1and rewrite (8.92) in the form(8.93)Ifis an initial approximation tothen successive approximations are generated according to Newton’s method by the iterationm = 0, l, .
. .where from (8.93)In this case there is no difficulty with convergence whenis large andnegative, which is the typical situation with stiff problems. Newton’smethod does, however, require the computation offor a singleequation and of the elements of the Jacobian matrix for a system ofequations. Subroutines for solving stiff differential equations can thereforebe expected to be somewhat complicated.For a system of linear equations of the formwhere A is a constantthe eigenvalues of thewidely separated, thensolving it by ordinaryy´ = Aymatrix, the stiffness of the problem is determined bymatrix A.
If the eigenvalues of A are negative andthe system is stiff and we can expect difficulty inmethods. For the example (8.88) the matrix A isand its eigenvalues are -1000, -1.For more general nonlinear systems of the formy´ = f(x, y)stiffness is determined by the eigenvalues of the Jacobian matrixThe reader is referred to Gear [30] for a more complete discussion of stiffproblems and for other methods for handling them.*8.13STIFF DIFFERENTIAL EQUATIONS405EXERCISES8.13-1 Try to solve the system (8.88) from x = 0 to x = 2 using the Runge-Kutta method oforder 4 for the step-sizes h = 0.001, 0.002, 0.003, 0.01.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
