Thompson - Computing for Scientists and Engineers (523188), страница 28
Текст из файла (страница 28)
Show that the trapezoid formula error is consistent with an h 3dependence, as in the estimate (4.44), and with a value of the second derivativeof about -0.5, just half the maximum-magnitude value of -1 (at x = 0).(c) Similarly, for the Simpson formula error show that there is an h 5 dependence, consistent with estimate (4.47), and show that an average fourth derivative of about 0.5 (again half the maximum value) would make a good matchbetween prediction and calculation. nTABLE 4.11 Trapezoid and Simpson integrals for the cosine function from -2 to 2 with variousstepsizes h.
The analytical integral has value 2 sin (2) = 1.8 1859.In Section 4.5 we investigated the numerical first derivative of the cosine function, and we showed (Table 4.8 and Figure 4.7) that an accuracy of a few parts in105 is attained for a method comparable to the Simpson integrator (central derivatives, CCD) only for h of order 10-2, whereas for integration h = 0.25 givesabout the same accuracy (Table 4.11). This contrast emphasizes the much greaterdifficulty of computing accurate numerical derivatives than numerical integrals.144NUMERICAL DERIVATIVES AND INTEGRALSFIGURE 4.9 Errors in the integral of the cosine from -2 to 2 for various stepsizes h.
Thetrapezoid-rule errors are indicated by the dashed line, and the Simpson-rule errors are shown by thesolid line.Higher-order polynomial integrationThe numerical-integration technique of approximating a function by a polynomial ofa given order, then evaluating the integral of this polynomial can be made quite general. If one does not require an estimate of the error involved, then a Lagrange interpolating polynomial, as described in Sections 6.1 and 7.1 of Maron’s book, may beuseful.
Various recipes for polynomial integration algorithms, with error estimatesof the kind that we derived, are derived in Chapter 4 of the text by Nakamura andlisted in the handbook of Abramowitz and Stegun. The spline method of integrationwe describe in Section 5.5 allows a low-order polynomial to be used for each stepsize, but accuracy is achieved by requiring continuity of values and derivativesacross adjacent steps.The technique of evaluating the integrals as a function of h, then taking the limitas h0 analytically, is called Romberg integration, and it can be very effective,especially if the integrands are difficult to calculate, so that the number of integrandvalues that one wishes to evaluate is to be minimized. Romberg integration is discussed in detail in Vandergraft’s text on numerical computation.If unequal steps of the integrand are allowed and the integrand can be evaluatedto arbitrary accuracy at any points, then the technique of Gaussian quadrature is capable of providing great accuracy efficiently.
Chapter 4 of Nakamura’s text on applied numerical methods has a description of such integration algorithms, and Chapter 4 of the numerical recipes book of Press et al. has a function that performsGaussian quadrature. A compendium of numerical-integration methods by Davisand Rabinowitz includes modem techniques (such as the use of splines) and verycomplete coverage of the whole spectrum of methods.4.7 ELECTROSTATIC POTENTIAL FROM A CHARGED WIRE1454.7 PROJECT 4B: ELECTROSTATIC POTENTIALFROM A CHARGED WIREIt is interesting to compute the electrostatic potential from a charged wire because itinvolves analysis, numerics, and application to a scientific problem. It is also simpleenough that you can understand the science without extensive physics background.Suppose that we have electric charge distributed along a finite segment of the y axis,as shown in the insert to Figure 4.10.In the figure insert the charge is indicated as being uniformly distributed alongthe wire, but we assume only that its distribution, (y), goes from y = -1 toy = 1 and is zero for y beyond this and for any nonzero x.
Our task is to computethe potential at any point in space, (x, y), not including the wire itself. Such pointsare called field points.From Coulomb’s law of electrostatic interactions, a small element of charge on asegment of lengthof wire located near y’ on the wire produces at (x, y) a potentialwhere r is the distance between the charge element andthe field point.
(In various systems of electrical units there will be an overall factorbetween charge and potential, which we ignore, but electricians should not.) Thetotal potential at this point is obtained by superimposing the potentials from all pointsalong the wire in the limit that the lengths of the segments shrink to zero. Thus, thepotential is the Riemann integral of (y‘)/r. Explicitly, the potential is given by(4.48)in which we used the theorem of Pythagoras to express r in terms of the coordinates.FIGURE 4.10 Potential from a wire having a uniform distribution of charge along the wire, asindicated in the insert. The potential is shown as a function of x distance from the wire for threedifferent heights, y, above the center of the wire.146NUMERICAL DERIVATIVES AND INTEGRALSNow that the physics problem has been set up, we should concern ourselveswith evaluating this integral.
Suppose that the charge density has a simple powerlaw dependence on position along the wire as(4.49)where p is any integer that is not negative. The potential formula can readily betransformed so that it is expressed in terms of some simpler integrals.Exercise 4.25(a) Make a change of variables in (4.48) and (4.49) in order to simplify the denominator in (4.48), namely, set(4.50)(x = 0 doesn’t occur in this problem), then show that the scaled potential,Vs, isgiven by(4.5 1)where the binomial expansion is used to expand (y + XZ)P. The units ofand therefore of Vs, depend upon the exponent p.
The integrals Ik are definedby(4.52)in which k is an integer. The potential problem has been replaced by the actualproblem of evaluating these integrals.(b) Show (by your own analytical skills, by consulting a table of integrals, orby using a computer-algebra system) that the indefinite integrals in (4.52) fork = 0 and 1 are given by(4.53)(4.54)(c) To derive the analytical form of the integral for higher values of k, write thek-1factor z k in (4.52) as z x z . Then use integration by parts, the result for I1in (4.54), and some algebraic manipulations, to derive the recurrence relation4.7 ELECTROSTATIC POTENTIAL FROM A CHARGED WIRE147(4.55)Check this result by differentiating both sides with respect to z in order to verifythat you recover the integrand in (4.52).
nFrom the results of this exercise we have the component integrals that can be combined by using (4.51), after inserting the limits in (4.52), to produce the total scaledpotential Vs at any field point (x, y).We now have a choice of methods for calculating the electrostatic potential, analytical or numerical. The former method, although it has the elegant results from Exercise 4.25 at its command, is restricted to integer powers p. On the other hand anypositive p may be used with numerical integration of (4.48), provided that |y'|,rather than y', is used in (4.49). However analytical integration is useful for testingthe accuracy of the numerical methods, so we start with it.Potentials by analytical integrationThe simplest analytical potentials are obtained if the electric charge is distributed uniformly along the wire.
Then in the above integrals we have k = p = 0 only. Thescaled potential is then given by(4.56)This potential is plotted in Figure 4.10 as a function of x for three values of y, halfway up the wire (y = 0.5), above the level of the wire (y = 1.5), and far above thewire (y = 2.5). As you would guess from Coulomb’s law of the inverse-distancedependence of the potential, Vs dies away rapidly as the field point is moved fartheraway from the source of the charge on the wire.If the electric charge is distributed linearly along the wire (which would have tobe made of a nonconducting material, so that the charges didn’t flow and neutralizeeach other), with positive charges along the top half and negative charges along thebottom half, then the scaled potential can be calculated from the integral expansion(4.5 1).Exercise 4.26(a) Write a small C function to evaluate the analytical potential, Vs (x,y), for thelinear charge distribution (p = 1) in terms of the integrals Ik obtained in Exercise 4.25 and the expression (4.51) for Vs.(b) Make a plot similar to Figure 4.10 for the potential due to this linear distribution of charge.
nWith two check cases in hand, we are ready to develop the numerical applications.148NUMERICAL DERIVATIVES AND INTEGRALSPotentials by numerical-integration methodsIn the preceding subsection we showed that analytical expressions for the potentialsfrom a line charge with various (positive-integer) power-law distributions for thecharge, given by (4.49), can be derived. In the numerical-integration methods weuse these as test cases, then explore the situation with noninteger power laws.First we need a way to handle the numerics.
It is straightforward to modify Program4.5, Trapezoid and Simpson Integrals, to do this. The resulting program, Electrostatic Potentials by Numerical Integration,is given in thefollowing program.PROGRAM 4.6 Electrostatic potentials by numerical integration, trapezoidal and Simpsonformulas, for a uniform charge distribution.#include#include<stdio.h><math.h>main(){/* Electrostatic Potentialsby Numerical Integration */double X,Y,zmin,zmax,h,VAnl,VByTrap,error,VBySimp;int n;doubleVs(),SimpInt(),TrapInt(),y();printf("Potentialsx = 1;while(X!=O)byNumericalIntegration\n");printf("\n\nInput X, Y, n (X=0 to end): ");scanf("%lf%lf%i",&X,&Y,&n) ;if (X == 0){printf("\nEnd Potentials by Numerical Integration");exit (0) ;n = 2*(n/2); /* truncate n to nearest even value */zmin = (X-1)/Y; zmax = (X+1)/Y;h = (zmax-zmin)/n;VAnl = Vs(X,Y);printf("Analytica1potential=%lg",VAnl);VByTrap = TrapInt(zmin,n,h);/* trapezoid rule */error = VAnl-VByTrap;printf("\nTrapezoidintegral=%lg",VByTrap);printf(" Error=%lg",error);VPySimp = SimpInt(zmin,n,h);/* Simpson rule */4.7 ELECTROSTATIC POTENTIAL FROM A CHARGED WIREerror = VAnl-VBySimp;printf("\nSimpsonintegral=%lg",VBySimp);printf(" Error=%lg",error);} /*end x loop*/double SimpInt(a,n,h)/* Simpson integral */double a,h;int n;double sum;double y();int n2,j;n2 = n/2;if ( ( 2*n2 != n ) || ( n2 == 1) ){printf("\n !! SimpInt has n=%i; not allowed",n);return 0;}elsesum= 0;for ( j = 1; j < n2; j++ )sum=sum+y(a+2*j*h)+2*y(a+(2*j-l)*h);return h*(y(a)+4*y(a+(n-l)*h)+y(a+n*h)+2*sum)/3;double TrapInt(a,n,h)/* Trapezoid integral */double a,h;int n;double sum;double y*();int j;sum = ( y(a) + y(a+n*h) )/2;for ( j = 1; j < n; j++ )sum = sum+y(a+j*h);}return h*sum;149150NUMERICAL DERIVATIVES AND INTEGRALSdouble y(x)/* Potential integrand value*/double x;double prod;prod = l/sqrt(x*x+l);return prod;double Vs(x,y)/* Analytical potential integral */double x,y;double sum;sum = -log(sqrt(pow((x+l)/y,2)+1)-(x+1)/y)-log(sqrt(pow((x-1)/y,2)+l)+(x-l)/y);return sum;In Program 4.6 there has been some renaming of variables in the main program,slight modifications of input, but no changes to integrator, functions SimpInt andTrapInt.
These two functions may therefore be copied from their parent, Program 4.5. The price paid for this is that the variable names do not agree with thoseused in formulating our electrostatic-potential problem, so care is needed if the programs are modified. Note that the method of signalling termination of program execution, input x = 0, is convenient rather than necessary. Provided that |y| > 1, sothe field point is not touching the wire (how shocking!), the potential is well-definedon the x axis.
If you don’t like this way of stopping execution, devise your own.The functions for the potential integrand, y (x) , and for the analytical potential,Vs (x , y ) , have been coded to agree with (4.48) and (4.56), respectively, for thecase p = 0. Now you should be ready to code and test the potential program.Exercise 4.27(a) Code and test Electrostatic Potentials by Numerical Integration.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
