c16-6 (779598)
Текст из файла
734Chapter 16.Integration of Ordinary Differential EquationsNote that for compatibility with bsstep the arrays y and d2y are of length 2n for asystem of n second-order equations. The values of y are stored in the first n elements of y,while the first derivatives are stored in the second n elements. The right-hand side f is storedin the first n elements of the array d2y; the second n elements are unused. With this storagearrangement you can use bsstep simply by replacing the call to mmid with one to stoermusing the same arguments; just be sure that the argument nv of bsstep is set to 2n. Youshould also use the more efficient sequence of stepsizes suggested by Deuflhard:(16.5.6)and set KMAXX = 12 in bsstep.CITED REFERENCES AND FURTHER READING:Deuflhard, P.
1985, SIAM Review, vol. 27, pp. 505–535.16.6 Stiff Sets of EquationsAs soon as one deals with more than one first-order differential equation, thepossibility of a stiff set of equations arises. Stiffness occurs in a problem wherethere are two or more very different scales of the independent variable on whichthe dependent variables are changing. For example, consider the following setof equations [1]:u0 = 998u + 1998vv0 = −999u − 1999v(16.6.1)with boundary conditionsu(0) = 1v(0) = 0(16.6.2)v = −y + z(16.6.3)By means of the transformationu = 2y − zwe find the solutionu = 2e−x − e−1000xv = −e−x + e−1000x(16.6.4)If we integrated the system (16.6.1) with any of the methods given so far in thischapter, the presence of the e−1000x term would require a stepsize h 1/1000 forthe method to be stable (the reason for this is explained below). This is so evenSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).n = 1, 2, 3, 4, 5, . . .73516.6 Stiff Sets of EquationsyFigure 16.6.1. Example of an instability encountered in integrating a stiff equation (schematic). Hereit is supposed that the equation has two solutions, shown as solid and dashed lines. Although the initialconditions are such as to give the solid solution, the stability of the integration (shown as the unstabledotted sequence of segments) is determined by the more rapidly varying dashed solution, even after thatsolution has effectively died away to zero.
Implicit integration methods are the cure.though the e−1000x term is completely negligible in determining the values of u andv as soon as one is away from the origin (see Figure 16.6.1).This is the generic disease of stiff equations: we are required to follow thevariation in the solution on the shortest length scale to maintain stability of theintegration, even though accuracy requirements allow a much larger stepsize.To see how we might cure this problem, consider the single equationy0 = −cy(16.6.5)where c > 0 is a constant.
The explicit (or forward) Euler scheme for integratingthis equation with stepsize h isyn+1 = yn + hyn0 = (1 − ch)yn(16.6.6)The method is called explicit because the new value yn+1 is given explicitly interms of the old value yn . Clearly the method is unstable if h > 2/c, for then|yn | → ∞ as n → ∞.The simplest cure is to resort to implicit differencing, where the right-hand sideis evaluated at the new y location. In this case, we get the backward Euler scheme:0yn+1 = yn + hyn+1oryn+1 =yn1 + ch(16.6.7)(16.6.8)The method is absolutely stable: even as h → ∞, yn+1 → 0, which is in fact thecorrect solution of the differential equation.
If we think of x as representing time,then the implicit method converges to the true equilibrium solution (i.e., the solutionat late times) for large stepsizes. This nice feature of implicit methods holds onlyfor linear systems, but even in the general case implicit methods give better stability.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).x736Chapter 16.Integration of Ordinary Differential EquationsOf course, we give up accuracy in following the evolution towards equilibrium ifwe use large stepsizes, but we maintain stability.These considerations can easily be generalized to sets of linear equations withconstant coefficients:(16.6.9)where C is a positive definite matrix. Explicit differencing givesyn+1 = (1 − Ch) · yn(16.6.10)Now a matrix An tends to zero as n → ∞ only if the largest eigenvalue of Ahas magnitude less than unity.
Thus yn is bounded as n → ∞ only if the largesteigenvalue of 1 − Ch is less than 1, or in other wordsh<2λmax(16.6.11)where λmax is the largest eigenvalue of C.On the other hand, implicit differencing givesyn+1 = yn + hy0n+1(16.6.12)yn+1 = (1 + Ch)−1 · yn(16.6.13)orIf the eigenvalues of C are λ, then the eigenvalues of (1 + Ch)−1 are (1 + λh)−1 ,which has magnitude less than one for all h. (Recall that all the eigenvaluesof a positive definite matrix are nonnegative.) Thus the method is stable for allstepsizes h.
The penalty we pay for this stability is that we are required to inverta matrix at each step.Not all equations are linear with constant coefficients, unfortunately! Forthe systemy0 = f(y)(16.6.14)yn+1 = yn + hf(yn+1 )(16.6.15)implicit differencing givesIn general this is some nasty set of nonlinear equations that has to be solved iterativelyat each step. Suppose we try linearizing the equations, as in Newton’s method:"#∂f · (yn+1 − yn )(16.6.16)yn+1 = yn + h f(yn ) +∂y ynHere ∂f/∂y is the matrix of the partial derivatives of the right-hand side (the Jacobianmatrix). Rearrange equation (16.6.16) into the form−1∂f· f(yn )yn+1 = yn + h 1 − h∂y(16.6.17)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).y0 = −C · y16.6 Stiff Sets of Equations737If h is not too big, only one iteration of Newton’s method may be accurate enoughto solve equation (16.6.15) using equation (16.6.17). In other words, at each stepwe have to invert the matrix1−h∂f∂y(16.6.18)In both the routines to be given in this section, we have explicitly carried out thisreplacement for you, so the routines can handle right-hand sides of the form f(y, x)without any special effort on your part.We now mention an important point: It is absolutely crucial to scale your variables properly when integrating stiff problems with automatic stepsize adjustment.As in our nonstiff routines, you will be asked to supply a vector yscal with whichthe error is to be scaled.
For example, to get constant fractional errors, simply setyscal = |y|. You can get constant absolute errors relative to some maximum valuesby setting yscal equal to those maximum values. In stiff problems, there are oftenstrongly decreasing pieces of the solution which you are not particularly interestedin following once they are small. You can control the relative error above somethreshold C and the absolute error below the threshold by settingyscal = max(C, |y|)(16.6.20)If you are using appropriate nondimensional units, then each component of C shouldbe of order unity.
If you are not sure what values to take for C, simply trysetting each component equal to unity. We strongly advocate the choice (16.6.20)for stiff problems.One final warning: Solving stiff problems can sometimes lead to catastrophicprecision loss. Be alert for situations where double precision is necessary.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).to find yn+1 . Solving implicit methods by linearization is called a “semi-implicit”method, so equation (16.6.17) is the semi-implicit Euler method. It is not guaranteedto be stable, but it usually is, because the behavior is locally similar to the case ofa constant matrix C described above.So far we have dealt only with implicit methods that are first-order accurate.While these are very robust, most problems will benefit from higher-order methods.There are three important classes of higher-order methods for stiff systems:• Generalizations of the Runge-Kutta method, of which the most usefulare the Rosenbrock methods.
The first practical implementation of theseideas was by Kaps and Rentrop, and so these methods are also calledKaps-Rentrop methods.• Generalizations of the Bulirsch-Stoer method, in particular a semi-implicitextrapolation method due to Bader and Deuflhard.• Predictor-corrector methods, most of which are descendants of Gear’sbackward differentiation method.We shall give implementations of the first two methods.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
