10.†  The n^th order Bessel function

arises naturally in many applications. (For example the natural frequencies of vibration of a circular drum skin are simply related to the roots of  J_n(x)  divided by the skin radius.) The above series converges for all  x. The aim of this question is to compare various different approximations for the zeroth order Bessel function  J₀(x). A skeleton program for (i) and (ii) can be obtained via the command

cpi  $mc3fp/bessel.inc  bessel.f90

Complete this program according to the specifications below. [If selected for handing in, then submit only part (i), or part (ii), that follow, but not both.]

(i)  Write a STATEMENT FUNCTION called SERIES8(X) that estimates  J₀(x)  by using the series (1) up to degree 8 terms. Write an external FUNCTION SUBPROGRAM called SERIES(X) that estimates  J₀(x)  by using the partial sum of  M+1  terms in (1), i.e. up to degree  2M. The subprogram should use a recurrence relation to generate successive terms, and choose M such that the partial sums up to the  M^th  and  (M+1)^th  terms agree to 4 leading digits. This requires stopping the summation loop when the magnitude of the last term added on, i.e. the  (M+1)^th, is less than 0.5E-4 times the magnitude of the last partial sum. Write and run a program – call the main program BESSEL – that includes SERIES8 and SERIES, and produces a 4-column table, with meaningful headings, showing the values of X, SERIES8(X), SERIES(X), and the  J₀ estimate returned by the subroutine BSSLJ in the NSWC (US Naval Surface Warfare Centre) scientific subroutine library. The table should display results for X = 0.,2.,4.,6.,...,30., using format G11.4 for the  J₀  estimates. Notes: BSSLJ does not use (1), but an alternative approximation. The formula in BSSLJ has a maximum error less than 10^-14 for the range of X being considered, but the accuracy actually obtained depends on the number of digits employed in calculating BSSLJ. Access the NSWC library via

f90o  bessel.f90  $nswclib

and CALL BSSLJ(Z,0,W) in BESSEL. Subroutine BSSLJ(Z,N,W) takes complex Z, integer N, and returns a complex value W for  J_N(Z) . So your main program must declare Z and W to be COMPLEX, and assign Z = CMPLX(X,0.) (or equivalently just Z = X using mixed mode) inside the X-loop before the CALL to BSSLJ. Since  J₀(X) is real for real X, your table should just display REAL(W) as the estimate for  J₀ . Note that X is real, but the X-loop index must be an integer.

(ii)^A Proceed as in Part (i) but with the following modifications.

(a)  Change the argument list of SERIES to X, M. Reference to SERIES(X,M) then returns the value of SERIES and also of M (2M being the highest power used to construct SERIES.

(b)  Copy your FUNCTION SERIES(X,M) subprogram to the end of the entire program file, change the subprogram copy's name to DSERIES(X,M), and modify it internally to reflect this change.

(c)  In your main program declare DSERIES to be DOUBLE PRECISION and introduce another double precision variable DX that will range through the values 0D0 to 30D0 (the double precision extension of the range 0 to 30) as X ranges from 0. to 30.0 (i.e. 0E0 to 30E0). In DSERIES, declare all REAL variables (including X and DSERIES) to be double precision - use IMPLICIT or explicit type declarations - see fortut8. Most constants in DSERIES must also be double precision, otherwise they might reduce the value of DSERIES to single precision accuracy. INTEGER constants are automatically converted to double precision in mixed mode expressions, but REAL constants like 1., if any, must be changed to 1D0 etc.

(d)  Your final program should produce a 7-column table, with meaningful headings, showing the values of X, SERIES8(X), SERIES(X,M1), the number of terms used in SERIES, DSERIES(DX,M2), the number of terms in DSERIES, and finally the J₀ estimate returned by BSSLJ.

[You should observe as follows:- (a) The SERIES8 approximation is inaccurate for x > 2 due to truncation error caused by stopping at degree 8. (b) The SERIES approximation converges, but is inaccurate for x > 12, due to roundoff error. Such errors will be discussed in Part 2 of the course. (g) The DSERIES approximation agrees with NSWC, and is accurate across the whole range considered. However, note that the use of higher precision as in DSERIES merely defers the effect of roundoff error to values of  x  higher than those considered here.]