Fortran code snippets

12

My earlier program also should have specified the kind argument in the cmplx function. Note to self: do not use the real and cmplx functions without a kind argument. A version of the code using double precision shows that with the complex step you can use very small step sizes and not have problems with roundoff error, as opposed to the conventional approach with dydx = (f(x+h)-f(x-h))/(2*h)

program xcomplex_diff
! illustrate numerical differentiation with a complex step approximation
implicit none
integer, parameter :: wp = kind(1.0d0)
real(kind=wp) :: h
integer :: i
write (*,"(3a18)") "h","err_cmplx","err_real"
do i=1,14
   h = 0.1_wp**i
   write (*,"(3f18.14)") h,abs(1.0_wp-[aimag(exp(cmplx(0.0_wp,h,kind=wp))/h),(exp(h)-exp(-h))/(2*h)])
end do
end program xcomplex_diff

output:

                 h         err_cmplx          err_real
  0.10000000000000  0.00166583353172  0.00166750019844
  0.01000000000000  0.00001666658333  0.00001666674999
  0.00100000000000  0.00000016666666  0.00000016666668
  0.00010000000000  0.00000000166667  0.00000000166689
  0.00001000000000  0.00000000001667  0.00000000001210
  0.00000100000000  0.00000000000017  0.00000000002676
  0.00000010000000  0.00000000000000  0.00000000052636
  0.00000001000000  0.00000000000000  0.00000000607747
  0.00000000100000  0.00000000000000  0.00000002722922
  0.00000000010000  0.00000000000000  0.00000008274037
  0.00000000001000  0.00000000000000  0.00000008274037
  0.00000000000100  0.00000000000000  0.00003338943111
  0.00000000000010  0.00000000000000  0.00024416632505
  0.00000000000001  0.00000000000000  0.00079927783736
1 Like