Contains reports on three research projects.Joint Services Electronics Programs (U.S. Army, U.S. Navy, and U.S. Air Force) under Contract DA 28-043-AMC-02536(E

Brown, Elaine S.

Isaacs, Elaine C.

Kleiman, Terry S.

River, Eleanor C.

DSpace@MIT

XV. COMPUTATION RESEARCHResearch StaffEleanor C. River Joan Harwitt Eloise KlebanElaine S. Brown Elaine C. Isaacs Terry S. KleimanA. READING BINARY TAPES IN UNDEFINED RECORD FORMATON THE 360/65In the Research Laboratory of Electronics there are several specialized data-takingdevices utilizing magnetic tape for permanent storage of digitized data. For instance,Potter Model 906 II-1 tape drives are used by both the flying-spot scanner and thePDP-4 to ensure physical compatibility with the IBM tape units used by the 360/65. Thedata format of such tapes is typically incompatible with Fortran IV (release 11 andhigher) standard binary tape format, which requires two full-words of control informa-tion at the beginning of each record.A general assembly language program which reads and writes the binary tapesdescribed above is now available at the M. I. T. Computation Center. The Job Controllanguage permits the user to set all valid OS/360 DD parameters, such as recordingdensity, record length, and block size, without program modification. Users oftendesire to unpack the data for easier internal processing; at present, our group hasadditional assembly language programs for packing and unpacking 12 bit and 18 bit datawords.Eleanor C. River, Elaine C. IsaacsB. GAUSSIAN QUADRATURE - A NUMERICAL TECHNIQUEFOR INTEGRATIONWe have developed a set of numerical integration programs in Fortran IV, whichutilize Gaussian Quadrature formulas to evaluate integrals which cannot be evaluatedusing more conventional methods, such as the Trapezoidal Rule or Simpson's Rule.Gaussian Quadratures are numerical approximations of the formnSf(x) dx = w f(a.) + E.j= 1The weights w. and abscissas aj are determined so that En = 0 for a polynomial of degreeless than or equal to 2n - 1. Since by the Weierstrass theorem on polynomial*This work was supported in part by the Joint Services Electronics Programs(U. S. Army, U. S. Navy, and U. S. Air Force) under Contract DA 28-043-AMC-02536(E).QPR No. 87 231(XV. COMPUTATION RESEARCH)approximation, any function f continuous on (a, b) can be approximated by a polynomialof sufficiently high degree, n can be chosen such that given [ > 0, IEni < p for any suchfunction f.There are numerous sets of weights and abscissas, proven mathematically to beappropriate to various types of intervals. For a function f continuous on a finite inter-val [a,b], we have written routines using the following formulas:nf(y) dy = -2 w f(Yi) + En (1)j=ly= x+jxj = the jth zero of the Legendre Polynomial P (x).2W. =b f(y) nSdy = a wf(y_) + E (2)j=1yj = (b-a)x. + a2 thx. = a. where a. is the j positive zero of the Legendre Polynomial Pzn(x).2w. = 2(1-aIPn(a )]2Sab f(y ) ndy =Tb wf(yj) + E (3)NT-b- y j=1yj = (a-b)x. + bx. and w. are as in (2).J jThe creation of these routines involved writing a program to evaluate the coefficientsof the Legendre Polynomials, and setting up other programs to calculate the weightsnwj and to form the sums K j w f(y.). A subprogram to compute f must be suppliedby the user.As an example of the use of (1) and (3), we have evaluatedQPR No. 87 232(XV. COMPUTATION RESEARCH)+sin- (2-1) +sin-1 (Z/l+ siny) 2 - siny sinx I sin xA = dy dx,/2 sin2 x - (2-siny sinx)z1 T 31r< <1, 0<x<r, -<y<2-,where "+ sin-" represents the greater value of the inverse sine in the range of x and ygiven. Let b = +sin (2~/(1+ siny)). The integrand is infinite at x = b. We can establishthat the integrand G(x,y) can be written as-sin-1(2-1) b F(x,y)/ dy dx.Rather than explicitly doing the algebra, we merely write2z_-sinxsiny sin x hJ]-xsin2 x - (Z-sinysinx) F(x,y)G(x,y) =Since the formula (3) does not evaluate F(x,y) at b, and sincelim =. ,x-b sinZ x - (2 -siny sinx)we can safely use (3) with f = F(x, y ) for a given yj. The procedure for the double inte-b F(x, y)gration is to use (1) with f = /ZNrb-x dx and at each yj defined in (1) apply (3) asexplained above.Elaine S. BrownC. BESSEL FUNCTIONS OF THE SECOND KINDA frequently occurring equation in physics is Bessel's equation,22 d dy 2 2x - + x x+ (x -n )y = 0.dxzHistorically, the development of solutions to this equation was slow and erratic.The earliest appearance of an equation of this type occurs in a paper published byJohn Bernouilli, in 1694. Seventy years later, an investigation by Euler on the vibra-tions of a stretched membrane reveals the earliest analysis of a Bessel function ofintegral order. Throughout the next half-century, the analysis of this equation wasQPR No. 87 233(XV. COMPUTATION RESEARCH)continued by such men as Lagrange, Laplace, and Fourier. Then, in a memoirwritten in 1824, the German mathematician, Friedrich Wilhelm Bessel, put fortha detailed discussion of the functions that now bear his name. His solution, theBessel coefficients of argument z and order n, has come to be denoted univer-sally by the symbol J (z).The desire for a second independent solution, and hence a fundamental sys-tem of solutions, led many mathematicians to work further on the problem. Theresult was several versions of the Bessel function of the second kind. No stand-ard form has yet been accepted universally. American authors commonly chooseas a standard particular solution Neumann's function of order n, but there remainssome discrepancy in the exact definition of this function. Some authors cite aNeumann function denoted by Y(n)(z)S(n) log- n-4 2 (n-m-1) n! J (z) ) (-1)m-1 (n+2m)(-) + Jnm(Z), (1)(nZ) = J I((z){logz - - m m(n+m) (n+2m)m=O (nm)(m!)(z m=lwhere1 1 1s s 1 +-+ ...- m= 1,2...o m 2 3 mThen, since any linear combination of J n(z) and Y(n)(z) is also a solution, these authorsdefine a second solution, Y(n) (z), the Weber function, asYn(Z) =2 Y(n)(z)+ (y-log 2)Jn(z)where y is Euler's constant = 0. 577216. Other authors prefer to substitute (1) into (2),combine terms, and to call the resulting function Neumann's function, denoted by Yn (z),wherenm m-1 s2m n1 2mn (-1) -s + ns z z-n (n-m-l)! z22 r z m n  zm-)Sn 2 m=0 2m+n )(m!)(m+n)! I (2m-n)(mz > 0; n = 0,1,..... (2)In order to be as complete and as helpful as possible, we felt that a subroutinewritten to compute both Yn(z), the Neumann function, and Y (z), the Weber function,would be desirable. Working with definitions (1) and (2) as stated above, the subroutineYNEU was written in Fortran IV for the IBM 360. For real arguments z of order n,real and positive, the routine returns values of Y, the Neumann function, and BY,the Weber function.QPR No. 87 234(XV. COMPUTATION RESEARCH)The subroutine calls two other subroutines, BESJ and BESY, supplied in version IIof the scientific subroutine package for the IBM 360.The routine has been written for five place accuracy and has been tested and checkedout for values of z between 0. 02 and 12. 0 and n = 0, 1,Z... 10.Terry S. KleimanQPR No. 87 235

Computation Research

http://dspace.mit.edu/bitstream/1721.1/55671/1/RLE_QPR_087_XV.pdf

Computation Research

Abstract

Similar works

Full text

Available Versions

DSpace@MIT