The Algorithm

we're concerned with an analytic function. that means the function is infinitely differentiable and can be smoothly extended into the complex plane. computationally, it probably means that it is defined by a single "one line" formula, not a more extensive piece of code with if statements and for loops

let F(z) be such a function, let x0 be a point on the real axis, and let h be a real parameter
expand F(z) in a Taylor series off the real axis:


      F(x0+ih) = F(x0) + ih ⋅ F'(x0) - h2 ⋅ F"(x0)/2! - ih3 ⋅ F'''/3! + .. 
take the imaginary part of both sides and divide by h :

    F'(x0) = Im [ F(x0+ih)) / h + O(h^2)

simply evaluating the function F at the imaginary argument x0+ih, and dividing by h, gives an approximation to the value of the derivative, F'(x0), that is accurate to order o(h^2)

метод комплексного шага

подумалось, что вещи тривиально-обыденные поступают в широкие (и даже научные) массы очень медленно

допустим, нет ничего проще как взять численно производную, верно? "дэ икс по дэ те", переносится в код машинально

естественно, в самых же первых курсах обучают что разностные схемы численно неустойчивы - числа можно слагать, умножать, делить, но с вычитанием надо быть осторожным - можно легко перейти от шести значащих цифр к ни одной в итоге. тем не менее все равно, не только в Джей Пи Моргане, но даже и среди физиков в обычае до сих пор писать конечные разности при взятии простых производных и при решении дифуравнений

между тем, легко и просто делать это правильно, методом комплексного шага :

    ∂F/∂x ≅ Im (F(x + ih)) / h 

вместо "минуса" имеем "плюс" и никаких тебе численных нестабильностей. вот так это выглядит например в octave :

classdef CSD < handle

  properties

    default_step = 1e-9 ;

  end
    
  methods

    function g = gradient ( this , f , x , step )

      num_var = length ( x ) ;
      g = zeros ( num_var , 1 ) ;

      if nargin < 4 dx = this.default_step ; else dx = step ; end

      for n = 1 : num_var
        [ fval ]   =  f ( x ( n ) + 1i * dx ) ;
        g ( n ) =  imag ( fval / dx ) ;
      end

    end

  end    

end 

проиграть исполнение можно посредством следующего:

clear all 

x = linspace ( 0 , 2 * pi , 500 ) ;
f = @(x) sin ( x ) .^ 2 ;

csd = CSD ;

d = csd.gradient ( f , x , 1e-20 ) ;

plot ( sin (x) .^ 2 ) ;
hold on ;

plot ( d , 'r' ) ;
hold off ;

pause ; 

результат работы этого скрипта :

для дифур в частных производных более высоких степеней - чуть менее прозрачно, но тоже применение не вызывает особых сложностей

есть еще пара методов, более быстрых, но этот вполне разумен и исключительно точен

следует помнить - любой, кто использует "конечные разности" для решения даже простых задач, рискует получить вполне заслуженное клеймо "идиота" и у него должны быть наготове очень веские аргументы - почему он делает именно так

ну и само собой такие вещи встречаются сплошь и рядом в экономике, финансах, статистике и т.д.


можно подумать, при вычислении в комплексных вычитать числа уже не придётся. примените этот ваш трюк к F(x) = 1/x и посмотрите, не найдётся ли там, часом, вычитаний

подумать конечно можно и даже нужно. нет, не придется

вообще-то данный трюк работает только тогда, когда функция F - комплексно дифференцируема, т.е. её компоненты - гармонические. таковы далеко не всякие "обычные" функции

верно - но только если подходить к вопросу примитивно "в лоб", по физтеховски

а "непримитивно" - это как? интересно, кто-нибудь там у вас доказывает в каждом конкретном случае, что это работает, или так - "По идее - должно работать, давай применим"?

это небесплатная информация, которая к тому же вам не нужна. достаточно того, что практически всегда это работает, редкие случаи - когда нет - требуют немного подумать головой

ну, конечно, как же иначе! "Вы не правы, а почему - я вам не скажу"

ищите сами, гуглем и по-английски. яндексом и по-русски вы вряд ли найдете

что искать-то?

subtractive cancellation errors


cancellation subtraction

subtraction of two numbers of similar magnitude gives a result that might have a smaller exponent and thus loses significant digits

    1.2345678 · 10⁵ − 1.2345663 · 10⁵ = 1.5 · 10¹  

thus it is a bad idea to compute a small quantity as difference of large quantities

so, as example, it can be better to evaluate

      x² − y² 
as
      (x + y) * (x − y)


спасибо, посмотрю

кроме "конечных разностей", миллион и других методов имеется. что ж всё в каменном веке алгоритмов-то сидеть?

да ну? правда?

правда. во-первых, производную не обязательно только по соседним точкам определять. а во-вторых, есть ведь и метод конечных элементов

этот ваш финт возможен только случае, если функция ф - комплексно дифференцируема, или я заблуждаюсь?

функция изначально вещественная, достаточно чтобы она была кусочно непрерывна

угу. особенно хорошо с F=|x| получается - функцией, комплексный вариант которой не имеет вообще мнимой части


James Delaney убедительно показал, что рекурсия

    I0 = 1 - 1 / e
    In+1 = 1 - n ⋅ In
не работает, из-за потери точности по причине "subtractive cancellation" :
     0     .6321205588285577
     1     .3678794411714423
     2     .2642411176571153
     3     .207276647028654
     4     .170893411885384
     5     .1455329405730801
     6     .1268023565615195
     7     .1123835040693635
     8     .1009319674450921
     9     9.161229299417073D-02
    10     .0838770700582927
    11     7.735222935878028D-02
    12     7.177324769463667D-02
    13     6.694777996972334D-02
    14     6.273108042387321D-02
    15     5.903379364190187D-02
    16     5.545930172957014D-02
    17     5.719187059730757D-02
    18     -2.945367075153627D-02
    19     1.559619744279189
    20     -30.19239488558378
    21     635.0402925972594
    22     -13969.88643713971
    23     321308.3880542133
    24     -7711400.313301118
    25     192785008.832528
    26     -5012410228.645727
    27     135335076174.4346
    28     -3789382132883.17
    29     109892081853612.9
    30     -3296762455608386 

уход в рациональные из флоатов - не помогает. в файле subc.hs (хаскель) :

main :: IO ()
main =
  let (n , d) = fun_rat 30
  in print $ (fromIntegral n) / (fromIntegral d)

type Rat = (Integer , Integer)

toRat :: Integer -> Rat
toRat n = (n , 1)

rev_e_rat :: Rat
rev_e_rat = (1000000000 , 2718281828)

sub_rat :: Rat -> Rat -> Rat
sub_rat (nom1 , denom1) (nom2 , denom2) =
  (nom1 * denom2 - nom2 * denom1 , denom1 * denom2)

mul_rat :: Rat -> Rat -> Rat
mul_rat (nom1 , denom1) (nom2 , denom2) =
  (nom1 * nom2 , denom1 * denom2)

unit_rat :: Rat
unit_rat = (1 , 1)

fun_rat :: Integer -> Rat 
fun_rat 0 = sub_rat unit_rat rev_e_rat
fun_rat n = sub_rat unit_rat $ mul_rat (toRat n) (fun_rat $ n - 1)  

результат выполнения:

λ> :l subc.hs
[1 of 1] Compiling Main             ( subc.hs, interpreted )
Ok, one module loaded.
λ> main
-1.6478838413564181e22