c16-5 (779597)
Текст из файла
732Chapter 16.Integration of Ordinary Differential Equations}CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§7.2.14. [1]Gear, C.W. 1971, Numerical Initial Value Problems in Ordinary Differential Equations (EnglewoodCliffs, NJ: Prentice-Hall), §6.2.Deuflhard, P. 1983, Numerische Mathematik, vol.
41, pp. 399–422. [2]Deuflhard, P. 1985, SIAM Review, vol. 27, pp. 505–535. [3]16.5 Second-Order Conservative EquationsUsually when you have a system of high-order differential equations to solve it is bestto reformulate them as a system of first-order equations, as discussed in §16.0. There isa particular class of equations that occurs quite frequently in practice where you can gainabout a factor of two in efficiency by differencing the equations directly. The equations aresecond-order systems where the derivative does not appear on the right-hand side:y00 = f (x, y),y(x0 ) = y0 ,y0 (x0 ) = z0(16.5.1)As usual, y can denote a vector of values.Stoermer’s rule, dating back to 1907, has been a popular method for discretizing suchsystems.
With h = H/m we havey1 = y0 + h[z0 + 12 hf (x0 , y0 )]yk+1 − 2yk + yk−1 = h2 f (x0 + kh, yk ),zm = (ym − ym−1 )/h +1hf (x02+ H, ym )k = 1, . . . , m − 1(16.5.2)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).dy[j]=yest[j];}else {for (k=1;k<iest;k++)fx[k+1]=x[iest-k]/xest;for (j=1;j<=nv;j++) {Evaluate next diagonal in tableau.v=d[j][1];d[j][1]=yy=c=yest[j];for (k=2;k<=iest;k++) {b1=fx[k]*v;b=b1-c;if (b) {b=(c-v)/b;ddy=c*b;c=b1*b;} elseCare needed to avoid division by 0.ddy=v;if (k != iest) v=d[j][k];d[j][k]=ddy;yy += ddy;}dy[j]=ddy;yz[j]=yy;}}free_vector(fx,1,iest);16.5 Second-Order Conservative Equations733Here zm is y0 (x0 + H). Henrici showed how to rewrite equations (16.5.2) to reduce roundofferror by using the quantities ∆k ≡ yk+1 − yk .
Start with∆0 = h[z0 + 12 hf (x0 , y0 )]y 1 = y 0 + ∆0(16.5.3)∆k = ∆k−1 + h2 f (x0 + kh, yk )yk+1 = yk + ∆k(16.5.4)Finally compute the derivative fromzm = ∆m−1 /h + 12 hf (x0 + H, ym )(16.5.5)Gragg again showed that the error series for equations (16.5.3)–(16.5.5) contains onlyeven powers of h, and so the method is a logical candidate for extrapolation à la Bulirsch-Stoer.We replace mmid by the following routine stoerm:#include "nrutil.h"void stoerm(float y[], float d2y[], int nv, float xs, float htot, int nstep,float yout[], void (*derivs)(float, float [], float []))Stoermer’s rule for integrating y00 = f (x, y) for a system of n = nv/2 equations. On inputy[1..nv] contains y in its first n elements and y0 in its second n elements, all evaluated atxs. d2y[1..nv] contains the right-hand side function f (also evaluated at xs) in its first nelements.
Its second n elements are not referenced. Also input is htot, the total step to betaken, and nstep, the number of substeps to be used. The output is returned as yout[1..nv],with the same storage arrangement as y. derivs is the user-supplied routine that calculates f .{int i,n,neqns,nn;float h,h2,halfh,x,*ytemp;ytemp=vector(1,nv);h=htot/nstep;Stepsize this trip.halfh=0.5*h;neqns=nv/2;Number of equations.for (i=1;i<=neqns;i++) {First step.n=neqns+i;ytemp[i]=y[i]+(ytemp[n]=h*(y[n]+halfh*d2y[i]));}x=xs+h;(*derivs)(x,ytemp,yout);Use yout for temporary storage of derivatives.h2=h*h;for (nn=2;nn<=nstep;nn++) {General step.for (i=1;i<=neqns;i++)ytemp[i] += (ytemp[(n=neqns+i)] += h2*yout[i]);x += h;(*derivs)(x,ytemp,yout);}for (i=1;i<=neqns;i++) {Last step.n=neqns+i;yout[n]=ytemp[n]/h+halfh*yout[i];yout[i]=ytemp[i];}free_vector(ytemp,1,nv);}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).Then for k = 1, . . . , m − 1, set734Chapter 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, . .
..
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
