the procedure diff calculates the first, second or third order derivative of a function by using neville's process to extrapolate from a sequence of simple polynomial approximations based on interpolating points distributed symmetrically about x0 (or lying only on one side of x0 should this be necessary). if the specified tolerance is non-zero then the procedure attempts to satisfy this absolute or relative accuracy requirement, while if it is unsuccessful or if the tolerance is set to zero then the result having the minimum achievable estimated error is returned instead.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(diff_func), | intent(inout) | :: | me | |||
| integer, | intent(in) | :: | iord | 1, 2 or 3 specifies that the first, second or third order derivative,respectively, is required. |
||
| real(kind=wp), | intent(in) | :: | x0 | the point at which the derivative of the function is to be calculated. |
||
| real(kind=wp), | intent(in) | :: | xmin | xmin, xmax restrict the interpolating points to lie in [xmin, xmax], which should be the largest interval including x0 in which the function is calculable and continuous. |
||
| real(kind=wp), | intent(in) | :: | xmax | xmin, xmax restrict the interpolating points to lie in [xmin, xmax], which should be the largest interval including x0 in which the function is calculable and continuous. |
||
| real(kind=wp), | intent(in) | :: | eps | denotes the tolerance, either absolute or relative. eps=0 specifies that the error is to be minimised, while eps>0 or eps<0 specifies that the absolute or relative error, respectively, must not exceed abs(eps) if possible. the accuracy requirement should not be made stricter than necessary, since the amount of computation tends to increase as the magnitude of eps decreases, and is particularly high when eps=0. |
||
| real(kind=wp), | intent(in) | :: | accr | denotes that the absolute (accr>0) or relative (accr<0) errors in the computed values of the function are most unlikely to exceed abs(accr), which should be as small as possible. if the user cannot estimate accr with complete confidence, then it should be set to zero. |
||
| real(kind=wp), | intent(out) | :: | deriv | the calculated value of the derivative |
||
| real(kind=wp), | intent(out) | :: | error | an estimated upper bound on the magnitude of the absolute error in the calculated result. it should always be examined, since in extreme case may indicate that there are no correct significant digits in the value returned for derivative. |
||
| integer, | intent(out) | :: | ifail | will have one of the following values on exit: 0 the procedure was successful. 1 the estimated error in the result exceeds the (non-zero) requested error, but the most accurate result possible has been returned. 2 input data incorrect (derivative and error will be undefined). 3 the interval [xmin, xmax] is too small (derivative and error will be undefined). |
subroutine diff(me,iord,x0,xmin,xmax,eps,accr,deriv,error,ifail)
implicit none
class(diff_func),intent(inout) :: me
integer,intent(in) :: iord !! 1, 2 or 3 specifies that the first, second or third order
!! derivative,respectively, is required.
real(wp), intent(in) :: x0 !! the point at which the derivative of the function is to be calculated.
real(wp), intent(in) :: xmin !! xmin, xmax restrict the interpolating points to lie in [xmin, xmax], which
!! should be the largest interval including x0 in which the function is
!! calculable and continuous.
real(wp), intent(in) :: xmax !! xmin, xmax restrict the interpolating points to lie in [xmin, xmax], which
!! should be the largest interval including x0 in which the function is
!! calculable and continuous.
real(wp), intent(in) :: eps !! denotes the tolerance, either absolute or relative. eps=0 specifies that
!! the error is to be minimised, while eps>0 or eps<0 specifies that the
!! absolute or relative error, respectively, must not exceed abs(eps) if
!! possible. the accuracy requirement should not be made stricter than
!! necessary, since the amount of computation tends to increase as
!! the magnitude of eps decreases, and is particularly high when eps=0.
real(wp), intent(in) :: accr !! denotes that the absolute (accr>0) or relative (accr<0) errors in the
!! computed values of the function are most unlikely to exceed abs(accr), which
!! should be as small as possible. if the user cannot estimate accr with
!! complete confidence, then it should be set to zero.
real(wp), intent(out) :: deriv !! the calculated value of the derivative
real(wp), intent(out) :: error !! an estimated upper bound on the magnitude of the absolute error in
!! the calculated result. it should always be examined, since in extreme case
!! may indicate that there are no correct significant digits in the value
!! returned for derivative.
integer, intent(out) :: ifail !! will have one of the following values on exit:
!! *0* the procedure was successful.
!! *1* the estimated error in the result exceeds the (non-zero) requested
!! error, but the most accurate result possible has been returned.
!! *2* input data incorrect (derivative and error will be undefined).
!! *3* the interval [xmin, xmax] is too small (derivative and error will be
!! undefined).
real(wp) :: acc,beta,beta4,h,h0,h1,h2, &
newh1,newh2,heval,hprev,baseh,hacc1,hacc2,nhacc1, &
nhacc2,minh,maxh,maxh1,maxh2,tderiv,f0,twof0,f1,f2,f3,f4,fmax, &
maxfun,pmaxf,df1,deltaf,pdelta,z,zpower,c0f0,c1,c2,c3,dnew,dprev, &
re,te,newerr,temerr,newacc,pacc1,pacc2,facc1,facc2,acc0, &
acc1,acc2,relacc,twoinf,twosup,s, &
d(10),denom(10),e(10),minerr(10),maxf(0:10),save(0:13),storef(-45:45),factor
integer :: i,j,k,n,nmax,method,signh,fcount,init
logical :: ignore(10),contin,saved
real(wp) :: dummy1,dummy2
integer,parameter :: eta = digits(1.0_wp) - 1 !! minimum number of significant binary digits (apart from the
!! sign digit) used to represent the mantissa of real(wp) numbers. it should
!! be decreased by one if the computer truncates rather than rounds.
integer,parameter :: inf = -minexponent(1.0_wp) - 2 !! the largest possible positive integers subject to
!! 2**(-inf) and -2**(-inf) being representable real(wp) numbers.
integer,parameter :: sup = maxexponent(1.0_wp) - 1 !! the largest possible positive integers subject to
!! 2**sup and -2**sup being representable real(wp) numbers.
real(wp),parameter :: sqrt2 = sqrt(2.0_wp) !! \( \sqrt(2) \)
real(wp),parameter :: sqrt3 = sqrt(3.0_wp) !! \( \sqrt(3) \)
if (iord<1 .or. iord>3 .or. xmax<=xmin .or. &
x0>xmax .or. x0<xmin .or. .not. associated(me%f)) then
ifail = 2
else
acc = accr
twoinf = 2.0_wp**(-inf)
twosup = 2.0_wp**sup
factor = 2.0_wp**(real((inf+sup),wp)/30.0_wp)
if (factor < 256.0_wp)factor=256.0_wp
maxh1 = xmax - x0
signh = 1
if (x0-xmin <= maxh1) then
maxh2 = x0 - xmin
else
maxh2 = maxh1
maxh1 = x0 - xmin
signh = -1
end if
relacc = 2.0_wp**(1.0_wp-eta)
maxh1 = (1.0_wp-relacc)*maxh1
maxh2 = (1.0_wp-relacc)*maxh2
s=128.0_wp*twoinf
if (abs(x0) > 128.0_wp*twoinf*2.0_wp**eta) s = abs(x0)*2.0_wp**(-eta)
if (maxh1 < s) then
! interval too small
ifail =3
return
end if
if (acc < 0.0_wp) then
if (-acc > relacc) relacc = -acc
acc = 0.0_wp
end if
! determine the smallest spacing at which the calculated
! function values are unequal near x0.
f0 = me%f(x0)
twof0 = f0 + f0
if (abs(x0) > twoinf*2.0_wp**eta) then
h = abs(x0)*2.0_wp**(-eta)
z = 2.0_wp
else
h = twoinf
z = 64.0_wp
end if
df1 = me%f(x0+signh*h) - f0
do
if (df1 /= 0.0_wp .or. z*h > maxh1) exit
h = z*h
df1 = me%f(x0+signh*h) - f0
if (z /= 2.0_wp) then
if (df1 /= 0.0_wp) then
h = h/z
z = 2.0_wp
df1 = 0.0_wp
else
if (z*h > maxh1) z = 2.0_wp
end if
end if
end do
if (df1 == 0.0_wp) then
! constant function
deriv = 0.0_wp
error = 0.0_wp
ifail = 0
return
end if
if (h > maxh1/128.0_wp) then
! minimum h too large
ifail = 3
return
end if
h = 8.0_wp*h
h1 = signh*h
h0 = h1
h2 = -h1
minh = 2.0_wp**(-min(inf,sup)/iord)
if (minh < h) minh = h
select case (iord)
case(1)
s = 8.0_wp
case(2)
s = 9.0_wp*sqrt3
case(3)
s = 27.0_wp
end select
if (minh > maxh1/s) then
ifail = 3
return
end if
if (minh > maxh2/s .or. maxh2 < 128.0_wp*twoinf) then
method = 1
else
method = 2
end if
! method 1 uses 1-sided formulae, and method 2 symmetric.
! now estimate accuracy of calculated function values.
if (method /= 2 .or. iord == 2) then
if (x0 /= 0.0_wp) then
dummy1 = 0.0_wp
dummy2 = -h1
call me%faccur(dummy1,dummy2,acc0,x0,twoinf,f0,f1)
else
acc0 = 0.0_wp
end if
end if
if (abs(h1) > twosup/128.0_wp) then
hacc1 = twosup
else
hacc1 = 128.0_wp*h1
end if
if (abs(hacc1)/4.0_wp < minh) then
hacc1 = 4.0_wp*signh*minh
else if (abs(hacc1) > maxh1) then
hacc1 = signh*maxh1
end if
f1 = me%f(x0+hacc1)
call me%faccur(hacc1,h1,acc1,x0,twoinf,f0,f1)
if (method == 2) then
hacc2 = -hacc1
if (abs(hacc2) > maxh2) hacc2 = -signh * maxh2
f1 = me%f(x0 + hacc2)
call me%faccur(hacc2,h2,acc2,x0,twoinf,f0,f1)
end if
nmax = 8
if (eta > 36) nmax = 10
n = -1
fcount = 0
deriv = 0.0_wp
error = twosup
init = 3
contin = .true.
do
n = n+1
if (.not. contin) exit
if (init == 3) then
! calculate coefficients for differentiation
! formulae and neville extrapolation algorithm
if (iord == 1) then
beta=2.0_wp
else if (method == 2) then
beta = sqrt2
else
beta = sqrt3
end if
beta4 = beta**4
z = beta
if (method == 2) z = z**2
zpower = 1.0_wp
do k = 1,nmax
zpower = z*zpower
denom(k) = zpower-1
end do
if (method == 2 .and. iord == 1) then
e(1) = 5.0_wp
e(2) = 6.3_wp
do i = 3,nmax
e(i) = 6.81_wp
end do
else if ((method /= 2 .and. iord == 1) .or. &
(method == 2 .and. iord == 2)) then
e(1) = 10.0_wp
e(2) = 16.0_wp
e(3) = 20.36_wp
e(4) = 23.0_wp
e(5) = 24.46_wp
do i = 6,nmax
e(i) = 26.0_wp
end do
if (method == 2 .and. iord == 2) then
do i = 1,nmax
e(i)=2.0_wp*e(i)
end do
end if
else if (method /= 2 .and. iord == 2) then
e(1) = 17.78_wp
e(2) = 30.06_wp
e(3) = 39.66_wp
e(4) = 46.16_wp
e(5) = 50.26_wp
do i = 6,nmax
e(i) = 55.0_wp
end do
else if (method == 2 .and. iord == 3) then
e(1) = 25.97_wp
e(2) = 41.22_wp
e(3) = 50.95_wp
e(4) = 56.4_wp
e(5) = 59.3_wp
do i = 6,nmax
e(i) = 62.0_wp
end do
else
e(1) = 24.5_wp
e(2) = 40.4_wp
e(3) = 52.78_wp
e(4) = 61.2_wp
e(5) = 66.55_wp
do i = 6,nmax
e(i) = 73.0_wp
end do
c0f0 = -twof0/(3.0_wp*beta)
c1 = 3.0_wp/(3.0_wp*beta-1.0_wp)
c2 = -1.0_wp/(3.0_wp*(beta-1.0_wp))
c3 = 1.0_wp/(3.0_wp*beta*(5.0_wp-2.0_wp*beta))
end if
end if
if (init >= 2) then
! initialization of steplengths, accuracy and other parameters
heval = signh*minh
h = heval
baseh = heval
maxh = maxh2
if (method == 1)maxh = maxh1
do k = 1,nmax
minerr(k) = twosup
ignore(k) = .false.
end do
if (method == 1) newacc = acc1
if (method == -1) newacc = acc2
if (method == 2) newacc = (acc1+acc2)/2.0_wp
if (newacc < acc) newacc = acc
if ((method /= 2 .or. iord == 2) .and. newacc < acc0) newacc = acc0
if (method /= -1) then
facc1 = acc1
nhacc1 = hacc1
newh1 = h1
end if
if (method /= 1) then
facc2 = acc2
nhacc2 = hacc2
newh2 = h2
else
facc2 = 0.0_wp
nhacc2 = 0.0_wp
end if
init = 1
j = 0
saved = .false.
end if
! calculate new or initial function values
if (init == 1 .and. (n == 0 .or. iord == 1) .and. &
.not.(method == 2 .and. fcount >= 45)) then
if (method == 2) then
fcount = fcount + 1
f1 = me%f(x0+heval)
storef(fcount) = f1
f2 = me%f(x0-heval)
storef(-fcount) = f2
else
j = j+1
if (j <= fcount) then
f1 = storef(j*method)
else
f1 = me%f(x0+heval)
end if
end if
else
f1 = me%f(x0+heval)
if (method == 2) f2 = me%f(x0-heval)
end if
if (n == 0) then
if (method == 2 .and. iord == 3) then
pdelta = f1-f2
pmaxf = (abs(f1)+abs(f2))/2.0_wp
heval = beta*heval
f1 = me%f(x0+heval)
f2 = me%f(x0-heval)
deltaf = f1-f2
maxfun = (abs(f1)+abs(f2))/2.0_wp
heval = beta*heval
f1 = me%f(x0+heval)
f2 = me%f(x0-heval)
else if (method /= 2 .and. iord >= 2) then
if (iord == 2) then
f3 = f1
else
f4 = f1
heval = beta*heval
f3 = me%f(x0+heval)
end if
heval = beta*heval
f2 = me%f(x0+heval)
heval = beta*heval
f1 = me%f(x0+heval)
end if
end if
! evaluate a new approximation dnew to the derivative
if (n > nmax) then
n = nmax
do i = 1,n
maxf(i-1) = maxf(i)
end do
end if
if (method == 2) then
maxf(n) = (abs(f1)+abs(f2))/2.0_wp
if (iord == 1) then
dnew = (f1-f2)/2.0_wp
else if (iord == 2) then
dnew = f1+f2-twof0
else
dnew = -pdelta
pdelta = deltaf
deltaf = f1-f2
dnew = dnew + 0.5_wp*deltaf
if (maxf(n) < pmaxf) maxf(n) = pmaxf
pmaxf = maxfun
maxfun = (abs(f1)+abs(f2))/2.0_wp
end if
else
maxf(n) = abs(f1)
if (iord == 1) then
dnew = f1-f0
else if (iord == 2) then
dnew = (twof0-3.0_wp*f3+f1)/3.0_wp
if (maxf(n) < abs(f3)) maxf(n) = abs(f3)
f3 = f2
f2 = f1
else
dnew = c3*f1+c2*f2+c1*f4+c0f0
if (maxf(n) < abs(f2)) maxf(n) = abs(f2)
if (maxf(n) < abs(f4)) maxf(n) = abs(f4)
f4 = f3
f3 = f2
f2 = f1
end if
end if
if (abs(h) > 1) then
dnew = dnew/h**iord
else
if (128.0_wp*abs(dnew) > twosup*abs(h)**iord) then
dnew = twosup/128.0_wp
else
dnew = dnew/h**iord
end if
end if
if (init == 0) then
! update estimated accuracy of function values
newacc = acc
if ((method /= 2 .or. iord == 2) .and. newacc < acc0) newacc = acc0
if (method /= -1 .and. abs(nhacc1) <= 1.125_wp*abs(heval)/beta4) then
nhacc1 = heval
pacc1 = facc1
call me%faccur(nhacc1,newh1,facc1,x0,twoinf,f0,f1)
if (facc1 < pacc1) facc1=(3.0_wp*facc1+pacc1)/4.0_wp
end if
if (method /= 1 .and. abs(nhacc2) <= 1.125_wp*abs(heval)/beta4) then
if (method == 2) then
f1 = f2
nhacc2 = -heval
else
nhacc2 = heval
end if
pacc2 = facc2
call me%faccur(nhacc2,newh2,facc2,x0,twoinf,f0,f1)
if (facc2 < pacc2) facc2 = (3.0_wp*facc2+pacc2)/4.0_wp
end if
if (method == 1 .and. newacc < facc1) newacc = facc1
if (method == -1 .and. newacc < facc2) newacc = facc2
if (method == 2 .and. newacc < (facc1+facc2)/2.0_wp) &
newacc = (facc1+facc2)/2.0_wp
end if
! evaluate successive elements of the current row in the neville
! array, estimating and examining the truncation and rounding
! errors in each
contin = n < nmax
hprev = abs(h)
fmax = maxf(n)
if ((method /= 2 .or. iord == 2) .and. fmax < abs(f0)) fmax = abs(f0)
do k = 1,n
dprev = d(k)
d(k) = dnew
dnew = dprev+(dprev-dnew)/denom(k)
te = abs(dnew-d(k))
if (fmax < maxf(n-k)) fmax = maxf(n-k)
hprev = hprev/beta
if (newacc >= relacc*fmax) then
re = newacc*e(k)
else
re = relacc*fmax*e(k)
end if
if (re /= 0.0_wp) then
if (hprev > 1) then
re = re/hprev**iord
else if (2.0_wp*re > twosup*hprev**iord) then
re = twosup/2.0_wp
else
re = re/hprev**iord
end if
end if
newerr = te+re
if (te > re) newerr = 1.25_wp*newerr
if (.not. ignore(k)) then
if ((init == 0 .or. (k == 2 .and. .not.ignore(1))) &
.and. newerr < error) then
deriv = d(k)
error = newerr
end if
if (init == 1 .and. n == 1) then
tderiv = d(1)
temerr = newerr
end if
if (minerr(k) < twosup/4.0_wp) then
s = 4.0_wp*minerr(k)
else
s = twosup
end if
if (te > re .or. newerr > s) then
ignore(k) = .true.
else
contin = .true.
end if
if (newerr < minerr(k)) minerr(k) = newerr
if (init == 1 .and. n == 2 .and. k == 1 .and. .not.ignore(1)) then
if (newerr < temerr) then
tderiv = d(1)
temerr = newerr
end if
if (temerr < error) then
deriv = tderiv
error = temerr
end if
end if
end if
end do
if (n < nmax) d(n+1) = dnew
if (eps < 0.0_wp) then
s = abs(eps*deriv)
else
s = eps
end if
if (error <= s) then
contin = .false.
else if (init == 1 .and. (n == 2 .or. ignore(1))) then
if ((ignore(1) .or. ignore(2)) .and. saved) then
saved = .false.
n = 2
h = beta * save(0)
heval = beta*save(1)
maxf(0) = save(2)
maxf(1) = save(3)
maxf(2) = save(4)
d(1) = save(5)
d(2) = save(6)
d(3) = save(7)
minerr(1) = save(8)
minerr(2) = save(9)
if (method == 2 .and. iord == 3) then
pdelta = save(10)
deltaf = save(11)
pmaxf = save(12)
maxfun = save(13)
else if (method /= 2 .and. iord >= 2) then
f2 = save(10)
f3 = save(11)
if (iord == 3) f4 = save(12)
end if
init = 0
ignore(1) = .false.
ignore(2) = .false.
else if (.not. (ignore(1) .or. ignore(2)) .and. n == 2 &
.and. beta4*factor*abs(heval) <= maxh) then
! save all current values in case of return to current point
saved = .true.
save(0) = h
save(1) = heval
save(2) = maxf(0)
save(3) = maxf(1)
save(4) = maxf(2)
save(5) = d(1)
save(6) = d(2)
save(7) = d(3)
save(8) = minerr(1)
save(9) = minerr (2)
if (method == 2 .and. iord == 3) then
save(10) = pdelta
save(11) = deltaf
save(12) = pmaxf
save(13) = maxfun
else if (method /= 2 .and. iord >= 2) then
save(10) = f2
save(11) = f3
if (iord == 3) save(12) = f4
end if
h = factor*baseh
heval = h
baseh = h
n = -1
else
init = 0
h = beta*h
heval = beta*heval
end if
else if (contin .and. beta*abs(heval) <= maxh) then
h = beta*h
heval = beta*heval
else if (method /= 1) then
contin = .true.
if (method == 2) then
init = 3
method = -1
if (iord /= 2) then
if (x0 /= 0.0_wp) then
dummy1 = 0.0_wp
dummy2 = -h0
call me%faccur(dummy1,dummy2,acc0,x0,twoinf,f0,f1)
else
acc0 = 0.0_wp
end if
end if
else
init = 2
method = 1
end if
n = -1
signh = -signh
else
contin = .false.
end if
end do
if (eps < 0.0_wp) then
s = abs(eps*deriv)
else
s = eps
end if
ifail = 0
if (eps /= 0.0_wp .and. error > s) ifail = 1
end if
end subroutine diff