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
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
as
x² − y²
(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