Test problem in detail
To test the convolution codes we need a test
problem. Consider a charged sphere of
radius 10 and charge density 1. Gauss’s
law tells us that
4πr2 E(r)=-(4/3)πr3 for r<R
and 4πr2 E(r)=-(4/3)πR3 for r>R
so that
V(r)=1000/(3r3) r>R
V(r)=50-r2/6 r<R
We also have
which is a convolution integral. In fact this is a special kind of convolution in which both terms are functions of |r|.
The Fortran and wpj files are
for\CONVOL.for for\FFT.FOR for\SINFT.for - for\tprob.wpj
The complete set of Fortran and wpj files for this is in for\tprob.zip
C *** CONVOLUTION BASED ON FFT
EXTERNAL EXPDEC,CSPHER
CALL CONVOL(EXPDEC,CSPHER,NPC)
STOP
END
SUBROUTINE CONVOL(FUN1,FUN2,NPC)
PARAMETER (NMAX=4096)
DIMENSION Y1(NMAX+2),Y2(NMAX+2)
CHARACTER*64 NA,STATUS
REAL*8 TPI
EXTERNAL FUN1,FUN2
DATA EPSILON/1.D-5/
DATA TPI/6.283185308D0/
C *** FIND THE BEST SET OF POINTS TO REPRESENT FUN1(X), AND FUN2(X)
C *** NOTE THAT FUN1 AND FUN2 ALREADY CONTAIN A POWER OF R
RMAX=40.
H=RMAX/NMAX
DO I=1,NMAX
R=H*(I-1)
Y1(I)=FUN1(R)/(2*TPI) fun1=expdec=1
Y2(I)=FUN2(R) fun2=cspher=r r<10
ENDDO =0 r>10
CALL SINFT(Y1,NMAX)
C *** MAKE THIS INTO A 3D SPHERICALLY SYMMETRIC TRANSFORM OF FUN/R
AKH=TPI/(2*RMAX)
NA='FKFUN1N.OUT'
CALL MAOPEN(1,NA,STATUS)
NA='FKFUN1A.OUT'
CALL MAOPEN(2,NA,STATUS)
DO I=2,NMAX
AK=(I-1)*AKH
Y1(I)=2*TPI*H*Y1(I)/AK
T1A=1/(AK*AK)
WRITE(1,'(2G16.4)')AK,Y1(I)
WRITE(2,'(2G16.4)')AK,T1A
C *** CORRECTING THE 1/R TRANSFORM
C Y1(I)=T1A
ENDDO
CLOSE(1)
CLOSE(2)
noting that we need to introduce an epsilon and take its limit to zero
in order to recover the 1/r → 1/k2 transform limit, we should not be surprised to see that the constant does not transform with tremendous
accuracy. The final result, however will surpise us. Remember that the transforms are also summation transforms and that their errors cancel to help us in the final back transforms.
CALL
SINFT(Y2,NMAX) 
NA='FKFUN2N.OUT'
CALL MAOPEN(1,NA,STATUS)
NA='FKFUN2A.OUT'
CALL MAOPEN(2,NA,STATUS)
DO I=2,NMAX
AK=(I-1)*AKH
Y2(I)=2*TPI*H*Y2(I)/AK
T2A=2*TPI*(SIN(10*AK)-10*AK*COS(10*AK))/AK**3
WRITE(1,'(2G16.4)')AK,Y2(I)
WRITE(2,'(2G16.4)')AK,T2A
C Y2(I)=T2A
ENDDO
CLOSE(1)
CLOSE(2)
The triangle form is easier to transform than the constant. Note that the range is to nmax/rmax=200 much further than shown here
there is essentially no error at the points transformed, but there is a sertain shortage of points in the (0,1) range -- and this with 4096 total points.
C *** FORM AND PLOT THE PRODUCT FORM
NA='FKPRODN.OUT'
CALL MAOPEN(1,NA,STATUS)
NA='FKPRODA.OUT'
CALL MAOPEN(2,NA,STATUS)
AKH=TPI/(2*RMAX)
DO I=2,NMAX
AK=(I-1)*AKH
Y1(I)=Y1(I)*Y2(I)*AK
WRITE(1,'(2G16.4)')AK,Y1(I)
AK10=AK*10
ANAL=0
ANAL=(2*TPI)*(SIN(AK10)-AK10*COS(AK10))/(AK**4)
WRITE(2,'(2G16.4)')AK,ANAL
C Y1(I)=ANAL
ENDDO
CLOSE(1)
CLOSE(2)
Y1(1)=2*Y1(2)-Y1(3)
CALL SINFT(Y1,NMAX)
C *** PLOT THE CONVOLUTION
NA='FRCONVN.OUT'
CALL MAOPEN(1,NA,STATUS)
NA='FRCONVA.OUT'
CALL MAOPEN(2,NA,STATUS)
DO I=2,NMAX
R=(I-1)*H
Y1(I)=Y1(I)*2*AKH/(R*TPI*TPI)
WRITE(1,'(2G16.4)')R,Y1(I)
IF(R.LT.10D0)THEN
PHI=50-R*R/6
ELSE
PHI=1000/(3*R)
ENDIF
WRITE(2,'(2G16.4)')R,PHI
ENDDO
CLOSE(1)
CLOSE(2)
and we see that the large r range is where the error in the k transform which was at small k appears. Making R larger will give better small k intervals and will shove this error out to larger values.
RETURN
END
FUNCTION EXPDEC(X)
DATA EPSILON/1.D-5/
EXPDEC=1
RETURN
END
FUNCTION CSPHER(X)
CSPHER=0
IF(X.LT.10.D0)CSPHER=X
RETURN
END
C$INCLUDE SINFT
C$INCLUDE MAOPEN