!==================================================================== ! Driver program: Calculates first- or second-order derivatives ! for a function defined in ni base points ! Method: divided differences for equally spaced points ! AG: October 2009 !==================================================================== implicit none integer, parameter :: ni=101 ! initial arrays size for base points integer, parameter :: nc=21 ! number of points where derivatives to be calc. double precision xmin, xmax ! interval double precision xi(ni), yi(ni) double precision x, yx, yxx, step, f, deriv4 double precision fx, fxx integer i, j ! open files (if needed) !open (unit=1, file='table01.dat') !open (unit=2, file='table02.dat') xmin = 0.0 xmax = 3.1415926 ! ! step 1: generate xi and yi from f(x) xmin, xmax, nint ! step = (xmax-xmin)/float(ni-1) do i=1,ni xi(i) = xmin + step*float(i-1) yi(i) = f(xi(i)) ! write (*,200) xi(i), yi(i) end do ! ! step 2: calculate derivatives ! step = (xmax-xmin)/float(nc-1) write (*,202) do i=1, nc x = xmin + step*float(i-1) yx = cos(x) ! f' fx = deriv4(x, xi, yi, ni, 4, 1) yxx = -sin(x) ! f" fxx= deriv4(x, xi, yi, ni, 4, 2) write (*,201) x, yx, fx, yxx, fxx end do 200 format (2f12.5) 201 format (7f12.5) 202 format (' x fx fx_num', & ' fxx fxx_num') stop end ! ! Function f(x) ! function f(x) implicit none double precision f, x f = sin(x) return end function deriv4(xx, xi, yi, ni, n, m) !==================================================================== ! Evaluate first- or second-order derivatives ! on three or four data points ! using interpolation based on divided differences ! written by: Alex Godunov (October 2009) !-------------------------------------------------------------------- ! input ... ! xx - the abscissa at which the derivative is to be evaluated ! xi() - the arrays of data abscissas ! yi() - the arrays of data ordinates ! ni - size of the arrays xi() and yi() ! n - number of points to evaluate derivatives ! m - order of a derivative (1 or 2) ! output ... ! deriv4 - the first- or second-order derivative !============================================================================*/ implicit none double precision deriv4, xx integer ni, n, m double precision xi(ni), yi(ni) double precision x(n), f(n) double precision d1, d2, d3, h, s integer i, j, k, ix ! exit if too high-order derivative was needed, ! or too many base points were needed for derivatives if (m > 2 .or. m>= n .or. n>5) then deriv4 = 0.0 return end if ! if x is ouside the xi(1)-xi(ni) interval set deriv4=0.0 if (xx < xi(1) .or. xx > xi(ni)) then deriv4 = 0.0 return end if ! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i) i = 1 j = ni do while (j > i+1) k = (i+j)/2 if (xx < xi(k)) then j = k else i = k end if end do ! shift i that will correspond to n-th order of interpolation ! the search point will be in the middle in x_i, x_i+1, x_i+2 ... i = i + 1 - n/2 ! check boundaries: if i is ouside of the range [1, ... n] -> shift i if (i < 1) i=1 if (i + n > ni) i=ni-n+1 ! old output to test i ! write(*,100) xx, i ! 100 format (f10.5, I5) ! just want to use index i ix = i ! initialization of f(n) and x(n) do i=1,n f(i) = yi(ix+i-1) x(i) = xi(ix+i-1) end do ! calculate divided difference coefficients d1 = f(2) - f(1) if (n > 2) d2 = f(3) - 2.0*f(2) + f(1) if (n > 3) d3 = f(4) - 3.0*f(3) + 3.0*f(2) - f(1) h = x(2) - x(1) s = (xx - x(1))/h ! calculate the first order derivative if (m == 1) then deriv4 = (1.0/h)*d1 if (n > 2) deriv4 = deriv4 + (1.0/h)*((2.0*s-1.0)/2.0)*d2 if (n > 3) deriv4 = deriv4 + (1.0/h)*((3.0*s*s-6.0*s+2.0)/6.0)*d3 end if ! calculate the second order derivative if (m == 2 .and. n > 2) then deriv4 = (1.0/h**2)*d2 if (n > 3) deriv4 = deriv4 + (1.0/h**2)*(s-1.0)*d3 end if end function deriv4