| Madhusudan Singh 2005-05-14, 7:13 pm |
| Hi
I wrote some code for divergence of a complex vector field. As a means to
that end, I needed code for differentiating the mth component of a vector
wrt to the nth coordinate. Basing my code on the NR routine dfridr, I
wrote :
function df_complex(func,m,x,n,h,err) result(dfc) ! Differentiate
the mth component of func with respect to variable n
real(dp), intent(in) :: h
real(dp), dimension(:), intent(in) :: x
real(dp), intent(out) :: err
integer(i4b), intent(in) :: n,m
complex(dp) :: dfc
interface
function func(x,m)
use nrtype
implicit none
real(dp), dimension(:), intent(in) :: x
integer(i4b), intent(in) :: m
complex(dp) :: func
end function func
end interface
integer(i4b), parameter :: NTAB=10
real(dp), parameter ::
CON=1.4_sp,CON2=CON*CON,BIG=huge(x),SAFE=2.0
integer(i4b) :: ierrmin,i,j
real(dp), dimension(:), allocatable :: hh
real(dp) :: errt,fac
complex(dp), dimension(NTAB,NTAB) :: a
if(h.eq.0) then
print*, 'Supply a non-zero value for h'
return
end if
i=size(x,1)
if(n.gt.i) then
print*,'Specified variable is not in the vector space'
return
else
allocate(hh(i))
hh=0.0_dp
end if
hh(n)=h
a(1,1)=(func(x+hh,m)-func(x-hh,m))/(2.0_dp*hh(n))
err=BIG
do i=2,NTAB
hh=hh/CON
a(1,i)=(func(x+hh,m)-func(x-hh,m))/(2.0*hh(n))
fac=CON2
do j=2,i
a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0_dp)
fac=CON2*fac
errt=max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
if (errt.le.err) then
err=errt
dfc=a(j,i)
end if
end do
if(abs(a(i,i)-a(i-1,i-1)).ge.SAFE*err) then
deallocate(hh)
return
end if
end do
deallocate(hh)
end function df_complex
The problem seems to be that the derivative is a very strong function of h.
I have chosen h to be smaller than the smallest length scale for noticeable
change in the vector component, and then also larger h. The results are
very different. That leaves me with an uneasy feeling about the reliability
of the results.
What is the general strategy for calculating reliable derivatives,
numerically ? I thought that this problem was a simple one (proceeding from
the definition - f(x+h)-f(x-h)/2h), but when I implement it, I get quite a
bit of variation in the results as a function of h.
What do you suggest ?
Thanks.
PS : Hard as it may seem, I have never had to deal with numerical
derivatives until now. Integrals are so much nicer !
|