c16-4 (779596), страница 3

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 3)

However, our feeling is that this view is guidedmore by the kinds of problems used for tests than by one method being actually“better.” Accordingly, we provide the optional routine rzextr for rational functionextrapolation, an exact substitution for pzextr above.#include "nrutil.h"extern float **d,*x;Defined in bsstep.void rzextr(int iest, float xest, float yest[], float yz[], float dy[], int nv)Exact substitute for pzextr, but uses diagonal rational function extrapolation instead of polynomial extrapolation.{int k,j;float yy,v,ddy,c,b1,b,*fx;fx=vector(1,iest);x[iest]=xest;if (iest == 1)for (j=1;j<=nv;j++) {yz[j]=yest[j];d[j][1]=yest[j];Save current independent variable.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).void pzextr(int iest, float xest, float yest[], float yz[], float dy[], int nv)Use polynomial extrapolation to evaluate nv functions at x = 0 by fitting a polynomial to asequence of estimates with progressively smaller values x = xest, and corresponding functionvectors yest[1..nv]. This call is number iest in the sequence of calls.

Extrapolated functionvalues are output as yz[1..nv], and their estimated error is output as dy[1..nv].{int k1,j;float q,f2,f1,delta,*c;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);.

Характеристики

Список файлов книги