数値微分で定義域を複素数に拡張して

実関数の数値微分を複素関数で

藤井文夫, 田中真人, 佐藤維美『Fortran90/95による有限要素法プログラミング』（https://www.amazon.co.jp/dp/4621087843 ）をパラパラ眺めていたところ、実数関数の数値微分の桁落ちによる誤差を避けるために、定義域を複素数に拡張して、虚軸側から数値微分する方法が載っていて興味を惹かれました。

真似て計算して見ます。

正則な複素関数の微分は複素平面のどの方向から近づいても同じ値に収束するので、正則関数ならこの方法が使えると思います。

計算結果

とりあえず、f(x) = sin(cos(x)) という、あんまり意味ない関数で試してみることにします。sin も cos もべき関数で定義されて正則なのでその合成関数も正則だったと思いますｗ
　微分は df/dx = -sin(x) * cos(cos(x)) となりますので、この値と数値微分値を比較することにします。

x = 0.5 [ f(0.5)=0.769196354841008, f'(0.5)= -0.306358909189995] の地点で微係数を求めます。数値微分は微小値 h に対して：実数 (f(x+h)-f(x-h)) / (2h) , 複素数 (f(x+ih) - f(x-ih)) / i(2h) で計算します。

ここで図の横軸は数値微分のための微小値 h、縦軸は微分の理論値と数値微分値の差の絶対値です。グラフは log-log なので y= x^n 型のべき関数が直線になります。

図中の橙色の線が実数関数での数値微分、青線が複素関数での数値微分です。赤線は h^2 に比例する直線、緑線は 1/h に比例する直線です。

数値微分のオレンジと青は、微小値 h が小さくなるにつれて（図の右側から左側に行くにつれ）理論値との差が小さくなってゆきます。これは数値微分の打ち切り誤差に比例する項で、テイラー展開すればわかるように、h^2 に比例する誤差を持ちます。この理論値通りに h^2（赤線）比例して小さくなっています。

しかし、実数関数微分の方は途中から誤差が増え始めます。これは、f(x+h) と f(x-h) の値が同じくらいになって引き算による桁落ちが起きるためです。この誤差は有効数字の数に比例するので、橙色の線は左側の領域で 1/h（緑線）に比例する傾きを持っています。

一方、複素関数微分の方は、桁落ちを起こすことなくずっと h^2 に比例した誤差を持っています。確かに、複素関数で虚軸方向から数値微分すると、桁落ちによる誤差を避けることができました。

このまま調子に乗って、カットラインのある Log などで計算するか、微分でないけど超関数でも計算して見るかと思いましたが、雪かきで疲れたのでやめにします。

プログラム

数値微分は、x,h と微分される関数を関数引数として持つ関数を、実数と複素数の場合に定義して、それをさらに総称関数として束ねています。

注：Fortran の sin, cos などの組み込み関数は定義域が複素数になっています。（Fortran2003あたりからかな？ちょっと失念）

    module m_mod
      implicit none
      integer, parameter :: kd = kind(0.0d0)
      real(kd), parameter :: pi = 4 * atan(1.0_kd) 
      interface df
        procedure :: df_r, df_c
      end interface  
      interface 
        pure real(kd) function p_rf(x)
          import
          real(kd), intent(in) :: x
        end function p_rf  
        pure complex(kd) function p_cf(z)
          import
          complex(kd), intent(in) :: z
        end function p_cf  
      end interface
    contains  
      pure real(kd) function rfunc(x) 
        real(kd), intent(in) :: x
        rfunc = sin(cos(x))
      end function rfunc   

      pure complex(kd) function zfunc(z) 
        complex(kd), intent(in) :: z
        zfunc = sin(cos(z))
      end function zfunc   

      pure elemental real(kd) function dfunc(x)
        real(kd), intent(in) :: x
        dfunc = -sin(x) * cos(cos(x))
      end function dfunc

      pure real(kd) function df_r(x, h, func) 
        real(kd), intent(in) :: x, h
        procedure(p_rf) :: func
        df_r = (func(x + h) - func(x - h)) / (2 * h)
      end function df_r

      pure real(kd) function df_c(x, h, func) 
        real(kd), intent(in) :: x, h
        procedure(p_cf) :: func
        complex(kd) :: c1, c0
        c1 = cmplx(x, h, kind = kd)
        c0 = cmplx(x,-h, kind = kd)
        df_c = real( (func(c1) - func(c0)) / cmplx(0.0_kd, 2 * h, kind = kd) )
      end function df_c
    end module m_mod

    program deriv
      use m_mod
      implicit none
      integer, parameter :: n = 30
      real(kd), parameter :: hh = 0.1_kd
      real(kd) :: h(n), x, y1(n), y2(n), d 
      integer :: i
      h = [(hh/2**i, i = 1, n)]
      x = 0.5_kd
      forall(i = 1:n) y1(i) = df(x, h(i), rfunc) 
      forall(i = 1:n) y2(i) = df(x, h(i), zfunc) 
      d = dfunc(x)
      print *, " x f(x) f'(x)", x, rfunc(x), d
      do i = 1, n
        print *, 2 * h(i), d - y1(i), d - y2(i)
      end do
    end program deriv

計算値

  x f(x) f'(x)  0.500000000000000       0.769196354841008
 -0.306358909189995
  0.100000000000000       2.471233679433027E-004 -2.479682338492051E-004
  5.000000000000000E-002  6.186000401675606E-005 -6.191280813439004E-005
  2.500000000000000E-002  1.546995069579005E-005 -1.547325095774266E-005
  1.250000000000000E-002  3.867797066847700E-006 -3.868003329576819E-006
  6.250000000000000E-003  9.669686005242539E-007 -9.669814948209954E-007
  3.125000000000000E-003  2.417433517809542E-007 -2.417441650748309E-007
  1.562500000000000E-003  6.043586719961525E-008 -6.043596573190868E-008
  7.812500000000000E-004  1.510899605428051E-008 -1.510898672840710E-008
  3.906250000000000E-004  3.777260504378432E-009 -3.777246349034868E-009
  1.953125000000000E-004  9.444687254500650E-010 -9.443116288920805E-010
  9.765625000000001E-005  2.356312922557890E-010 -2.360779349785957E-010
  4.882812500000000E-005  5.941669378728420E-011 -5.901945598907332E-011
  2.441406250000000E-005  1.621569545307011E-011 -1.475486399726833E-011
  1.220703125000000E-005  2.573274926476188E-012 -3.688771510468314E-012
  6.103515625000000E-006 -6.521672091253095E-012 -9.222067554048863E-013
  3.051757812500000E-006 -6.521672091253095E-012 -2.305378110634138E-013
  1.525878906250000E-006  2.985811597966403E-011 -5.762057497804562E-014
  7.629394531250000E-007 -4.290146016217022E-011 -1.443289932012704E-014
  3.814697265625000E-007  1.026176921214983E-010 -3.552713678800501E-015
  1.907348632812500E-007 -1.884206124458387E-010 -9.436895709313831E-016
  9.536743164062501E-008 -1.352573830715187E-009 -2.220446049250313E-016
  4.768371582031250E-008  9.757326058235094E-010 -5.551115123125783E-017
  2.384185791015625E-008 -3.680880267253883E-009  0.000000000000000E+000
  1.192092895507813E-008 -3.680880267253883E-009  0.000000000000000E+000
  5.960464477539063E-009  5.632345478900902E-009  0.000000000000000E+000
  2.980232238769531E-009 -1.299410601340867E-008  0.000000000000000E+000
  1.490116119384766E-009  2.425879697121047E-008  0.000000000000000E+000
  7.450580596923829E-010  2.425879697121047E-008  0.000000000000000E+000
  3.725290298461914E-010  2.425879697121047E-008  0.000000000000000E+000
  1.862645149230957E-010  2.425879697121047E-008  0.000000000000000E+000

個人 memo 帳 お絵かき命令

title("derivative of sin(cos(x))")
xlabel("h")
ylabel("|difference|")
loglog(x[1:25],y3[1:25], x[1:25], y4[1:25], x[1:25], 1.0e-16 ./ x[1:25], x[1:25], x[1:25].^2)
savefig("deriv.png")