c6-3 (779502)
Текст из файла
222Chapter 6.Special FunctionsCITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapters 6, 7, and 26.Pearson, K. (ed.) 1951, Tables of the Incomplete Gamma Function (Cambridge: CambridgeUniversity Press).The standard definition of the exponential integral isZ∞En (x) =1e−xtdt,tnx > 0,n = 0, 1, .
. .(6.3.1)The function defined by the principal value of the integralZEi(x) = −∞−xe−tdt =tZx−∞etdt,tx>0(6.3.2)is also called an exponential integral. Note that Ei(−x) is related to −E1 (x) byanalytic continuation.The function En (x) is a special case of the incomplete gamma functionEn (x) = xn−1 Γ(1 − n, x)(6.3.3)We can therefore use a similar strategy for evaluating it. The continued fraction —just equation (6.2.6) rewritten — converges for all x > 0:n1 n+1 21···(6.3.4)En (x) = e−xx+ 1+ x+ 1+ x+We use it in its more rapidly converging even form,1·n2(n + 1)1−x···En (x) = ex+n− x+n+2− x+n+4−(6.3.5)The continued fraction only really converges fast enough to be useful for x >∼ 1.For 0 < x <∼ 1, we can use the series representationEn (x) =(−x)n−1[− ln x + ψ(n)] −(n − 1)!∞Xm=0m6=n−1(−x)m(m − n + 1)m!(6.3.6)The quantity ψ(n) here is the digamma function, given for integer arguments byψ(1) = −γ,ψ(n) = −γ +n−1Xm=11m(6.3.7)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).6.3 Exponential Integrals6.3 Exponential Integrals223where γ = 0.5772156649 .
. . is Euler’s constant. We evaluate the expression (6.3.6)in order of ascending powers of x:xx2(−x)n−21−+−···+En (x) = −(1 − n) (2 − n) · 1 (3 − n)(1 · 2)(−1)(n − 2)!n−1nn+1(−x)(−x)(−x)[− ln x + ψ(n)] −++···+(n − 1)!1 · n!2 · (n + 1)!(6.3.8)e−xx1,En (0) =n−1E0 (x) =(6.3.9)n>1The routine expint allows fast evaluation of En (x) to any accuracy EPSwithin the reach of your machine’s word length for floating-point numbers. Theonly modification required for increased accuracy is to supply Euler’s constant withenough significant digits. Wrench [2] can provide you with the first 328 digits ifnecessary!#include <math.h>#define MAXIT 100#define EULER 0.5772156649#define FPMIN 1.0e-30#define EPS 1.0e-7Maximum allowed number of iterations.Euler’s constant γ.Close to smallest representable floating-point number.Desired relative error, not smaller than the machine precision.float expint(int n, float x)Evaluates the exponential integral En (x).{void nrerror(char error_text[]);int i,ii,nm1;float a,b,c,d,del,fact,h,psi,ans;nm1=n-1;if (n < 0 || x < 0.0 || (x==0.0 && (n==0 || n==1)))nrerror("bad arguments in expint");else {if (n == 0) ans=exp(-x)/x;Special case.else {if (x == 0.0) ans=1.0/nm1;Another special case.else {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).The first square bracket is omitted when n = 1. This method of evaluation has theadvantage that for large n the series converges before reaching the term containingψ(n). Accordingly, one needs an algorithm for evaluating ψ(n) only for small n,n <∼ 20 – 40.
We use equation (6.3.7), although a table look-up would improveefficiency slightly.Amos [1] presents a careful discussion of the truncation error in evaluatingequation (6.3.8), and gives a fairly elaborate termination criterion. We have foundthat simply stopping when the last term added is smaller than the required toleranceworks about as well.Two special cases have to be handled separately:224Chapter 6.Special Functions}}}return ans;}A good algorithm for evaluating Ei is to use the power series for small x andthe asymptotic series for large x.
The power series isEi(x) = γ + ln x +x2x++···1 · 1! 2 · 2!(6.3.10)where γ is Euler’s constant. The asymptotic expansion isEi(x) ∼exx1+1!2!+ 2 +···xx(6.3.11)The lower limit for the use of the asymptotic expansion is approximately | ln EPS|,where EPS is the required relative error.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).if (x > 1.0) {Lentz’s algorithm (§5.2).b=x+n;c=1.0/FPMIN;d=1.0/b;h=d;for (i=1;i<=MAXIT;i++) {a = -i*(nm1+i);b += 2.0;d=1.0/(a*d+b);Denominators cannot be zero.c=b+a/c;del=c*d;h *= del;if (fabs(del-1.0) < EPS) {ans=h*exp(-x);return ans;}}nrerror("continued fraction failed in expint");} else {Evaluate series.ans = (nm1!=0 ? 1.0/nm1 : -log(x)-EULER);Set first term.fact=1.0;for (i=1;i<=MAXIT;i++) {fact *= -x/i;if (i != nm1) del = -fact/(i-nm1);else {psi = -EULER;Compute ψ(n).for (ii=1;ii<=nm1;ii++) psi += 1.0/ii;del=fact*(-log(x)+psi);}ans += del;if (fabs(del) < fabs(ans)*EPS) return ans;}nrerror("series failed in expint");}6.3 Exponential Integrals#include <math.h>#define EULER 0.57721566#define MAXIT 100#define FPMIN 1.0e-30#define EPS 6.0e-8225Euler’s constant γ.Maximum number of iterations allowed.Close to smallest representable floating-point number.Relative error, or absolute error near the zero of Ei atx = 0.3725.if (x <= 0.0) nrerror("Bad argument in ei");if (x < FPMIN) return log(x)+EULER;Special case: avoid failure of convergenceif (x <= -log(EPS)) {test because of underflow.sum=0.0;Use power series.fact=1.0;for (k=1;k<=MAXIT;k++) {fact *= x/k;term=fact/k;sum += term;if (term < EPS*sum) break;}if (k > MAXIT) nrerror("Series failed in ei");return sum+log(x)+EULER;} else {Use asymptotic series.sum=0.0;Start with second term.term=1.0;for (k=1;k<=MAXIT;k++) {prev=term;term *= k/x;if (term < EPS) break;Since final sum is greater than one, term itself approximates the relative error.if (term < prev) sum += term;Still converging: add new term.else {sum -= prev;Diverging: subtract previous term andbreak;exit.}}return exp(x)*(1.0+sum)/x;}}CITED REFERENCES AND FURTHER READING:Stegun, I.A., and Zucker, R.
1974, Journal of Research of the National Bureau of Standards,vol. 78B, pp. 199–216; 1976, op. cit., vol. 80B, pp. 291–311.Amos D.E. 1980, ACM Transactions on Mathematical Software, vol. 6, pp. 365–377 [1]; alsovol. 6, pp. 420–428.Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 5.Wrench J.W. 1952, Mathematical Tables and Other Aids to Computation, vol. 6, p. 255.
[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).float ei(float x)Computes the exponential integral Ei(x) for x > 0.{void nrerror(char error_text[]);int k;float fact,prev,sum,term;226Chapter 6.Special Functions1(0.5,5.0)incomplete beta function Ix (a,b)(8.0,10.0)(1.0,3.0).6(0.5,0.5).4.2(5.0,0.5)00.2.4.6.81xFigure 6.4.1. The incomplete beta function Ix (a, b) for five different pairs of (a, b).
Notice that the pairs(0.5, 5.0) and (5.0, 0.5) are related by reflection symmetry around the diagonal (cf. equation 6.4.3).6.4 Incomplete Beta Function, Student’sDistribution, F-Distribution, CumulativeBinomial DistributionThe incomplete beta function is defined byZ x1Bx (a, b)≡Ix (a, b) ≡ta−1 (1 − t)b−1 dtB(a, b)B(a, b) 0(a, b > 0)(6.4.1)It has the limiting valuesI0 (a, b) = 0I1 (a, b) = 1(6.4.2)and the symmetry relationIx (a, b) = 1 − I1−x(b, a)(6.4.3)If a and b are both rather greater than one, then Ix (a, b) rises from “near-zero” to“near-unity” quite sharply at about x = a/(a + b).
Figure 6.4.1 plots the functionfor several pairs (a, b).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)..8.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
