Thompson - Computing for Scientists and Engineers (523188), страница 75
Текст из файла (страница 75)
Further, for large v, F appears to be about 1/(2 v). Appropriateseries can be found as follows.Exercise 10.40(a) Consider the first-order nonlinear differential equation for F, (10.75). Assume a series solution of odd powers of v, since we know that F (v) is an oddfunction of v. Thus show that(10.77)(b) Express the first derivative of F with respect to v as a first derivative withrespect to 1/v, then use (10.74) to show that this quantity is finite asonly if(10.78)which provides the asymptotic value of Dawson’s integral.
nAlthough the series in v, (10.77), is appropriate for very small v, its convergence istoo slow for much practical use, since v < 2 is a fairly common argument.Numerical integration from the formula (10.73) is a practical method for v values that are not large enough for (10.78) to be accurate. Some care needs to betaken because both large and small exponentials appear in the integration.
Computerunderflow and overflow will therefore be a problem for large enough v. One way toreduce this problem is to rewrite (10.73) as(10.79)The range of exponentials that appear is then exp(-v 2/2) to exp(+v2 /2).10.4 COMPUTING AND APPLYING THE VOIGT PROFILE413Exercise 10.41Show that if the maximum power of 10 that the computer can handle is PMAX,and that a safety factor is 10 is allowed for, then the maximum value of v thatcan be used with the form (10.79) is VLIM =nThis condition is checked in the function FDawson given in Program 10.2.Another consideration in the numerical integration is that if v < 0, the integration can either not be started (v = 0) or should range over decreasing values of u.To avoid the awkward coding involved for the latter case, one may use the fact thatF is an odd function of v, so that F (v) = (v /|v|) F (|v|).
When v is nearly zero,the integral also needs to be set to zero in order to avoid underflow in the exponentialfunction.The program for Dawson’s integral can therefore be written compactly as shownin function FDawson in the complete Voigt Profi1e program in the next subsection. For simplicitly, the trapezoidal rule, (4.43), is used for numerical integration.Exercise 10.42(a) Code and run as a stand-alone program FDawson, with PMAX adjusted foryour computer. Check your programming against numerical values for bothsmall and large v from the series expansions (10.77) and (10.78), respectively.Table 7.5 in Abramowitz and Stegun may also be used for this purpose.(b) Plot Dawson’s integral, F (v), against v for, say, 0 < v < 2.
Check theappearance of your graph against the curve for F (v) in Figure 10.10. nNow that the numerics of Dawson’s integral has been covered, the Hn (v) functions are relatively straightforward to compute.Program for series expansion of profileThe series expansion of the Voigt profile H (a,v) in powers of the ratio of Lorentzian to Gaussian widths, a, is given by (10.66). The functions Hn (v) for n = 0 3 are specified by (10.69), (10.71), (10.70), and (10.76), and for n = 1 and 3 require Dawson’s F function that we have just considered.Exercise 10.43(a) Program the formulas for the Hn (v), using the functions given in Voigtprofi1e below. Check values for even n against hand calculations and thosefor odd n against spot checks for values of F(v) from Exercise 10.42.(b) Plot the Hn (v) functions against v for say 0 < v < 2.
Compare yourcurves with those in Figure 10.10. nNow that we are convinced that the component functions for the Voigt functionare correct, it is time to assemble them to form the series expansion (10.66). Given414FOURIER INTEGRAL TRANSFORMSthat you are using the divide-and-conquer approach to program writing and verification, this is now straightforward (especially if you copy out the Voigt Profilecode that we provide). Note that for FDawson, since each value is used twice, wesignificantly improve the program efficiency by invoking this function only once,then we save the value.
The programming for H1 and H3 must be slightly differentthan in Exercise 10.42, in order to accommodate this sharing of F. We also includethe code section for computing the Voigt function by direct numerical integration ofits Fourier transform as described in the next subsection.The program Voigt Profile includes writing to an output file, "VOIGT" , inorder to save the Voigt function values for graphing.
Also included are comparisonvalues of the appropriate Gaussian and Lorentzian functions that are also output tothe file. Our choice of units, which also reduces the complexity of coding and input,is to set the Gaussian standard deviationPROGRAM 10.2 Voigt profiles by series and by direct integration.#include <stdio.h>#include <math.h>#define MAX 102main()/* Voigt Profile; Series,Integral,Gaussian,Lorentzian*/Gaussian standard deviation set to 1/sqrt(2) *//*FILE *fout;FILE *fopen();double gLseries[MAX],gLtrap[MAX],g[MAX],L[MAX];double pi,sp,gma,vmax,dv,dx,vmin,a,v,even,odd,FD;int Nx,Niv,iv;char wa;double FDawson(),Hzero(),Hone(),Htwo(),Hthree(),trap();pi = 4.0*atan(l.0);sp = sqrt(pi);printf("Voigt Profiles\n");gma = 1.0; /* gma is Lorentzian FWHM */printf("Write over output (w) or Add on (a): ");scanf("%s",&wa) ; fout = fopen("VOIGT",&wa);while ( gma > 0 )printf ("input gma,vmax,dv,dx,Nx (gma<=0 to end):\n");scanf("%le%lf%le%le%i",&gma,&vmax,&dv,&dx,&Nx);if ( gma > 0 )vmin = -vmax;10.4COMPUTING AND APPLYING THE VOIGT PROFILENiv = (vmax-vmin)/dv+l.l;if ( Niv > MAX+1 )printf("!! Niv=%i > array sizes %i\n",Niv,MAX-1);exit (1);}a = gma/2;v = vmin;for ( iv = 1; iv <= Niv; iv++ ){ /* loop over v values */even = Hzero(v)+a*a*Htwo(v);/* even series */FD = FDawson(v);/* odd series has Dawson */odd = a*(Hone(sp,FD,v)+a*a*Hthree(sp,FD,v));gLseries[iv] = (even+odd)/sp; /* series value *//* Trapezoid integral */gLtrap[iv]=trap(a,v,dx,Nx)/pi;g[iv] = exp(-v*v)/sp;/* Gaussian; sigma=l/sqrt(2) */L[iv] = (a/pi)/(v*v+a*a);/* normed Lorentzian;a=gma/2 */printf("%6.2lf %10.6lf %10.6lf %10.6lf %10.6lf\n",v,gLseries[iv],gLtrap[iv],g[iv],L[iv]);fprintf(fout, "%6.2lf %10.6lf %10.6lf %10.6lf %10.6lf\n",v,gLseries[iv],gLtrap[iv],g[iv],L[iv]);v=v+dv;} /*end v loop*/} /*end gma loop */printf("\nEnd Voigt Profile");double FDawson(v)/* Dawson's integral approximated by trapezoid rule */double v;{double PMAX,VLIM,absv,vsq2,expv,trap,u,du,signv;int Nu,iu;absvdu =Nu =PMAXVLIMif (= fabs(v);0.005; /* gives 5 sig.fig.
accuracy for v <= 2 */du = absv/Nu; /* more-exact du */absv/du;= 10; /* maximum power of 10 for computer */= sqrt(2*log(lO)*(PMAX-1)); /* max v */absv < pow(l0,-PMAX) )trap = 0;return trap;415416FOURIER INTEGRAL TRANSFORMS}elsesignv = v/fabs(v);if ( absv > VLIM )printf("!! |v|=%8.2lf set to VLIM=%8.2lf",absv,VLIM);absv = VLIM;vsq2 = absv*absv/2; expv = exp(-vsq2);trap = expv/2;u = du;for ( iu = 1; iu < Nu; iu++ )trap = trap+exp(u*u-vsq2); u = u+du;trap = signv*du* (expv*trap+0.5);return trap;double Hzero(v) /* Zeroth-order H function */double v;double h;h = exp(-v*v);return h;double Hone(sp,FD,v) /* First-order H function */double sp,FD,v;double h;h = 2*(2*v*FD-l)/sp;return h;}double Htwo(v) /* Second-order H function */double v;double h;h = (l-2*v*v)*exp(-v*v);return h;10.4 COMPUTING AND APPLYING THE VOIGT PROFILEdouble Hthree(sp,FD,v) /* Third-order H function */double sp,FD,v;{double h;h = 4*( (3-2 *v*v)*v*FD+v*v-1)/ (3*sp);return h;double trap(a,v,dx,Nx)/* Trapezoid rule direct integral for Voigt profile */double a,v,dx;int Nx;double x, sum;int ix;x = 0;sum = 0.5;for ( ix = 1; ix <= Nx; ix++ ){x = x+dx;sum = sum+cos(v*x)/exp(x*(a+x/4));}sum = dx*sum;return sum;}FIGURE 10.11 The Voigt profile.
(10.65). with ratio-of-widths parameter a = 0.25.417418FOURIER INTEGRAL TRANSFORMSFigure 10.11 illustrates the Voigt profile with a normalized Lorentzian for theparameter a = 0.25, corresponding to a ratio of Lorentzian to Gaussian widths forwhich convergence of the series expansion to about 1% accuracy is expected.
Asanticipated from our discussion of the convolution of a Gaussian with a Lorentzianand from Figure 10.11, the Voigt profile resembles a Gaussian with wings broadened by the Lorentzian tails. Since the convolution preserves area, according to(10.43), the peak value must be slightly reduced by convolution, as you see.Now that you are convinced that the analysis and program are correct, it is agood time for you to develop a complete program for the Voigt profile.Exercise 10.44(a) Adapt Voigt Profile for your computer system.
Note that the parameterPMAX in function FDawson is computer-dependent, as discussed in Exercise 10.41. For check values of the functions Hn consult Exercise 10.43.(b) Plot the output function H (a,v) against v for interesting values of a that aresomewhat less than unity, in order to assure reasonable convergence of the series expansion (10.66). Compare with the curves given in Figure 10.11.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
