SUBROUTINE NEUMAN(IOUT,IDBG,IE,BDYFCN,CURVE,B,DPTS,NPTS,IBEG,EPS0,
+ WORK,IWORK,NWORK,MWORK,MAXFFT,UVEC,ERROR,IER)
C -----------------
C
C This program will solve the interior and exterior Neumann
C problems for Laplace's equation on a simply-connected region D
C with smooth boundary C. It is assumed that C is at least two
C times continuously differentiable. To define a boundary, the user
C must supply a parametric form of the external routine CURVE with
C parameter range [0,B]. The routine CURVE must also produce the
C first and second derivatives on the curve as defined below. The
C boundary condition must be given by the routine BDYFCN.
C
C SUBROUTINE PARAMETERS:
C IOUT Input.
C The output unit number.
C IDBG Input.
C The debugging parameter.
C =0 produces a standard output.
C =1 produces a debugging output.
C IE Input
C =0 for the interior Neumann problem.
C =1 for the exterior Neumann problem.
C BDYFCN External function.
C Inputs : X,Y
C Outputs: BDYFCN
C This is a function of two variables X,Y. It gives
C the value of the normal derivative of the unknown
C harmonic function U on the boundary curve C. This
C is a user supplied function.
C CURVE External subroutine.
C Inputs : S
C Outputs: X,Y,DX,DY,D2X,D2Y
C This is a subroutine which defines the boundary
C curve C. The calling sequence for it is
C CALL CURVE(S,X,Y,DX,DY,D2X,D2Y)
C The variable S is the parameterization variable for
C the curve C. (X,Y) is the point on the curve
C corresponding to the variable S. DX,DY,D2X,D2Y are
C the first and second derivatives of X(S) and Y(S),
C respectively. The routine CURVE is supplied by the user.
C B Input.
C The limits of the variable S in defining the curve C
C are 0 .LE. S .LE. B
C DPTS Input array.
C This is a two dimensional array specifying the points
C at which the harmonic function U is to be evaluated,
C with each column denoting a distinct point. The
C dimension statement for DPTS is
C DIMENSION DPTS(4,NPTS)
C For point #J, its coordinates have the following
C meaning:
C X=DPTS(1,J), Y=DPTS(2,J)
C S=DPTS(3,J), IBD=DPTS(4,J)
C If (X,Y) is not on the boundary C, S need not be
C set and the user must set IBD=0.
C If (X,Y) is on the boundary, we need S to be the
C parameter on [0,B] corresponding to (X,Y); and the
C user must set IBD=1.
C NPTS Input.
C This is the number of points in DPTS. It must be
C greater than zero.
C IBEG Input.
C =0 means this is a first call on NEUMAN for this
C particular curve C, boundary function BDYFCN, and
C error tolerance EPS.
C =1 means that NEUMAN has been called previously for
C this choice of parameters, and the solution U is
C desired at a new set of points in DPTS. For such a
C call on NEUMAN, change only DPTS, NPTS, and IBEG.
C DO NOT CHANGE any other input variable.
C EPS0 Input.
C The desired error tolerance in the solution U.
C WORK Input.
C Temporary work space. This vector should have a
C dimension NWORK of at least 10,000. If some points of
C DPTS are close to the boundary, then the dimension of
C WORK should be increased further. A dimension of
C NWORK=300,000 will allow a wide variety of problems
C to be treated very accurately.
C NWORK Input.
C Dimension of WORK. See the preceding discussion of WORK
C and the general discussion of NWORK given below.
C IWORK Input.
C Temporary work space for integer variables. It is used
C as the pivot array in LU fcatorization by LAPACK
C subroutines in the subroutine INTEQN. It is also used
C in the subroutine EVALU for the FFT routines.
C MWORK Input.
C Dimension of IWORK.
C MWORK .GE. (SQRT(NWORK+49) + 8 )
C MAXFFT Input.
C Maximum dimension of FFT.
C MAXFFT .GE. (2*SQRT(NWORK + 49) - 14 )
C is desirable to solve many problems accurately.
C UVEC Output.
C This is an output vector. On exit, component #I will
C contain the approximate value of the solution function
C U(X(I),Y(I)) at the points (X(I),Y(I)) given in DPTS.
C ERROR Output.
C This is an output vector containing predicted error
C bounds for the solutions given in the vector UVEC.
C The element UVEC(J) is the approximate solution at
C the # J point of DPTS, and ERROR(J) is a predicted
C error bound. If the desired error EPS was not
C attained, then the sign of ERROR(J) is made negative
C for such a solution value. Otherwise the sign of
C ERROR(J) is positive.
C IER Output.
C = 1 means some or all of the solution values in UVEC
C may not satisfy the error tolerance EPS.
C = 0 means the program was completed satisfactorily.
C =-1 means EPS0 < 0.
C =-2 means NWORK < 0 or MWORK < (SQRT(NWORK+49)+8).
C =-3 means B .LE. 0.
C =-4 means IE or IBEG are out of range.
C =-5 means NPTS<=0, an error.
C =-6 means MAXFFT < M_0, where M_0 is defined below.
C
C *** Defining NWORK, the dimension of WORK ***
C Introduce variables MAXSYS and MAXMES. MAXSYS represents
C the maximum order of the linear system to be solved in
C subroutine INTEQN; and MAXMES represents the maximum number
C of mesh points of C at which the single layer density
C function is to be evaluated. MAXFFT represents the maximum
C degree of the Fourier expansion to be produced for the
C approximate density function. As examples,
C MAXSYS=512, MAXMES=8192, MAXFFT = 1024
C will solve many problems very accurately. We define NWORK by
C NWORK=MAX(MAXSYS**2+14*MAXSYS,8*MAXMES+4*MAXFFT)
C The program assumes
C MAXSYS .GE. 64, MAXMES .GE. 128, MAXFFT .GE. 128
C and you should set NWORK accordingly. These defaults can be
C changed by re-setting the respective parameters N_0 and M_0,
C as described below. In particular,
C MAXSYS .GE. 4*N_0, MAXMES .GE. 4*M_0, MAXFFT .GE. M_0
C
C *** THE PROGRAM PARAMETERS K_0, M_0 and N_0 ***
C N_0 Denotes the initial number of quadrature points used in the
C approximation of the integral equation in the subroutine INTEQN.
C K_0 and M_0 denote the initial degree of the FFT and the initial
C number of integration nodes in the subroutine EVALU. Always set
C M_0 and N_0 to be a power of 2, and K_0=N_0 is desirable.
C
C *** SOURCES OF INFORMATION ***
C
C For additional information on this program, see the paper
C
C K. Atkinson & Y. Jeon, "Automatic boundary integral equation
C programs for the planar Laplace equation", ACM Transactions on
C Mathematical Software.
C
C The authors can be contacted at the respective e_mail addresses:
C Kendall-Atkinson@uiowa.edu
C yjeon@madang.ajou.ac.kr
C
C The web site for this program is located at the URL
C http://www.math.uiowa.edu/~atkinson/laplace.html
C
C DATE OF LAST CHANGES TO THIS CODE: 7 April 1998
C
C ================= END OF INTRODUCTORY REMARKS ===================
C
INTEGER IBEG,IDBG,IE,IER,IOUT,MAXFFT,MWORK,NPTS,NWORK
DOUBLE PRECISION DPTS(4,NPTS),ERROR(NPTS),UVEC(NPTS),WORK(NWORK)
INTEGER IWORK(MWORK)
DOUBLE PRECISION B,BDYFCN,EPS0
EXTERNAL BDYFCN,CURVE
C
C LOCAL VARIABLES
DOUBLE PRECISION LBASE(11),NBASE(11)
DOUBLE PRECISION ZERO,TWO,SEVEN,EGHT
DOUBLE PRECISION D1MACH,EP,EPS,FIW,FL,FN,RHOERR,U100
INTEGER I,J,K,II,ID2,ID3,ID4,ID5,ID6,ID7,ID8,ID9,ID10,
+ ID11,IEE,IER1,IER2,K_0,LB,LD1,LD2,LD3,LD4,LD5,
+ LD6,LD7,LD8,LD9,LD10,LD11,LU,M_0,N_0,NB,NFINAL,
+ NFFT,NU
EXTERNAL D1MACH,EVALU,INTEQN
INTRINSIC FLOAT,LOG,MAX,SIGN,SQRT
COMMON /MACHCN/U100
COMMON /BLKEVL/K_0,M_0
COMMON /BLKINT/N_0
DATA ZERO/0.0D0/,TWO/2.0D0/,SEVEN/7.D0/,EGHT/8.0D0/
K_0 = 16
M_0 = 32
N_0 = 16
C
C TEST THE INPUT PARAMETERS.
C
IF (EPS0 .LE. 0) THEN
IER = -1
RETURN
END IF
IF (NWORK .LE. 0 .OR. MWORK .LT.
+ SQRT(FLOAT((NWORK+49)+8)) ) THEN
IER = -2
RETURN
END IF
IF (B .LE. ZERO) THEN
IER = -3
RETURN
END IF
IF ( (IE .LT. 0 .OR. IE .GT. 1) .OR.
+ (IBEG .LT. 0 .OR. IBEG .GT. 1) ) THEN
IER = -4
RETURN
END IF
IF( NPTS .LE. 0) THEN
IER = -5
RETURN
END IF
IF (MAXFFT .LT. M_0) THEN
IER = -6
RETURN
END IF
C SET MACHINE DEPENDENT CONSTANT. U100 IS 100 TIMES
C MACHINE UNIT ROUND.
U100 = 100*D1MACH(4)
EPS = EPS0
C
IF(IBEG .EQ. 1) GO TO 10
EP = EPS/TWO
C OBTAIN VALUE OF NUPPER FOR USE IN INTEQN
FIW = NWORK
FN = SQRT(FIW + SEVEN*SEVEN) - SEVEN
IEE = LOG(FN)/LOG(TWO) + U100
NU = 2**IEE
C BREAK UP WORK INTO VECTORS AND MATRICES FOR USE IN INTEQN.
C PRODUCE RELATIVE ADDRESSES WITHIN WORK OF RHO,X,...,D2Y,OLDRHO.
NBASE(1) = 1
NBASE(2) = 1 + NU*NU
DO I = 3,11
NBASE(I) = NBASE(I-1) + NU
END DO
ID2 = NBASE(2)
ID3 = NBASE(3)
ID4 = NBASE(4)
ID5 = NBASE(5)
ID6 = NBASE(6)
ID7 = NBASE(7)
ID8 = NBASE(8)
ID9 = NBASE(9)
ID10 = NBASE(10)
ID11 = NBASE(11)
C
CALL INTEQN(IOUT,IDBG,IE,CURVE,B,EP,BDYFCN,NU,WORK(ID2),RHOERR,
+ NFINAL,WORK(ID3),WORK(ID4),WORK(ID5),WORK(ID6),
+ WORK(ID7),WORK(ID8),WORK(ID9),WORK(ID10),
+ WORK(ID11),IWORK(16),WORK(1),IER1)
IF (IDBG .EQ. 1) THEN
WRITE (IOUT,9000) NFINAL,RHOERR,IER1
WRITE (IOUT,9010)
DO I = 1,NFINAL
WRITE (IOUT,9020) WORK(ID3+I-1),WORK(ID4+I-1),WORK(ID2+I-1)
END DO
END IF
C
IF (IER1 .EQ. 1) EP = RHOERR
C OBTAIN LUPPER FOR USE IN EVALU.
NFFT =MIN(4*NFINAL,MAXFFT)
FL = (FIW-4.D0*NFFT)/EGHT
IEE = LOG(FL)/LOG(TWO) + U100
LU = 2**IEE
C OBTAIN RELATIVE ADDRESSES FOR BREAKING UP WORK FOR USE
C IN EVALU.
LBASE(1) = 1
DO I = 2,9
LBASE(I) = LBASE(I-1) + LU
END DO
LBASE(10)=LBASE(9) + NFFT
LBASE(11)=LBASE(10) + NFFT
LD1 = LBASE(1)
LD2 = LBASE(2)
LD3 = LBASE(3)
LD4 = LBASE(4)
LD5 = LBASE(5)
LD6 = LBASE(6)
LD7 = LBASE(7)
LD8 = LBASE(8)
LD9 = LBASE(9)
LD10 = LBASE(10)
LD11 = LBASE(11)
C MOVE RHO,X,Y,...,D2Y,SPD AROUND IN WORK, LENGTHENING EACH OF
C THEM.
DO I = 1,8
IF(LBASE(I)+NFINAL-1 .GE. NBASE(I+2)) THEN
DO K = I,8
II = I + 8 - K
LB = LBASE(II) - 1
NB = NBASE(II+1) - 1
DO J = 1,NFINAL
WORK(LB+J) = WORK(NB+J)
END DO
END DO
GO TO 10
END IF
LB = LBASE(I) - 1
NB = NBASE(I+1) - 1
DO J = 1,NFINAL
WORK(LB+J) = WORK(NB+J)
END DO
END DO
C
10 CALL EVALU(IOUT,IDBG,IE,CURVE,B,BDYFCN,NFINAL,WORK(LD1),EP,
+ WORK(LD2),WORK(LD3),WORK(LD4),WORK(LD5),WORK(LD6),
+ WORK(LD7),WORK(LD8),NFFT,WORK(LD9),WORK(LD10),
+ WORK(LD11),IWORK(1),LU,DPTS,NPTS,IBEG,UVEC,ERROR,
+ IER2)
C
IER = MAX(IER1,IER2)
DO I = 1,NPTS
ERROR(I) = ERROR(I) + SIGN(RHOERR,ERROR(I))
END DO
IF(IER1 .EQ. 0) RETURN
DO I = 1,NPTS
IF (ERROR(I) .GT. ZERO) ERROR(I) = -ERROR(I)
END DO
RETURN
9000 FORMAT (/,' FROM NEUMAN: RESULTS FROM SUBROUTINE INTEQN.',
+ /,' NFINAL=',I3,5X,'RHOERR=',1P,E8.2,5X,'IER1=',I1,/)
9010 FORMAT (6X,'X',14X,'Y',19X,'RHO')
9020 FORMAT (1P,D12.4,D15.4,D25.12)
END
C
SUBROUTINE EVALU(IOUT,IDBG,IE,CURVE,B,BDYFCN,N,RHO,EPS,
+ X,Y,DX,DY,D2X,D2Y,SPD,NF,RFFT,OLDFFT,W,
+ IFAC,LU,DPTS,NP,IBEG,U,ERROR,IER)
C ----------------
C
C This program evaluates the single layer potential U at the
C given points in DPTS. The input is the density function RHO,
C and DPTS at which U is evaluated. RHO is evaluated in the
C subroutine INTEQN. These results are stored in U, along
C with the predicted error bound in ERROR. The desired error
C tolerance is EPS. If the desired error bound is not attained,
C then the corresponding entry in ERROR is made negative. Its
C magnitude is still an estimated error bound.
C
C IOUT Input.
C The output unit number.
C IDBG Input.
C The debugging parameter.
C = 0 produces a standard ouput.
C = 1 produces a debugging parameter.
C IE Input.
C =0 for the interior Neumann problem.
C =1 for the exterior Neumann problem.
C CURVE External subroutine.
C This program defines the curve of the boundary C.
C B Input.
C For the parameterization of C defined in CURVE,
C the parameterization interval is [0,B].
C BDYFCN Input.
C This program defines the Neumann data on the
C boundary.
C N Input.
C This is NFINAL as output from the subroutine INTEQ.
C RHO Input.
C An array which contains the value of the single
C layer density function defining U.
C EPS Input.
C The user-supplied absolute error tolerance.
C X,Y Input.
C Two arrays containing a sequence of points
C (X(I),Y(I)) produced by calling the subroutine
C NEWCUR. They correspond to an even subdivision
C of the parameterization interval [0,B].
C DX,DY Input.
C The derivative values corresponding to the points
C given in X,Y.
C D2X,D2Y Input.
C The second derivative values corresponding to the
C points given in X,Y.
C SPD Input.
C An array which contains SQRT(DX*DX+DY*DY).
C NF Input.
C The size of the array RFFT. NF .GE. 2*N.
C RFFT Output.
C An array which contains the discrete Fourier
C coefficients of RHO.
C W Input.
C A work array which is needed in the fft subroutine.
C IFAC Input.
C An integer array which is needed in the fft subroutine.
C LU Input.
C The upper bound on the size of the arrays X,Y,DX,DY,
C D2X,D2Y,SPD,RHO.
C DPTS Input.
C This is a two dimensional array which supplies the
C points at which U is to be evaluated.
C NP Input.
C This is the number of points in DPTS.
C IBEG Input.
C =0 means this is a first call on NEUMAN for this
C particular curve C, boundary function BDYFCN, and
C error tolerance EPS.
C =1 means that NEUMAN has been called previously for
C this choice of parameters, and the solution U is
C desired at a new set of points in DPTS. For such a
C call on NEUMAN, change only DPTS, NPTS, and IBEG.
C Do not change any other input variable.
C U Output.
C An output array which contains U(P) for points P
C given in DPTS.
C ERROR Output.
C An output array which contains the predicted error
C bound for the corresponding entries in U.
C IER Output.
C =0 means the program was completed satisfactorily.
C =1 means some or all of the solution values in U do
C not satisfy the error tolerance EPS.
C
INTEGER IBEG,IE,IER,IOUT,LU,N,NF,NP
DOUBLE PRECISION D2X(LU),D2Y(LU),DPTS(4,NP),DX(LU),DY(LU),
+ ERROR(NP),OLDFFT(NF),RFFT(NF),RHO(LU),SPD(LU),
+ U(NP),W(2*NF),X(LU),Y(LU)
INTEGER IFAC(15)
DOUBLE PRECISION D1MACH,B,BDYFCN,EPS
EXTERNAL BDYFCN,CURVE
C
C LOCAL VARIABLES
DOUBLE PRECISION BDCON,DIFF,DIST,ERR,ERRFFT,FCNK,FCNM,FCNU,FUEVAL,
+ H,HH,OLDIFF,OLDU,ORHO,PASTRT,PI,PX,PY,
+ R,RATIO,RATE,RTLOW,RTUP,S,SINE,COEF,COSN,SING,
+ SUM,T1,T2,THETA,TT,U100
DOUBLE PRECISION ZERO,ONE,TWO,FOUR
INTEGER I,IBD,J,K,KH,KK,KSTEP,K_0,L,LB,LD,LDM1,LOOP,M_0,
+ IDBG, JH, JLOOP
EXTERNAL FUCOEF,NEWCUR,FUEVAL
INTRINSIC ABS,LOG,MAX,MIN,SIN,SQRT
COMMON /MACHCN/U100
COMMON /BLKEVL/K_0,M_0
DATA ZERO/0.0D0/,ONE/1.0D0/,FOUR/4.D0/,RTUP/0.5D0/,
+ RTLOW/.1D0/,TWO/2.0D0/
C
C K_0 IS THE INITIAL DEGREE OF THE FOURIER SERIES EXPANSION TO BE
C USED IN THE APPROXIMATION OF THE DENSITY FUNCTION.
C M_0 DENOTES THE INITIAL NUMBER OF INTEGRATION NODES TO BE USED
C IN THE EVALUATION OF THE SINGLE LAYER POTENTIAL. THESE OPERATIONS
C WILL ALSO BE PERFORMED WITH THE PARAMETER 2*K_0 and 2*M_0, SO THAT
C THE VALUE OF NWORK IN THE CALLING PROGRAM NEEDS TO BE SET
C ACCORDINGLY. ALWAYS SET K_0 AND M_0 TO BE A POWER OF 2.
C
C STAGE 1: EVALUATE THE FOURIER COEFFICIENTS OF RHO.
C INITIALIZE. LD IS THE RUNNING DIMENSION OF RHO,X,...,D2Y. LB
C IS THE RUNNING DIMENSION OF THE FOURIER COEFFICIENT VECTOR RFFT.
C LB CAN BE AS BIG AS NF, WHERE NF IS INPUT FROM NEUMAN AND IS
C SET TO BE MIN(4*NFINAL,MAXFFT) WITH MAXFFT SET BY THE USER ON
C CALLING SUBROUTINE NEUMAN.
C
C DATA 'PI'
PI = FOUR*ATAN(ONE)
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,9000)
WRITE(IOUT,9010)
END IF
IF (IBEG .GT. 0) GO TO 40
RATE = RTUP
LD = N
LB = K_0
LOOP = 1
DO I = 1,LD
RHO(I) = RHO(I)*SPD(I)
END DO
PASTRT=RTUP
C EVALUATE THE FOURIER COEFFICIENTS OF RHO*SPD FOR LATER USE.
C FIRST ASSIGN THE INTERPOLATING POINTS TO RFFT. AFTER
C CALLING FUCOEF, RFFT CONTAINS THE THE FOURIER COEFFICIENTS.
10 KSTEP = LD/LB
DO J = 1, LB-1
JH = J*KSTEP
RFFT(J+1) = RHO(JH)
END DO
RFFT(1) = RHO(LD)
CALL FUCOEF(LB,RFFT,W,IFAC)
DIFF = ZERO
IF (LOOP .EQ. 1) GO TO 20
C ESTIMATE THE ERROR IN THE FFT
DO I = 1, LB/2
COEF = ONE/I
IF (I .LT. LB/4) THEN
COSN = RFFT(2*I) -OLDFFT(2*I)
SINE = RFFT(2*I+1) -OLDFFT(2*I+1)
ELSE IF (I. EQ. LB/4) THEN
COSN = RFFT(2*I) - OLDFFT(2*I)/2
SINE = RFFT(2*I+1)
ELSE IF (I .LT. LB/2) THEN
COSN = RFFT(2*I)
SINE = RFFT(2*I+1)
ELSE
COSN = RFFT(LB)/2
SINE = ZERO
END IF
DIFF = DIFF + TWO*COEF*SQRT(COSN*COSN + SINE*SINE)
END DO
IF (LOOP .EQ. 2) GO TO 20
C UPDATE THE RATE OF CONVERGENCE OF THE FFT
RATE = MAX(PASTRT,RTLOW,MIN(RTUP,ABS(DIFF/OLDIFF)))
PASTRT = MIN(RTUP,ABS(DIFF/OLDIFF))
ERRFFT = RATE/(ONE-RATE)*DIFF
IF (IDBG .EQ. 1) THEN
WRITE(IOUT, 9020) LB, LD,DIFF, ERRFFT, RATE
END IF
IF (ERRFFT .LT. EPS/20) GO TO 30
20 OLDIFF = MAX(D1MACH(4), DIFF)
DO I = 1, LB
OLDFFT(I) = RFFT(I)
END DO
LB = 2*LB
LOOP = LOOP + 1
IF ((LB .LE. LD) .AND. (LB .LE. NF)) GO TO 10
C AN INSUFFICIENT NUMBER OF VALUES IN RHO. THUS VALUES OF RHO ON
C A FINER MESH MUST BE CREATED.
LD = 2*LD
IF (LD .GT. NF) THEN
LD = LD/2
LB = LB/2
LOOP = LOOP - 1
GO TO 30
END IF
C THERE IS SUFFICIENT SPACE IN RHO,X,Y,...,D2Y FOR AN
C INCREASED SUB-DIVISION OF (0,B). MOVE OLD VALUES OF
C RHO,...,D2Y TO MAKE ROOM FOR NEW VALUES.
DO J = 2,LD,2
K = LD + 2 - J
KH = K/2
RHO(K) = RHO(KH)
X(K) = X(KH)
Y(K) = Y(KH)
DX(K) = DX(KH)
DY(K) = DY(KH)
D2X(K) = D2X(KH)
D2Y(K) = D2Y(KH)
SPD(K) = SPD(KH)
END DO
C PRODUCE NEW CURVE PARAMETERS FOR FINER SUBDIVISION.
H = B/LD
LDM1 = LD - 1
DO J = 1, LDM1, 2
S = J*H
CALL NEWCUR(IE,CURVE,S,X(J),Y(J),DX(J),DY(J),D2X(J),
+ D2Y(J))
SPD(J) = SQRT(DX(J)*DX(J)+DY(J)*DY(J))
END DO
C PRODUCE NEW VALUES OF RHO.
HH = B/N
KSTEP = LD/N
DO J = 1,LDM1,2
SUM = ZERO
DO K = KSTEP,LD,KSTEP
T1 = X(J) - X(K)
T2 = Y(J) - Y(K)
FCNK = (DY(J)*T1-DX(J)*T2)/ ((T1*T1+T2*T2)*SPD(J))
SUM = SUM + FCNK*RHO(K)
END DO
IF(IE .EQ. 1) THEN
C FOR THE EXTERIOR PROBLEM
BDCON = BDYFCN(J*H)
ELSE
C FOR THE INTERIOR PROBLEM
R = X(J)**2 + Y(J)**2
BDCON = -ONE/R*BDYFCN(J*H)
END IF
ORHO = - (BDCON+HH*SUM)/PI
RHO(J) = ORHO*SPD(J)
END DO
GO TO 10
30 CONTINUE
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,9030)
WRITE(IOUT,9040) RFFT(1)
WRITE(IOUT,9050) (I, RFFT(2*I), RFFT(2*I +1),I=1,LB/2-1)
WRITE(IOUT,9060) LB/2, RFFT(LB)
WRITE(IOUT,9070)
END IF
C
C STAGE2: BEGIN LOOP TO EVALUATE U AT POINTS P IN DPTS.
C
40 IER = 0
DO I = 1,NP
PX = DPTS(1,I)
PY = DPTS(2,I)
IBD = DPTS(4,I)
IF (IBD.NE.0) THETA = DPTS(3,I)
C IF IT IS AN INTERIOR PROBLEM. CHANGE THE DPTS BY USING
C KELVIN TRANSFORMATION. IF DPTS=(0,0), THEN ASSIGN SOME BIG
C NUMBERS FOR (PX,PY).
IF(IE .EQ. 0) THEN
R = PX*PX + PY*PY
IF(R .LT. U100) THEN
R = U100
PX = 1.D0/R
PY = 1.D0/R
ELSE
PX = PX/R
PY = PY/R
END IF
END IF
IF (IDBG .EQ. 1) THEN
WRITE (IOUT,9080) I,DPTS(1,I),DPTS(2,I),THETA,IBD,LD
END IF
C
C NOW BEGIN EVALUATION OF U(PX,PY) USING NUMERICAL INTEGRATION.
C INITIALIZE, AND BEGIN WITH M_0 NODES.
RATE = RTUP
L = M_0
LOOP = 1
JLOOP = 1
PASTRT = RTUP
C CALCULATE NUMERICAL INTEGRAL WITH L SUBDIVISIONS OF (0,B).
50 SUM = ZERO
KSTEP = LD/L
H = B/L
DO K = KSTEP,LD,KSTEP
KK = K/KSTEP
S = KK*H
T1 = X(K) - PX
T2 = Y(K) - PY
DIST = SQRT(T1*T1+T2*T2)
C IF DPTS IS A BOUNDARY POINT(IBD=1), THE INTEGRAND HAS A
C SINGULARITY AT S=THETA. DIVIDE THE INTEGRAND INTO TWO PARTS.
C FCNM IS THE SMOOTH PART, AND SING IS THE SINGULAR PART.
IF(IBD .EQ. 1) THEN
IF(THETA .EQ. S) THEN
TT = SQRT(DX(K)*DX(K)+DY(K)*DY(K))
FCNM = LOG(TT)
ELSE
TT = DIST/ABS(TWO*SIN((THETA-S)/TWO))
FCNM = LOG(TT)
END IF
ELSE
FCNM = LOG(DIST)
END IF
FCNM = -FCNM*RHO(K)
SUM = SUM + FCNM
END DO
IF(IBD .EQ. 1) THEN
SING = PI*FUEVAL(B,THETA,RFFT,LB)
ELSE
SING = ZERO
END IF
FCNU = H*SUM + SING
IF (IDBG. EQ. 1) THEN
WRITE (IOUT,9090) FCNU,L,LD
END IF
IF(LOOP .EQ. 1) GO TO 80
C ESTIMATE ERROR IN FCNU.
DIFF = ABS(FCNU-OLDU)
IF(LOOP .EQ. 2) GO TO 60
C UPDATE RATE OF CONVERGENCE OF NUMERICAL INTEGRATION.
RATIO = ABS(DIFF/OLDIFF)
C RATE=MAX(PASTRT,RTLOW,MIN(RTUP,RATIO))
PASTRT = MIN(RTUP,RATIO)
60 ERR = (RATE/ (ONE-RATE))*DIFF
IF (IDBG. EQ. 1) THEN
WRITE (IOUT,9100) DIFF,ERR,RATE
END IF
IF (ERR .GT. EPS) GO TO 70
C THE VALUE OF FCNU IS SUFFICIENTLY ACCURATE.
U(I) = FCNU
ERROR(I) = ERR
IF (IBD .EQ. 1) THEN
ERROR(I) = ERR + ERRFFT
IF (ERROR(I) .GT. EPS) THEN
ERROR(I) = -ERROR(I)
IER = 1
END IF
END IF
GO TO 100
C FCNU IS NOT SUFFICIENTLY ACCURATE.
C RE-INITIALIZE FOR ANOTHER NUMERICAL INTEGRATION.
70 OLDIFF = MAX(DIFF,D1MACH(4))
80 OLDU = FCNU
LOOP = LOOP + 1
L = 2*L
IF(L .LE. LD) GO TO 50
C NOT SUFFICIENT VALUES IN RHO. THUS VALUES OF RHO ON A FINER
C MESH MUST BE CREATED.
LD = 2*LD
IF(LD .GT. LU) GO TO 90
C THERE IS SUFFICIENT SPACE IN RHO,X,Y,...,D2Y FOR AN
C INCREASED SUB-DIVISION OF (0,B). MOVE OLD VALUES OF
C RHO,...,D2Y TO MAKE ROOM FOR NEW VALUES.
DO J = 2,LD,2
K = LD + 2 - J
KH = K/2
RHO(K) = RHO(KH)
X(K) = X(KH)
Y(K) = Y(KH)
DX(K) = DX(KH)
DY(K) = DY(KH)
D2X(K) = D2X(KH)
D2Y(K) = D2Y(KH)
SPD(K) = SPD(KH)
END DO
C PRODUCE NEW CURVE PARAMETERS FOR FINER SUBDIVISION.
H = B/LD
LDM1 = LD - 1
DO J = 1,LDM1,2
S = J*H
CALL NEWCUR(IE,CURVE,S,X(J),Y(J),DX(J),DY(J),D2X(J),D2Y(J))
SPD(J) = SQRT(DX(J)*DX(J)+DY(J)*DY(J))
END DO
C PRODUCE NEW VALUES OF RHO.
HH = B/N
KSTEP = LD/N
DO J = 1,LDM1,2
SUM = ZERO
DO K = KSTEP,LD,KSTEP
T1 = X(J) - X(K)
T2 = Y(J) - Y(K)
FCNK = (DY(J)*T1-DX(J)*T2)/ ((T1*T1+T2*T2)*SPD(J))
SUM = SUM + FCNK*RHO(K)
END DO
IF(IE .EQ. 1) THEN
C FOR THE EXTERIOR PROBLEM
BDCON = BDYFCN(J*H)
ELSE
C FOR THE INTERIOR PROBLEM
R = X(J)**2 + Y(J)**2
BDCON = -ONE/R*BDYFCN(J*H)
END IF
ORHO = - (BDCON+HH*SUM)/PI
RHO(J) = ORHO*SPD(J)
END DO
GO TO 50
C THE UPPER LIMITS FOR RHO,X,Y,...,D2Y HAVE BEEN REACHED.
C MARK ERROR BOUND ACCORDINGLY AND CONTINUE ONTO NEXT POINT P.
90 ERROR(I) = -ERR
U(I) = FCNU
IER = 1
LD = LD/2
ERROR(I) = - ERR
IF (IBD .EQ. 1) ERROR(I) = - (ERR + ERRFFT)
100 END DO
RETURN
9000 FORMAT (/,' FROM SUBROUTINE EVALU.',/)
9010 FORMAT (/,' STAGE 1: FOURIER COEFFICIENTS.',/)
9020 FORMAT (' LB =', I4, 3X, 'LD = ', I4, 3X, 'DIFF =', 1PD9.2, 3X,
* 'ERROR =', D9.2, 3X, 'RATE =', 1D9.2)
9030 FORMAT (/,10X,'COEFF. OF COSINE',11X,'COEFF. OF SINE')
9040 FORMAT (/,2X,' 0',3X,E20.12)
9050 FORMAT (1X, I3,3X,E20.12,5X,E20.12)
9060 FORMAT (1X, I3,3X,E20.12)
9070 FORMAT (/,'STAGE2: EVALUATE U(P)', /)
9080 FORMAT (/,' I=',I2,3X,'PX=',1PD11.4,3X,'PY=',D11.4,3X,
+ 'THETA=',D11.4,3X,'IBD=',I1,3X,'LD=',I5)
9090 FORMAT (' NUM INT =',1P,E20.12,5X,'L=',I6,5X,'LD=',I6)
9100 FORMAT (5X,'DIFF=',1PE9.2,5X,' ERROR=',D9.2,5X,'RATE=',D11.4)
END
C
SUBROUTINE FUCOEF(N,R,W,IFAC)
C -----------------
C
C This generates approximate Fourier coefficients for RHO,
C with the calculations done with an FFT program.
C
C The input is R(1),R(2),...R(N), which contain RHO at N evenly
C space points in the interval [0,2*pi], with R(i) = i*2*pi/N.
C The output is R(1)...R(2*N), which contains the Fourier
C coefficients. The array W is workspace for the FFT.
C
C Inputs: N, R, W, IFAC
C Outputs: R
C
INTEGER J,N
DOUBLE PRECISION R(N),W(2*N)
INTEGER IFAC(15)
EXTERNAL DRFFTF,DRFFTI
C
CALL DRFFTI(N,W,IFAC)
CALL DRFFTF(N,R,W,IFAC)
DO 20 J = 1,N
R(J) = R(J)/FLOAT(N)
20 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION FUEVAL(B,S,R,N)
C ---------------------------------
C
C Evaluate the singular part of the integral for the single
C layer potential, evaluated at points on the boundary curve C.
C The input are the Fourier coefficients R(1)...R(N) which is
C evaluated in subroutine FUCOEF.
C
C Inputs: B, S, R, N
C Output: FUEVAL
C
INTEGER I,M,N
DOUBLE PRECISION R(N)
DOUBLE PRECISION B,COEF,COSN,S,SINE,PI
INTRINSIC COS,SIN,ATAN
C
PI = 4.D0*ATAN(1.D0)
M = N/2
SINE = 0.D0
COSN = 0.D0
DO 10 I = 1,M-1
COEF = 1.D0/I
COSN = COSN + COEF*R(2*I)*COS(I*S*2*PI/B)
SINE = SINE + COEF*R(2*I+1)*SIN(I*S*2*PI/B)
10 CONTINUE
FUEVAL= 2.D0*(-SINE + COSN) + R(N)*COS(M*S*2*PI/B)/M
RETURN
END
C
SUBROUTINE INTEQN(IOUT,IDBG,IE,CURVE,B,EPS,BDYFCN,NUPPER,RHO,
+ ERROR,NFINAL,X,Y,DX,DY,D2X,D2Y,SPD,OLDRHO,
+ WORK,IWORK,KERMAT,IER)
C -----------------
C
C This program solves the second kind boundary integral equation
C which arises from solving the exterior Neumann problem as a
C single layer potential.
C
C The integral equation is solved using Nystrom's method with
C the rectangular rule as the quadrature rule. The resulting
C linear system is solved directly using LINPACK routines.
C
C The output is the single layer density function RHO. This
C is to be found with such accuracy that the resulting harmonic
C function has an accuracy of EPS.
C
C This routine assumes the boundary C is at least two times
C continuously differentiable. The boundary C is defined by
C the subroutine CURVE.
C
C The present routine calculates with the rectangular rule for
C N=4,8,16,... until a sufficiently accurate value of RHO is
C obtained. This is subject to N .LE. NUPPER, with the latter
C based on the size of the vector WORK supplied by the user in
C calling NEUMAN.
C
C IOUT Input.
C The output unit number, to the file NEUMAN.ANS
C IDBG Input.
C The debugging parameter.
C = 0 produces a standard ouput.
C = 1 produces a debugging parameter.
C IE Input.
C =0 for the interior Neumann problem.
C =1 for the exterior Neumann problem.
C CURVE External subroutine.
C This program defines the curve of the boundary C.
C B Input.
C For the parameterization of C defined in CURVE,
C the parameterization interval is [0,B].
C EPS Input.
C The user-supplied absolute error tolerance.
C BDYFCN External function.
C This program defines the Neumann data on the boundary.
C NUPPER Input.
C This is the upper bound for the size of linear
C system that can be constructed and solved.
C RHO Output.
C An array which contains the value of the single
C layer density function defining U.
C ERROR Output.
C This is the predicted error estimate for RHO.
C NFINAL Output.
C This is the dimension of the final linear system
C constructed in solving for RHO.
C X,Y Output.
C Two arrays containing a sequence of points
c (X(I),Y(I)) produced by calling the subroutine
C NEWCUR. They correspond to an even subdivision
C of the parameterization interval [0,B].
C DX,DY Output.
C The derivative values corresponding to the points
C given in X,Y.
C D2X,D2Y Output.
C The second derivative values corresponding to the
C points given in X,Y.
C SPD Output.
C An array which contains SQRT(DX*DX+DY*DY).
C OLDRHO Output.
C An array containing the preceding value of RHO,
C also produced in this program.
C WORK Output.
C This is a work array used in the LAPACK routine.
C IWORK Output.
C This is an array for pivoting used in the LAPACK routine.
C KERMAT Output.
C This is array contains the linear system associated
C with the Nystrom method.
C IER Output.
C =0 means the program was completed satisfactorily.
C =1 means the desired error uniform error bound EPS
C for the solution RHO was not attained.
C
INTEGER IE,IER,INFO,IOUT,NFINAL,NUPPER
DOUBLE PRECISION D2X(NUPPER),D2Y(NUPPER),DX(NUPPER),DY(NUPPER),
+ KERMAT(NUPPER,NUPPER),OLDRHO(NUPPER),
+ RHO(NUPPER),SPD(NUPPER),
+ X(NUPPER),Y(NUPPER)
DOUBLE PRECISION D1MACH,B,BDYFCN,EPS,ERROR
EXTERNAL BDYFCN,CURVE,NEWCUR
C
C VARIABLES FOR THE LAPACK ROUTINES.
DOUBLE PRECISION WORK(4*NUPPER)
INTEGER IWORK(NUPPER), NRHS, IDBG
C
C LOCAL VARIABLES
DOUBLE PRECISION DIFF,DIST,H,OLDIFF,PI,R,RATE,RCOND,
+ RTLOW,RTUP,SUM,SUMAX,T1,T2
DOUBLE PRECISION ZERO, ONE, TWO, FOUR
INTEGER I,J,JH,N,NM1,N_0
INTRINSIC ABS,LOG,MAX,MIN,SQRT
COMMON /BLKINT/N_0
DATA RTLOW/.1D0/,RTUP/.5D0/,ZERO/0.0D0/,ONE/1.0D0/,TWO/2.0D0/,
+ FOUR/4.D0/, NRHS/1/
C
C THE PARAMETER N_0 GIVES THE INITIAL NUMBER OF QUADRATURE POINTS
C USED IN THE APPROXIMATION OF THE INTEGRAL EQUATION. THE EQUATION
C WILL ALSO BE SOLVED WITH 2*N_0 NODES. IN THE PROGRAM CALLING
C NEUMAN, THE PARAMETER NWORK SHOULD BE SET ACCORDING. ALWAYS SET
C N_0 TO BE A POWER OF 2.
C
C INITIAL CASE, N=N_0. INITIALIZE PARAMETERS.
IF (IDBG .EQ. 1) THEN
WRITE (IOUT,9000)
END IF
C DATA 'PI'
PI = FOUR*ATAN(ONE)
N = N_0
RATE = RTUP
C DEFINE STEPSIZE AND POINTS ON CURVE.
H = B/N
DO I = 1,N
CALL NEWCUR(IE,CURVE,I*H,X(I),Y(I),DX(I),DY(I),D2X(I),
+ D2Y(I))
SPD(I) = SQRT(DX(I)*DX(I)+DY(I)*DY(I))
END DO
GO TO 20
C DEFINE H AND POINTS ON CURVE FOR CONTINUING LOOP ON N.
10 H = B/N
DO I = 2,N,2
J = N + 2 - I
JH = J/2
X(J) = X(JH)
Y(J) = Y(JH)
DX(J) = DX(JH)
DY(J) = DY(JH)
D2X(J) = D2X(JH)
D2Y(J) = D2Y(JH)
SPD(J) = SPD(JH)
END DO
NM1 = N - 1
DO I = 1,NM1,2
CALL NEWCUR(IE,CURVE,I*H,X(I),Y(I),DX(I),DY(I),D2X(I),
+ D2Y(I))
SPD(I) = SQRT(DX(I)*DX(I)+DY(I)*DY(I))
END DO
C SET UP MATRIX EQUATION, AND EVALUATE THE MAXIMUM NORM OF
C SINGLE LAYER POTENTIAL TO BE USED IN ERROR ESTIMATE.
C HERE 'SUM' REPRESENTS THE NORM.
20 SUMAX = ZERO
DO I = 1,N
IF(IE .EQ. 1) THEN
C BOUNDARY CONDITION FOR THE EXTERIOR PROBLEM
RHO(I) = BDYFCN(I*H)
ELSE
C BOUNDARY CONDITION FOR THE INTERIOR PROBLEM
R = X(I)**2 + Y(I)**2
RHO(I) = -ONE/R*BDYFCN(I*H)
END IF
SUM = ZERO
DO J = 1,N
IF(I .EQ. J) THEN
C DEFINE KERNEL FOR T(I) = T(J).
T1 = DX(I)
T2 = DY(I)
DIST = T1*T1 + T2*T2
KERMAT(I,I) = -PI - H* (T1*D2Y(I)-T2*D2X(I))/
+ (TWO*DIST)
ELSE
C DEFINE KERNEL FOR T(I) .NE. T(J)
T1 = X(I) - X(J)
T2 = Y(I) - Y(J)
DIST = T1*T1 + T2*T2
KERMAT(I,J) = -H* (DY(I)*T1-DX(I)*T2)*SPD(J)/
+ (SPD(I)*DIST)
SUM = SUM + ABS(LOG(DIST))*SPD(J)
END IF
SUM = H*SUM
END DO
SUMAX = MAX(SUMAX,SUM)
END DO
C SOLVE LINEAR SYSTEM.
CALL DGESV(N,NRHS,KERMAT,NUPPER,IWORK,RHO,NUPPER,INFO)
CALL DGECON('I',N,KERMAT,NUPPER,SUMAX,RCOND,WORK,
+ IWORK,INFO)
IF(N .EQ. N_0) GO TO 40
C CALCULATE NORM OF RHO-OLDRHO.
DIFF = ZERO
DO I = 2,N,2
DIFF = MAX(DIFF,ABS(RHO(I)-OLDRHO(I/2)))
END DO
IF(N .EQ. 2*N_0) GO TO 30
C MEASURE RATE OF CONVERGENCE.
RATE = DIFF/OLDIFF
RATE = MAX(RTLOW,MIN(RATE,RTUP))
C ESTIMATE ERROR IN RHO.
30 ERROR = (ONE/RCOND)*SUM*DIFF*RATE/ (ONE-RATE)
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,9010) N,DIFF,ERROR
END IF
IF(ERROR .LE. EPS) THEN
C EXIT FOR SUCCESSFUL RETURN.
NFINAL = N
IER = 0
RETURN
ELSE IF (2*N .GT. NUPPER) THEN
C EXIT FOR UNSUCCESSFUL RETURN.
IER = 1
NFINAL = N
RETURN
END IF
C PREPARE FOR ANOTHER LOOP ON N.
OLDIFF = MAX(DIFF,D1MACH(4))
40 DO I = 1,N
OLDRHO(I) = RHO(I)
END DO
N = 2*N
GO TO 10
9000 FORMAT (/,' FROM SUBROUTINE INTEQN')
9010 FORMAT (' N=',I3,5X,'DIFF=',1PD8.2,5X,'ERROR=',D8.2)
END
C
SUBROUTINE NEWCUR(IE,CURVE,S,X,Y,DX,DY,D2X,D2Y)
C -----------------
C
C Define a new curve if the given problem is an interior
C Neumann problem(IE=0), using the Kelvin transformation.
C
C Inputs: IE, S
C External subroutine: CURVE
C Outputs: X, Y, X, Y, DX, DY, D2X, D2Y
C
DOUBLE PRECISION D2X,D2Y,DX,DY,S,X,Y
INTEGER IE
EXTERNAL CURVE, KVTRNF
DOUBLE PRECISION D2TX,D2TY,DTX,DTY,TX,TY
CALL CURVE(S,X,Y,DX,DY,D2X,D2Y)
IF (IE .EQ. 0) THEN
CALL KVTRNF(X,Y,DX,DY,D2X,D2Y,TX,DTX,D2TX)
CALL KVTRNF(Y,X,DY,DX,D2Y,D2X,TY,DTY,D2TY)
X = TX
Y = TY
DX = DTX
DY = DTY
D2X = D2TX
D2Y = D2TY
END IF
RETURN
END
C
SUBROUTINE KVTRNF(X,Y,DX,DY,D2X,D2Y,TX,TDX,TD2X)
C -----------------
C
C Define the Kelvin transformation. If we disposition
C X's and Y's, we have TY, TDY, TD2Y which are the
C transformed values of Y, DY, D2Y.
C
C Inputs: X, Y, DX, DY, D2X, D2Y
C Outputs: TX, TDX, TD2X
C
DOUBLE PRECISION D2X,D2Y,DX,DY,TD2X,TDX,TX,X,Y
DOUBLE PRECISION DIST,T1,T2,XS,YS
XS = X*X
YS = Y*Y
DIST = XS + YS
TX = X/DIST
TDX = (DX* (YS-XS)-2.D0*X*Y*DY)/ (DIST*DIST)
T1 = D2X* (YS*YS-XS*XS) - 2.D0*X* (XS+YS)* (DY*DY+DX*DX+Y*D2Y)
T2 = -4.D0* (X*DX+Y*DY)* (DX* (YS-XS)-2.D0*X*Y*DY)
TD2X = (T1+T2)/DIST**3
RETURN
END