program test_fdcoef implicit none c..tests the finite difference coefficient generator integer ngrid parameter (ngrid=41) integer i,iseed double precision dx(ngrid),xgrid(ngrid),f(ngrid),df(ngrid), 1 ddf(ngrid),fp(ngrid),fpp(ngrid),diff1,diff2, 2 period,ran1 double precision pi parameter (pi = 3.1415926535897932384d0) c..create a grid with random cell sizes iseed = -14134 dx(1) = 0.0d0 xgrid(1) = 0.0d0 period = 2.0d0*pi/float(ngrid-1) do i=2,ngrid dx(i) = 2.0d0 * period * ran1(iseed) xgrid(i) = xgrid(i-1) + dx(i) end do c..compute a function and its exact first and second derivatives on the grid do i=1,ngrid f(i) = sin(xgrid(i)) df(i) = cos(xgrid(i)) ddf(i) = -f(i) enddo c..compute the numerical second order first derivative call fss002(xgrid,ngrid,f,fp) c..compute the numerical fourth order second derivative call fss004(xgrid,ngrid,f,fpp) c..compare exact and numerical values write(6,05) 'i','x','dx','f','df/dx exact','df/dx numer','error', 1 '2nd exact','2nd numer','error' 05 format(1x,t3,a,t9,a,t23,a,t37,a,t51,a,t65,a,t79,a, 1 t93,a,t107,a,t121,a) do i=1,ngrid diff1 = 0.0d0 if (df(i) .ne. 0.0) diff1 = (df(i) - fp(i))/df(i) diff2 = 0.0d0 if (ddf(i) .ne. 0) diff2 = (ddf(i) - fpp(i))/ddf(i) write(6,10) i,xgrid(i),dx(i),f(i),df(i),fp(i),diff1, 1 ddf(i),fpp(i),diff2 10 format(1x,i4,1p12e14.6) end do c..thats all end subroutine fss002(grid,ngrid,u,du) implicit none save c..this routine computes second order accurate first derivatives c..on an arbitrary grid c..input: c..ngrid = number of points in the grid, c..grid(ngrid) = array of independent values c..u(ngrid) = function values at the grid points c..output: c..du(ngrid) = first derivative values at the grid points c..declare the pass integer ngrid double precision grid(ngrid),u(ngrid),du(ngrid) c..local variables integer m,n parameter (m=2, n=3) integer i double precision coef(n) c..first point; one sided call fdcoef(m,n,grid(1),grid(1),coef) du(1) = coef(1)*u(1) + coef(2)*u(2) + coef(3)*u(3) c..middle of the grid; central differences do i=2,ngrid-1 call fdcoef(m,n,grid(i),grid(i-1),coef) du(i) = coef(1)*u(i-1) + coef(2)*u(i) + coef(3)*u(i+1) enddo c..last point; one sided call fdcoef(m,n,grid(ngrid),grid(ngrid-2),coef) du(ngrid) = coef(1)*u(ngrid-2)+coef(2)*u(ngrid-1)+coef(3)*u(ngrid) return end subroutine fss004(grid,ngrid,u,du) implicit none save c..this routine computes fourth order accurate second derivatives c..on an arbitrary grid c..input: c..ngrid = number of points in the grid, c..grid(ngrid) = array of independent values c..u(ngrid) = function values at the grid points c..output: c..du(ngrid) = first derivative values at the grid points c..declare the pass integer ngrid double precision grid(ngrid),u(ngrid),du(ngrid) c..local variables integer m,n parameter (m=3, n=5) integer i,j double precision coef(n) c..first point; one sided call fdcoef(m,n,grid(1),grid(1),coef) du(1) = 0.0d0 do j=1,5 du(1) = du(1) + coef(j)*u(j) enddo c..second point; one sided call fdcoef(m,n,grid(2),grid(1),coef) du(2) = 0.0d0 do j=1,5 du(2) = du(2) + coef(j)*u(j) enddo c..middle of the grid do i=3,ngrid-2 call fdcoef(m,n,grid(i),grid(i-2),coef) du(i) = 0.0d0 do j=1,5 du(i) = du(i) + coef(j)*u(i + j - 3) enddo enddo c..second to last point; one sided call fdcoef(m,n,grid(ngrid-1),grid(ngrid-4),coef) du(ngrid-1) = 0.0d0 do j=1,5 du(ngrid-1) = du(ngrid-1) + coef(j)*u(ngrid + j - 5) enddo c..last point; one sided call fdcoef(m,n,grid(ngrid),grid(ngrid-4),coef) du(ngrid) = 0.0d0 do j=1,5 du(ngrid) = du(ngrid) + coef(j)*u(ngrid + j - 5) enddo return end subroutine fdcoef(mord,nord,x0,grid,coef) implicit none save c..this routine implements simple recursions for calculating the weights c..of finite difference formulas for any order of derivative and any order c..of accuracy on one-dimensional grids with arbitrary spacing. c..from bengt fornberg's article c..generation of finite difference formulas on arbitrary spaced grids. c..math. comp., 51(184):699-706, 1988. c..input: c..mord = the order of the derivative c..nord = order of accuracy n c..x0 = point at which to evaluate the coefficients c..grid(nord) = array containing the grid c..output: c..coef(nord) = coefficients of the finite difference formula. c..declare the pass integer mord,nord double precision x0,grid(nord),coef(nord) c..local variables integer nu,nn,mm,nmmin,mmax,nmax parameter (mmax=8, nmax=10) double precision weight(mmax,nmax,nmax),c1,c2,c3,c4,pfac c..zero the weight array do nu=1,nord do nn=1,nord do mm=1,mord weight(mm,nn,nu) = 0.0d0 enddo enddo enddo weight(1,1,1) = 1.0d0 c1 = 1.0d0 nmmin = min(nord,mord) do nn = 2,nord c2 = 1.0d0 do nu=1,nn-1 c3 = grid(nn) - grid(nu) c2 = c2 * c3 c4 = 1.0d0/c3 pfac = grid(nn) - x0 weight(1,nn,nu) = c4 * ( pfac * weight(1,nn-1,nu) ) do mm=2,nmmin weight(mm,nn,nu) = c4 * ( pfac * weight(mm,nn-1,nu) 1 - dfloat(mm-1)*weight(mm-1,nn-1,nu) ) enddo enddo pfac = (grid(nn-1) - x0) weight(1,nn,nn) = c1/c2 * (-pfac*weight(1,nn-1,nn-1)) c4 = c1/c2 do mm=2,nmmin weight(mm,nn,nn) = c4 * (dfloat(mm-1)*weight(mm-1,nn-1,nn-1) - 1 pfac*weight(mm,nn-1,nn-1)) enddo c1 = c2 enddo c..load the coefficients do nu = 1,nord coef(nu) = weight(mord,nord,nu) enddo return end double precision function ran1(idum) implicit none save c..minimal random number generator of park and miller with bays and durham c..shuffle and added safeguards. returns a uniform deviate between 0.0 and 1.0 c..(exclusive of the endpoint values). set or reset idum to any negative c..integer value to initialize. idum must not be altered between calls for c..successive deviates in a sequence. c.. c..period is about 10**8, but passes more statistical tests than ran0 c..declare integer ntab parameter (ntab=32) integer idum,k,j,iv(ntab),iy,ia,im,iq,ir,ndiv double precision am,eps,rnmx parameter (ia=16807, im=2147483647, am=1.0d0/im, 1 iq=127773, ir=2836, ndiv=1+(im-1)/ntab, 2 eps=1.0d-13, rnmx=1.0d0-eps) data iv/ntab*0/, iy/0/ c..initialize the shuffle table if (idum .le. 0 .or. iy .eq. 0) then idum = max(-idum,1) do j=ntab+8,1,-1 k = idum/iq idum = ia*(idum -k*iq) - ir*k if (idum .lt. 0) idum = idum + im if (j .le. ntab) iv(j) = idum enddo iy = iv(1) end if c..start here when not initializing; idum by schrage's method k = idum/iq idum = ia*(idum -k*iq) - ir*k if (idum .lt. 0) idum = idum + im j = 1 + iy/ndiv iy = iv(j) iv(j) = idum ran1 = min(am*iy,rnmx) return end