7  Subroutine H2LC

In this section the Fortran subroutine H2LC is described. The subroutine computes the discrete form of the two-dimensional Helmholtz integral operators. Details of the methods employed in the subroutine are given. The subroutine is applied to the Burton and Miller integral equation for the problem where the boundary is a circle and the results are given.

7.1  Background Methods

The regular integrals that arise are approximated by a standard quadrature rule such as a Gauss-Legendre rule which is specified in the parameter list to the subroutines. Tables of Gauss-Legendre rules are given in Stroud and Secrest (1966) and can also generated from the NAG library [30]. The non-regular integrals that arise in the formulae (39) and (40) are computed via the following methods. See Jaswon and Symm [11] and Kirkup [14] for the background to these methods.

The M0 and M0t operators have regular kernels, hence the aim is to find expressions for:

{ L0 ~
e
 
}DG(p) =
ó
õ
DG 
 G0(p,q)   dSq  ,
(49)

{ N0 ~
e
 
}DG(p; vp) =
vp

ó
õ
DG 
  G0
nq
(p,qdSq   ,
(50)
where DG is a straight line element, p Î DG (though not on an edge or corner of the element). Let it be assumed that the element DG has length a+b with q = q(x) and p = q(51) for x Î [-a,b]. This gives the following formulae for (50) and (51).

{ L0 ~
e
 
}DG(p) = 1
2 p
[ a + b - a loga -b logb ] ,
(52)

{ N0 ~
e
 
}DG(p; vp) = - 1
2 p
[ 1
a
+ 1
b
] .
(53)

7.2  Test Problem and Results

The test problem is in the file H2LC.FOR . The test problem consists of a boundary of a circle of unit radius, centred at the origin. The exterior solution and Neumann boundary condition are determined by a point source at (0,0.5) and a point sink at (0,-0.5). The circle is approximated by 32 straight line elements of equal length, the vertices of the regular polygon lying on the original circle. Eight point Gauss-Legendre quadrature rules are used for all but the computation of the diagonal components of the matrices where 2 ×8 point composite Gaussian-Legendre quadrature rules are used.

The functions FNSQRT and FNLOG simply call the corresponding standard Fortran functions. The subroutine FNHANK was constructed through the use of NAG routines S17AEF, S17AFF, S17ACF and S17ADF and the subroutine is only suitable for real values of k. Note that routines for computing the Hankel functions for complex arguments are now available in Mark 14 of the NAG Fortran library [30]. The subroutine FNHANK is listed in this section.

C SUBROUTINE TO COMPUTE HANKEL FUNCTIONS
C FOR THE MOMENT ONLY Z REAL IS CATERED FOR

      SUBROUTINE FNHANK(Z,H)
      COMPLEX*16 Z,H(0:1)
      COMPLEX*16 J(0:1),Y(0:1)
      COMPLEX*16 CIMAG
      REAL*8 ZERO,ONE
      REAL*8 X
      REAL*8 S17AEF,S17AFF,S17ACF,S17ADF
      INTEGER IFAIL
      IF (ABS(AIMAG(Z)).GT.1.0E-6) STOP
      IFAIL=0
      ZERO=0.0D0
      ONE=1.0D0
      CIMAG=CMPLX(ZERO,ONE)
      X=REAL(Z)
      J(0)=S17AEF(X,IFAIL)
      J(1)=S17AFF(X,IFAIL)
      Y(0)=S17ACF(X,IFAIL)
      Y(1)=S17ADF(X,IFAIL)
      H(0)=J(0)+CIMAG*Y(0)
      H(1)=J(1)+CIMAG*Y(1)
      END

In the example shown here, the circle has radius 1.0, the circle is divided into 16 uniform elements and p* = (0.0,0.5). Since the problem is symmetric about the y axis and antisymmetric about the x axis, the computed and exact solutions are listed at eight collocation points. The computed and exact solutions are also compared at the exterior point (0,2). the wavenumbers considered are k = 0.0 (with a = 1 and b = 0) and k = 1.0 (with a = 1 and b = i). Results are given in tables 1A and 1B.

Table 1A: Solution for the circle with k = 0.0, a = 1, b = 0
Point Exact Solution Computed Solution
1 0.173161 0.173408
2 0.160666 0.160913
3 0.139776 0.140010
4 0.114977 0.115174
5 0.089028 0.089176
6 0.063136 0.063236
7 0.037647 0.037704
8 0.012506 0.012524
(0,2) 0.081300 0.081008

Table 1B: Solution for the circle with k = 1.0, a = 1, b = i
Point Exact Solution Computed Solution
1 0.204820 + i0.106145 0.210022 + i0.106022
2 0.190535 + i0.102052 0.195224 + i0.101934
3 0.166397 + i0.094027 0.170345 + i0.093919
4 0.137368 + i0.082387 0.140541 + i0.082294
5 0.106643 + i0.067589 0.109068 + i0.067513
6 0.075745 + i0.050206 0.077453 + i0.050150
7 0.045200 + i0.030909 0.046214 + i0.030875
8 0.015018 + i0.010435 0.015355 + i0.010424
(0,2) 0.028905 + i0.140053 0.028868 + i0.141126