SUBROUTINE DRCHLT(IOUT,IDBG,IBEG,IE,BDYFCN,CURVE,B,DPTS,NPTS,
+ EPS0,R_FORM,WORK,IWORK,NWORK,MWORK,UVEC,
+ ERROR,IER)
C -----------------
C
C This program will solve the interior and exterior Dirichlet
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 INPUT PARAMETERS:
C
C IOUT Input.
C The output unit number.
C IDBG Input.
C The debugging parameter.
C =0 produces a shortcut output.
C =1 produces a debugging output.
C IBEG Input.
C =0 means this is a first call on DRCHLT for this
C particular curve C, boundary function BDYFCN, and
C error tolerance EPS.
C =1 means that DRCHLT 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 DRCHLT, change only DPTS, NPTS, AND IBEG. Do
C not change any other input variable.
C IE Input.
C =0 for the interior Dirichlet problem.
C =1 for the exterior Dirichlet 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 harmonic function U on the
C boundary curve C.
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.
C B Input.
C The limits of the variable S in defining the curve C
C are 0 .LE. S .LE. B
C DPTS User supplied array.
C This is a two dimensional array specifying the points
C at which the harmonic function U is to be evaluated.
C The dimension statement for DPTS is
C DIMENSION DPTS(2,NPTS)
C The point #J is given by
C X(J)=DPTS(1,J), Y(J)=DPTS(2,J)
C NPTS Input.
C This is the number of points in DPTS. It must be
C greater than zero.
C EPS0 Input.
C The desired absolute error tolerance in the solution U.
C R_FORM Input.
C This specifies the way in which the variable RATE
C is to be defined in the routines INTEQN and EVALU.
C =0 means we use the "normal" way to define RATE
C based on estimating the rate of convergence in the
C approximates calculated to date.
C =1 means we use a "conservative" error test
C in which RATE=0.5 and the approximates are assumed to
C have a very slow rate of convergence. Use this for
C more ill-behaved problems and boundaries.
C WORK Real work array
C Temporary work space. This vector should have a
C dimension NWORK of at least 5,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 great many problems to be
C treated very accurately. For more details on
C computing the needed dimension NWORK for WORK, see
C the discussion in the driver program for DRCHLT or
C the accompanying paper.
C IWORK Integer work array
C Temporary integer work space. It will be used as
C pivot and work array in LAPACK.
C NWORK Input.
C Dimension of WORK. NWORK = 300,000 will solve a lot of
C problems very accurately. See the preceding discussion
C of WORK and the general discussion of NWORK given below.
C MWORK Input.
C Dimension of IWORK. It is dependendent of NWORK.
C MWORK .GE. (SQRT(NWORK+36)-6).
C UVEC Output vector.
C On exit, it will contain the approximate value of the
C solution function U(X(I),Y(I)) at the points (X(I),Y(I))
C 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 EPSO < 0.
C =-2 means NWORK < 0 or MWORK < SQRT(NWORK*36) - 6.
C =-3 means B .LE. 0.
C =-4 means IE or IBEG are out of range.
C =-5 means NPTS<=0, an error.
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 double 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
C will solve many problems very accurately. We define NWORK by
C NWORK=MAX(MAXSYS**2+12*MAXSYS,7*MAXMES)
C The program assumes
C MAXSYS .GE. 64, MAXMES .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 in the routines EVALU and INTEQN, respectively. As can be
C noted from the numbers given, these parameters should be
C chosen as powers of 2.
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,MWORK,NPTS,NWORK,R_FORM
DOUBLE PRECISION DPTS(2,NPTS),ERROR(NPTS),UVEC(NPTS),WORK(NWORK)
INTEGER IWORK(MWORK)
DOUBLE PRECISION B,BDYFCN,EPS0
EXTERNAL BDYFCN,CURVE
C
C LOCAL VARIABLES
DOUBLE PRECISION D1MACH,EP,EPS,FIW,FL,FN,RHOERR,U100
DOUBLE PRECISION ZERO,TWO,SIX,SEVEN
INTEGER I,ID2,ID3,ID4,ID5,ID6,ID7,ID8,ID9,ID10,IEE,IER1,IER2,II,
+ J,K,LB,LD1,LD2,LD3,LD4,LD5,LD6,LD7,LU,NB,NFINAL,NU
INTEGER LBASE(7),NBASE(10)
EXTERNAL D1MACH,EVALU,INTEQN
INTRINSIC FLOAT,LOG,MAX0,SIGN,SQRT
DATA ZERO/0.0D0/, TWO/2.0D0/, SIX/6.0D0/, SEVEN/7.0D0/
C
C TEST THE INPUT PARAMETERS.
C
IF (EPS0 .LE. 0) THEN
IER = -1
RETURN
END IF
IF (NWORK .LE. 0 .OR.
+ MWORK .LE. (SQRT(NWORK+SIX*SIX) - SIX)) 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
C
C SET MACHINE DEPENDENT CONSTANT. U100 IS 100 TIMES
C MACHINE UNIT ROUND.
U100 = 100*D1MACH(4)
EPS = EPS0
IF (IBEG .EQ. 1) GO TO 10
EP = EPS/TWO
C OBTAIN VALUE OF NUPPER FOR USE IN INTEQN
FIW = NWORK
FN = SQRT(FIW+SIX*SIX) - SIX
IEE = LOG(FN)/LOG(TWO)
NU = 2**IEE
C BREAK UP WORK INTO VECTORS AND MATRIX FOR USE IN INTEQ.
C PRODUCE RELATIVE ADDRESSES WITHIN WORK OF RHO,X,...,D2Y,
C OLDRHO,Z.
NBASE(1) = 1
NBASE(2) = 1 + NU*NU
DO I = 3,10
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)
C
CALL INTEQN(IOUT,IDBG,IE,B,EP,R_FORM,BDYFCN,CURVE,NU,WORK(ID2),
+ RHOERR,NFINAL,WORK(ID3),WORK(ID4),WORK(ID5),
+ WORK(ID6),WORK(ID7),WORK(ID8),WORK(ID9),IWORK,
+ WORK(ID10),WORK(1),IER1)
IF (IDBG .EQ. 1) THEN
WRITE (IOUT,FMT=9000) NFINAL,RHOERR,IER1
WRITE (IOUT,FMT=9010)
DO I = 1,NFINAL
WRITE (IOUT,FMT=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.
FL = FIW/SEVEN
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,7
LBASE(I) = LBASE(I-1) + LU
END DO
LD1 = LBASE(1)
LD2 = LBASE(2)
LD3 = LBASE(3)
LD4 = LBASE(4)
LD5 = LBASE(5)
LD6 = LBASE(6)
LD7 = LBASE(7)
C MOVE RHO,X,Y,...,D2Y AROUND IN WORK, LENGTHEN EACH OF THEM.
DO I = 1,7
IF (LBASE(I)+NFINAL-1 .GE. NBASE(I+2)) THEN
DO K = I,7
II = I + 7 - 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,IBEG,IE,B,BDYFCN,CURVE,NFINAL,WORK(LD1),
+ EP,R_FORM,WORK(LD2),WORK(LD3),WORK(LD4),WORK(LD5),
+ WORK(LD6),WORK(LD7),LU,DPTS,NPTS,UVEC,ERROR,IER2)
C
IER = MAX0(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 SUBROUTINE DRCHLT: SUBROUTINE INTEQN RESULTS.',
+ /,' 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
SUBROUTINE EVALU(IOUT,IDBG,IBEG,IE,B,BDYFCN,CURVE,N,RHO,EPS,
+ R_FORM,X,Y,DX,DY,D2X,D2Y,LU,DPTS,NP,
+ U,ERROR,IER)
C ----------------
C
C This program evaluates the double 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, to the file DRCHLT.ANS
C IDBG Input.
C The dubugging parameter.
C =0 produces a shortcut output file.
C =1 produces a debugging output.
C IBEG Input.
C =0 means this is a first call on DRCHLT for this
C particular curve C, boundary function BDYFCN, and
C error tolerance EPS.
C =1 means that DRCHLT 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 DRCHLT, change only DPTS, NPTS, and IBEG.
C Do not change any other input variable.
C IE Input.
C =0 for the interior Dirichlet problem.
C =1 for the exterior Dirichlet 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 External subroutine.
C This program defines the Dirichlet 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 double
C layer density function defining U.
C EPS Input.
C The user-supplied absolute error tolerance.
C R_FORM Input.
C This specifies the way in which the variable RATE
C is to be defined.
C If R_FORM=0, we use the "normal" way to define RATE
C based on estimating the rate of convergence in the
C approximates calculated to date.
C If R_FORM=1, then we use a "conservative" error test
C in which RATE=0.5 and the approximates are assumed to
C have a very slow rate of convergence. Use this for
C more ill-behaved problems and boundaries.
C X,Y Inputs.
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 Inputs.
C The derivative values corresponding to the points
C given in X,Y.
C D2X,D2Y Inputs.
C The second derivative values corresponding to the
C points given in X,Y.
C LU Input.
C The upper bound of the size of the arrays X,Y,DX,DY,
C D2X,D2Y,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 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,IDBG,IE,IER,IOUT,LU,N,NP,R_FORM
DOUBLE PRECISION DPTS(2,NP),X(LU),Y(LU),DX(LU),DY(LU),D2X(LU),
+ D2Y(LU),ERROR(NP),RHO(LU),U(NP)
DOUBLE PRECISION D1MACH,B,BDYFCN,EPS
EXTERNAL BDYFCN,CURVE
C
C LOCAL VARIABLES
DOUBLE PRECISION AIT,BIT,CMIN,CUT,D2F,DF,DIFF,DISTPQ,DNORM,ERR,
+ FCNK,FCNM,FCNU,H,MESH,OLDIFF,OLDU,ONE,PASTRT,
+ PI,PROD,PX,PY,Q,QD2X,QD2Y,QDX,QDY,QSPD,QX,QY,R,
+ RATE,RTLOW,RTUP,S,SB,SLOPE,SONE,SUM,SZERO,T1,T2,
+ TQX,TQY,TR,TWO,TX,TY,VALUE,ZERO,PARM
INTEGER I,ITR,J,JL,JM,JMIN,JSTEP,JU,K,KH,KSTEP,L,LD,LDM1,LOOP,
+ LOOP1,LOOP2,M_0,M2
PARAMETER (M_0=32)
C
C M_0/2 DENOTES THE INITIAL NUMBER OF INTEGRATION NODES TO BE USED
C IN THE EVALUATION OF THE DOUBLE LAYER POTENTIAL, AND IT WILL ALSO
C BE PERFORMED WITH M_0 NODES, SO THAT THE VALUE OF NWORK IN THE
C CALLING PROGRAM NEEDS TO BE SET ACCORDINGLY. ALWAYS SET M_0 TO
C BE A POWER OF 2.
C
DOUBLE PRECISION SQUARE(M_0)
EXTERNAL NEWCUR
INTRINSIC ABS,MAX,MIN,MIN0
DATA ZERO/0.0D0/,ONE/1.0D0/,RTUP/0.5D0/,RTLOW/.1D0/,CUT/.01D0/,
+ TWO/2.0D0/
C
C DATA 'PI'
PI=4.D0*ATAN(1.D0)
C
C INITIALIZE. LD IS THE RUNNING DIMENSION OF RHO,X,...,D2Y.
IF (IBEG .EQ. 0) LD = N
IER = 0
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9000)
C BEGIN LOOP TO EVALUATE U AT POINTS P IN DPTS.
DO I = 1,NP
PX = DPTS(1,I)
PY = DPTS(2,I)
C
IF (IDBG .EQ. 1) THEN
WRITE (IOUT,FMT=9010) I,DPTS(1,I),DPTS(2,I),LD
END IF
C
C IF IT IS THE EXTERIOR PROBLEM, CHANGE THE DPTS BY USING
C THE KELVIN TRANSFORMATION.
IF (IE .EQ. 1) THEN
R = PX*PX + PY*PY
PX = PX/R
PY = PY/R
END IF
C
C BEGIN THE CALCULATION OF THE POINT Q OF C WHICH IS CLOSEST TO P.
C SEARCH USING M2 EVENLY SPACED POINTS OF C, GIVEN IN X,Y, AT
C STEPS OF JSTEP. INITIALLY CALCULATE SQUARES OF APPROPRIATE
C ANGELS.
M2 = MIN0(LD,M_0)
JSTEP = LD/M2
DO J = 1,M2
JM = J*JSTEP
T1 = X(JM) - PX
T2 = Y(JM) - PY
SQUARE(J) = T1*T1 + T2*T2
END DO
CMIN = 1.0D50
DO J = 1,M2
IF (SQUARE(J) .LT. CMIN) THEN
CMIN = SQUARE(J)
JMIN = J
END IF
END DO
C THE POINT (X(K),Y(K)),K=JSTEP*JMIN, IS CLOSEST TO AMONG ALL
C M2 POINTS SEARCHED.
JM = JMIN*JSTEP
T1 = X(JM) - PX
T2 = Y(JM) - PY
PROD = ((T1*DX(JM)+T2*DY(JM))**2)/
+ (CMIN* (DX(JM)**2+DY(JM)**2))
IF (PROD .GT. CUT) THEN
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,*) 'STEP 1 FOR FINDING Q FAILED. PROD=',
+ PROD
END IF
GO TO 10
END IF
C THIS POINT IS ACCEPTABLY CLOSE, AND WILL BE CALLED Q.
C EVALUATE RHO AT Q AND SAVE.
VALUE = RHO(JSTEP*JMIN)
Q = B/LD*(JSTEP*JMIN)
C FOR AN EXTERIOR PROBLEM, OBTAIN THE CORRESPONDING POINT ON THE
C GIVEN CURVE USING THE KELVIN TRANSFORMATION.
R = X(JM)*X(JM) + Y(JM)*Y(JM)
TR = ONE
IF (IE .EQ. 1) TR = ONE/R
QX = TR*X(JM)
QY = TR*Y(JM)
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9020) QX,QY,VALUE
GO TO 60
C LOOK MORE CAREFULLY AT THE POINTS IN X,Y FOR A CLOSEST POINT.
C INITIALLY, SELECT AN INTERVAL IN WHICH TO SEARCH.
10 IF (JMIN .EQ. M2) THEN
JMIN = JMIN*JSTEP
JL = JMIN - JSTEP + 1
JU = JMIN
DO J = JL,JU
T1 = X(J) - PX
T2 = Y(J) - PY
DNORM = T1*T1 + T2*T2
IF (DNORM .LT. CMIN) THEN
CMIN = DNORM
JMIN = J
END IF
END DO
JL = 1
JU = JSTEP
ELSE
JMIN = JMIN*JSTEP
JL = JMIN - JSTEP + 1
JU = JMIN + JSTEP - 1
END IF
DO J = JL,JU
T1 = X(J) - PX
T2 = Y(J) - PY
DNORM = T1*T1 + T2*T2
IF (DNORM .LT. CMIN) THEN
CMIN = DNORM
JMIN = J
END IF
END DO
T1 = X(JMIN) - PX
T2 = Y(JMIN) - PY
PROD = ((T1*DX(JMIN)+T2*DY(JMIN))**2)/
+ (CMIN* (DX(JMIN)**2+DY(JMIN)**2))
IF (PROD .GT. CUT) THEN
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,*) 'STEP 2 FOR FINDING Q FAILED.
+ PROD=', PROD
END IF
GO TO 20
END IF
C THIS POINT Q=(X,Y) IS ACCEPTABLY CLOSE TO P. EVALUATE RHO AT Q.
VALUE = RHO(JMIN)
Q = B/LD*JMIN
C FOR AN EXTERIOR PROBLEM, OBTAIN A CORRESPONDING POINTS ON THE
C GIVEN CURVE USING THE KELVIN TRANSFORMATION.
TR = ONE
R = X(JMIN)*X(JMIN) + Y(JMIN)*Y(JMIN)
IF (IE .EQ. 1) TR = ONE/R
QX = TR*X(JMIN)
QY = TR*Y(JMIN)
IF(IDBG .EQ. 1) WRITE (IOUT,FMT=9030) QX,QY,VALUE
GO TO 60
C NO ACCEPTABLE POINT Q FOUND ON C USING THE VECTORS X,Y.
C BEGIN ITERATION METHOD.
20 H = B/LD
SB = JMIN*H
SZERO = SB
LOOP = 0
QX = X(JMIN)
QY = Y(JMIN)
QDX = DX(JMIN)
QDY = DY(JMIN)
QD2X = D2X(JMIN)
QD2Y = D2Y(JMIN)
T1 = QX - PX
T2 = QY - PY
C ITERATION LOOP.
30 DF = TWO* (T1*QDX+T2*QDY)
D2F = TWO* (QDX**2+QDY**2+T1*QD2X+T2*QD2Y)
SONE = SZERO - DF/D2F
LOOP = LOOP + 1
CALL NEWCUR(IE,CURVE,SONE,QX,QY,QDX,QDY,QD2X,QD2Y)
T1 = QX - PX
T2 = QY - PY
PROD = ((T1*QDX+T2*QDY)**2)/ ((T1*T1+T2*T2)*
+ (QDX*QDX+QDY*QDY))
IF (IDBG .EQ. 1) THEN
WRITE(IOUT,*)'STEP 3 FOR FINDING Q. PROD=', PROD
END IF
IF (PROD .LE. CUT) THEN
Q=SONE
GO TO 50
END IF
C THE NEW POINT (QX,QY) IS NOT SUFFICIENTLY CLOSE. CHECK FOR
C POSSIBLE DIVERGENCE OF ITERATION.
IF ((SB-H .GT. SONE) .OR. (SB+H .LT. SONE)) GO TO 40
C CONTINUE ITERATION.
SZERO = SONE
GO TO 30
C PRIMARY ITERATION IS DIVERGING. GO TO A METHOD GUARANTEED
C TO CONVERGE, THE BISECTION METHOD
40 AIT = SB - H
BIT = SB + H
C BEGINNING OF ITERATION LOOP. WE LIMIT THE NUMBER OF ITERATION
C TO A CERTAIN NUMBER( HERE 10)
DO ITR = 1,10
Q = (AIT+BIT)/2.0
LOOP = LOOP + 1
CALL NEWCUR(IE,CURVE,Q,QX,QY,QDX,QDY,QD2X,QD2Y)
T1 = QX - PX
T2 = QY - PY
SLOPE = T1*QDX + T2*QDY
PROD = (SLOPE*SLOPE)/ ((T1*T1+T2*T2)* (QDX*QDX+QDY*QDY))
IF (IDBG .EQ. 1) WRITE(IOUT,*)'STEP 4 FOR FINDING Q.
+ PROD=', PROD
IF (PROD .LE. CUT) THEN
PARM=Q
GO TO 50
END IF
SLOPE = 2.0*SLOPE
IF (SLOPE .LT. ZERO) THEN
AIT = Q
ELSE
BIT = Q
END IF
END DO
C THE ITERATION IS IN A TIGHT LOOP. SOMETHING IS WRONG ABOUT
C CURVE. FOR THE EXTERIOR PROBLEM, OBTAIN A CORRESPONDING POINT
C ON GIVEN CURVE.
IF (IE .EQ. 1) THEN
R = QX*QX + QY*QY
QX = QX/R
QY = QY/R
END IF
IF(IDBG.EQ.1) WRITE (IOUT,FMT=9040) Q,QX,QY
C A SUFFICIENTLY ACCURATE CLOSEST POINT TO P HAS BEEN FOUND, AND
C NOW EVALUATE RHO AT THIS POINT, USING THE NYSTROM INTERPOLATION
C FORMULA.
50 SUM = ZERO
KSTEP = LD/N
DO K = KSTEP,LD,KSTEP
T1 = X(K) - QX
T2 = Y(K) - QY
SUM = SUM + RHO(K)* (DY(K)*T1-DX(K)*T2)/ (T1*T1+T2*T2)
END DO
TQX = QX
TQY = QY
C FOR AN EXTERIOR PROBLEM, USE THE KELVIN TRANSFORMATION TO OBTAIN
C THE CORRESPONDING POINT ON THE GIVEN CURVE.
IF (IE .EQ. 1) THEN
R = QX**2 + QY**2
TQX = QX/R
TQY = QY/R
END IF
VALUE = - (BDYFCN(TQX,TQY)+KSTEP*H*SUM)/PI
C IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9050) TQX,TQY,VALUE,LOOP
C CLOSEST POINT Q AND THE VALUE OF RHO AT Q HAVE BEEN CALCULATED.
C NOW BEGIN EVALUATION OF U(PX,PY) USING NUMERICAL INTEGRATION.
C INITIALIZE, AND BEGIN WITH M_0/2 NODES.
C
60 CALL NEWCUR(IE,CURVE,Q,QX,QY,QDX,QDY,QD2X,QD2Y)
DISTPQ = SQRT((PX-QX)*(PX-QX) + (PY-QY)*(PY-QY))
QSPD = SQRT(QDX*QDX+QDY*QDY)
RATE = RTUP
PASTRT = RTUP
L = M_0/2
LOOP1 = 1
LOOP2 = 1
C CALCULATE NUMERICAL INTEGRAL WITH L SUBDIVISIONS OF (0,B).
70 SUM = ZERO
KSTEP = LD/L
DO K = KSTEP,LD,KSTEP
T1 = X(K) - PX
T2 = Y(K) - PY
FCNM = (DY(K)*T1-DX(K)*T2)/ (T1*T1+T2*T2)
SUM = SUM + FCNM* (RHO(K)-VALUE)
END DO
FCNU = -TWO*PI*VALUE - (B/L)*SUM
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9060) FCNU,L,LD
IF (LOOP1 .EQ. 1) GO TO 100
C ESTIMATE ERROR IN FCNU.
DIFF = ABS(FCNU-OLDU)
IF (LOOP1 .EQ. 2) GO TO 80
C UPDATE RATE OF CONVERGENCE OF NUMERICAL INTEGRATION.
IF(R_FORM .EQ. 0) THEN
C THE FOLLOWING IS A SOPHISTICATED ERROR ESTIMATOR,
C USUALLY REASONABLY ACCURATE.
RATE = MAX(PASTRT,RTLOW,MIN(RTUP,ABS(DIFF/OLDIFF)))
ELSE
C THE FOLLOWING IS A CONSERVATIVE ERROR ESTIMATOR.
RATE = RTUP
END IF
PASTRT = MIN(RTUP,ABS(DIFF/OLDIFF))
80 ERR = (RATE/ (ONE-RATE))*DIFF
IF(IDBG .EQ. 1) WRITE (IOUT,FMT=9070) ERR,RATE
C
IF (ERR .GT. EPS) GO TO 90
IF (ERR .EQ. ZERO) THEN
ERR = D1MACH(4)
OLDIFF = ERR
END IF
C FOR A POINT CLOSE TO THE BOUNDARY, ITERATE TWO MORE TIMES.
MESH = QSPD*B/L
IF ((MESH .GT. DISTPQ) .AND. (LOOP2 .LE. 2)) THEN
LOOP2 = LOOP2 + 1
GO TO 90
END IF
C THE VALUE OF FCNU IS SUFFICIENTLY ACCURATE.
U(I) = FCNU
ERROR(I) = ERR
GO TO 120
C FCNU IS NOT SUFFICIENTLY ACCURATE.
C RE-INITIALIZE FOR ANOTHER NUMERICAL INTEGRATION.
90 OLDIFF = DIFF
IF(OLDIFF .EQ. ZERO) OLDIFF = D1MACH(4)
100 OLDU = FCNU
LOOP1 = LOOP1 + 1
L = 2*L
IF (L .LE. LD) GO TO 70
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 110
C THERE IS SUFFICIENT SPACE IN RHO,X,Y,...,D2Y FOR AN INCREASED SUB-
C DIVISION OF (0,B).
C MOVE OLD VALUES OF 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)
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))
END DO
C PRODUCE NEW VALUES OF RHO.
H = B/N
KSTEP = LD/N
DO J = 1,LDM1,2
SUM = ZERO
DO K = KSTEP,LD,KSTEP
T1 = X(K) - X(J)
T2 = Y(K) - Y(J)
FCNK = (DY(K)*T1-DX(K)*T2)/ (T1*T1+T2*T2)
SUM = SUM + FCNK*RHO(K)
END DO
C FOR AN EXTERIOR PROBLEM, USE THE KELVIN TRANSFORMATION.
TX = X(J)
TY = Y(J)
IF (IE .EQ. 1) THEN
R = X(J)*X(J) + Y(J)*Y(J)
TX = X(J)/R
TY = Y(J)/R
END IF
RHO(J) = - (BDYFCN(TX,TY)+H*SUM)/PI
END DO
GO TO 70
C THE UPPER LIMITS FOR RHO,X,Y,...,D2Y HAVE BEEN REACHED.
C MARK ERROR BOUND ACCORDINGLY AND CONTINUE ONTO NEXT POINT P.
110 ERROR(I) = -ERR
U(I) = FCNU
IER = 1
LD = LD/2
120 END DO
RETURN
9000 FORMAT (/,' FROM SUBROUTINE EVALU.',/)
9010 FORMAT (/,' I=',I3,4X,'PX=',1P,D11.4,3X,'PY=',D11.4,5X,'LD=',I6)
9020 FORMAT (' QSTAGE1. QX=',1P,D11.4,3X,'QY=',D11.4,3X,'RHO=',
+ D20.12)
9030 FORMAT (' QSTAGE2. QX=',1P,D11.4,3X,'QY=',D11.4,3X,'RHO=',
+ D20.12)
9040 FORMAT (' PROD IS NOT CONVERGING TO ZERO IN LOOP BEGINNING AT
+ 112',/,' Q=',1P,D11.4,5X,'QX=',D11.4,5X,'QY=',/)
C9050 FORMAT (' QSTAGE3. QX=',1P,D11.4,3X,'QY=',D11.4,3X,'RHO=',
C + D20.12,3X,'LOOPS=',I1)
9060 FORMAT (' NUM INT =',1P,E20.12,5X,'L=',I6,5X,'LD=',I6)
9070 FORMAT (' ERROR=',1P,D20.12,5X,'RATE=',D11.4)
END
SUBROUTINE INTEQN(IOUT,IDBG,IE,B,EPS,R_FORM,BDYFCN,CURVE,NUPPER,
+ RHO,ERROR,NFINAL,X,Y,DX,DY,D2X,D2Y,
+ OLDRHO,IWORK,WORK,KERMAT,IER)
C -----------------
C
C This program solves the second kind boundary integral equation
C which arises from solving the interior Dirichlet problem as a
C double 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 LAPACK routines.
C
C The output is the double 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 DRCHLT.
C
C IOUT Input.
C The output unit number, to the file DRCHLT.ANS
C IDBG Input.
C The debugging parameter.
C =0 produces a shortcut output.
C =1 produces a debugging output.
C IE Input.
C =0 for the interior Dirichlet problem.
C =1 for the exterior Dirichlet problem.
C CURVE Input.
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 error tolerance.
C R_FORM Input.
C This specifies the way in which the variable RATE
C is to be defined.
C If R_FORM=0, we use the "normal" way to define RATE
C based on estimating the rate of convergence in the
C approximates calculated to date.
C If R_FORM=1, then we use a "conservative" error test
C in which RATE=0.5 and the approximates are assumed to
C have a very slow rate of convergence. Use this for
C more ill-behaved problems and boundaries.
C BDYFCN Input.
C This program defines the Dirichlet data on the
C 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 double
C layer density function defining U.
C ERROR Output.
C An output array which contains the predicted error
C bound for the corresponding entries in U.
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 OLDRHO Output.
C An array containing the preceding value of RHO,
C also produced in this program.
C IWORK Integer work space
C This is an array for pivoting used for the subroutines
C in LAPACK.
C WORK Real work space.
C A work array for the subroutines in LAPACK.
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 some or all of the solution values in U do
C not satisfy the error tolerance EPS.
C
INTEGER IDBG,IE,IER,INFO,IOUT,NFINAL,NUPPER,R_FORM
DOUBLE PRECISION D2X(NUPPER),D2Y(NUPPER),DX(NUPPER),DY(NUPPER),
+ KERMAT(NUPPER,NUPPER),OLDRHO(NUPPER),
+ RHO(NUPPER),X(NUPPER),Y(NUPPER),WORK(4*NUPPER)
INTEGER IWORK(NUPPER)
DOUBLE PRECISION D1MACH,RTLOW,B,BDYFCN,EPS,ERROR
EXTERNAL BDYFCN,CURVE
C
C LOCAL VARIABLES.
DOUBLE PRECISION DIFF,DIST,H,OLDIFF,ONE,PI,R,RATE,RCOND,
+ RTUP,SUM,SUMAX,T1,T2,TWO,TX,TY,ZERO
INTEGER I,J,JH,N,NM1,N_0,NRHS
EXTERNAL NEWCUR
INTRINSIC ABS,MAX
PARAMETER (N_0 = 32)
DATA RTLOW/.1D0/,RTUP/.5D0/,ZERO/0.0D0/,ONE/1.0D0/,TWO/2.0D0/,
+ 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 DATA 'PI'
PI = 4.D0*ATAN(1.D0)
C INITIAL CASE, N=N_0. INITIALIZE PARAMETERS.
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9000)
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))
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)
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))
END DO
C SET UP MATRIX EQUATION.
20 DO I = 1,N
RHO(I) = BDYFCN(X(I),Y(I))
IF (IE .EQ. 1) THEN
C FOR AN EXTERIOR DIRICHLET PROBLEM, EVALUATE THE BOUNDARY DATA
C ON THE ORIGINALLY GIVEN CURVE BY USING THE KELVIN
C TRANSFORMATION.
R = X(I)*X(I) + Y(I)*Y(I)
TX = X(I)/R
TY = Y(I)/R
RHO(I) = BDYFCN(TX,TY)
END IF
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(J) - X(I)
T2 = Y(J) - Y(I)
DIST = T1*T1 + T2*T2
KERMAT(I,J) = -H* (DY(J)*T1-DX(J)*T2)/DIST
END IF
END DO
END DO
IF (N .EQ. N_0) THEN
SUMAX = 1.D0
GO TO 30
END IF
C CALCULATE PI+NORM(INTEGRAL OPERATOR).
SUMAX = ZERO
DO I = 1,N
SUM = ZERO
DO J = 1,N
IF (I .EQ. J) THEN
SUM = SUM + ABS(PI+KERMAT(I,I)) + PI
ELSE
SUM = SUM + ABS(KERMAT(I,J))
END IF
END DO
SUMAX = MAX(SUMAX,SUM)
END DO
30 CALL DGESV(N,NRHS,KERMAT,NUPPER,IWORK,RHO,NUPPER,INFO)
CALL DGECON('I',N,KERMAT,NUPPER,SUMAX,RCOND,WORK,IWORK,INFO)
IF(RCOND .EQ. ZERO) RCOND = D1MACH(4)
IF (IDBG .EQ. 1) WRITE(IOUT,9030) ONE/RCOND
IF (N .EQ. N_0) GO TO 60
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 40
C MEASURE RATE OF CONVERGENCE.
RATE = DIFF/OLDIFF
IF(RATE .GE. ONE) THEN
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9020) N,ERROR,EPS
IF (2*N .GT. NUPPER) THEN
C EXIT FOR UNSUCCESSFUL RETURN.
IER = 1
NFINAL = N
RETURN
ELSE
GO TO 50
END IF
END IF
IF(R_FORM .EQ. 0) THEN
C THE FOLLOWING IS A SOPHISTICATED ERROR ESTIMATOR,
C USUALLY REASONABLY ACCURATE.
RATE=MAX(RTLOW,MIN(RATE,RTUP))
ELSE
C THE FOLLOWING USES A RATE THAT ASSUMES THE RATE OF
C CONVERGENCE IS PROPORTIONAL TO 1/N. IT IS A CONSERVATIVE,
C BUT SAFER, CHOICE.
RATE = RTUP
END IF
C ESTIMATE ERROR IN RHO.
40 ERROR = (ONE/RCOND)*(SUMAX)*DIFF*RATE/ (ONE-RATE)
IF (IDBG .EQ. 1) WRITE (IOUT,FMT=9010) N,RATE,ERROR,DIFF,SUMAX
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.
50 OLDIFF = DIFF
IF(OLDIFF .EQ. ZERO) OLDIFF = D1MACH(4)
60 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,3X,'RATE=',1P,D8.2,3X,'ERROR=',D8.2,3X,'DIFF=',
+ D8.2,3X,'SUMAX=',D8.2)
9020 FORMAT (' N=',I3,3X,'RATE > 1',25X,'ERROR=',1PD8.2,3X,
+ 'EPS=',D8.2)
9030 FORMAT (/,' CONDITION NUMBER = ',1PD8.2)
END
SUBROUTINE KVTRNF(X,Y,DX,DY,D2X,D2Y,T,DT,D2T)
C -----------------
C
C Define the Kelvin transformation.
C
C INPUTS: X, Y, DX, DY, D2X, D2Y
C OUTPUTS: T, DT, D2T
C
DOUBLE PRECISION D2T,D2X,D2Y,DIST,DT,DX,DY,T,X,Y
DIST = X*X + Y*Y
T = X/DIST
DT = (DX* (Y*Y-X*X)-2*X*Y*DY)/ (DIST*DIST)
D2T = D2X* (Y**4-X**4) - 2*X*DIST* (DY*DY+DX*DX+Y*D2Y)
D2T = D2T - 4* (X*DX+Y*DY)* (DX* (Y*Y-X*X)-2*X*Y*DY)
D2T = D2T/DIST**3
RETURN
END
SUBROUTINE NEWCUR(IE,CURVE,S,X,Y,DX,DY,D2X,D2Y)
C -----------------
C
C If IE=0, the resulting curve C will be the same as that
C defined in the subroutine CURVE.
C If IE=1, the resulting curve will be that produced by
C applying the Kelvin transformation to the original
C boundary curve C.
C
C INPUTS: IE, S
C EXTERNAL SUROUTINE: CURVE
C OUTPUTS: X, Y, DX, DY, D2X, D2Y
C
DOUBLE PRECISION D2TX,D2TY,D2X,D2Y,DTX,DTY,DX,DY,S,TX,TY,X,Y
INTEGER IE
EXTERNAL CURVE,KVTRNF
CALL CURVE(S,X,Y,DX,DY,D2X,D2Y)
IF (IE .EQ. 1) 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