在Haskell中逼近导数时的困惑行为

如何解决在Haskell中逼近导数时的困惑行为

我已经定义了一个类型类<gupdate xmlns="http://www.google.com/update2/response" protocol="1.0"> <app appid="xxxxxxxxxxx"> <updatecheck codebase="https://extensions.site.com/crx/build.crx" version="1.9.6"/> </app> </gupdate> ，该类型类可以由可以对无穷小操作的任何类型实现。 这是一个示例：

Differentiable

我还定义了一个简单的class Fractional a => Differentiable a where dif :: (a -> a) -> (a -> a) difs :: (a -> a) -> [a -> a] difs = iterate dif instance Differentiable Double where dif f x = (f (x + dx) - f(x)) / dx where dx = 0.000001 func :: Double -> Double func = exp 函数来区分。

但是当我在Double -> Double中测试时，会发生这种情况：

ghc

... $ ghci GHCi,version 8.8.4: https://www.haskell.org/ghc/ :? for help Prelude> :l testing [1 of 1] Compiling Main ( testing.hs,interpreted ) Ok,one module loaded. *Main> :t func func :: Double -> Double *Main> derivatives = difs func *Main> :t derivatives derivatives :: [Double -> Double] *Main> terms = map (\f -> f 0) derivatives *Main> :t terms terms :: [Double] *Main> take 5 terms [1.0,1.0000004999621837,1.000088900582341,-222.0446049250313,4.440892098500626e8] *Main> 的n阶导数的近似值为：

e^x|x=0

在给定设置的情况下，一阶和二阶导数是完全合理的近似值，但是突然之间，[1.0,4.440892098500626e8]上的func的三阶导数是... 0！ 如何！？

解决方法

此处使用的方法是一阶精度的 finite difference method 。

Layman的翻译：它可以工作，但从数字上来说还是很垃圾。具体来说，因为它只有一阶精度，所以即使使用精确实数运算，也需要那些很小的步骤才能获得良好的精度。您确实选择了一个较小的步长，这样很好，但是较小的步长会带来另一个问题：舍入错误。您需要将差异f (x+δx) - f x与较小的δx进行比较，这意味着差异较小，而各个值可能较大。总是会产生浮点误差-例如考虑

Prelude> (1 + pi*1e-13) - 1
3.141931159689193e-13

这实际上可能并没有那么严重，但是由于您随后需要除以δx，因此会增加错误。

当您使用高阶导数时，这个问题只会变得更糟/更复杂，因为现在f' x和f' (x+δx)中的每一个都已经（不相同！）增强了错误，因此差异和再次振作无疑是灾难的根源。

解决问题的最简单方法是切换到二阶精确方法，显而易见的是中心差异。然后，您可以将步骤扩大很多，从而在很大程度上避免舍入问题：

Prelude> let dif f x = (f (x + δx) - f(x - δx)) / (2*δx) where δx = 1e-3
Prelude> take 8 $ ($0) <$> iterate dif exp
[1.0,1.0000001666666813,1.0000003333454632,1.0000004990740052,0.9999917560676863,0.9957312752106873,8.673617379884035,7806.255641895632]

您看到现在的前几对导数很好，但是最终它也变得不稳定–在您迭代 any FD方法时会发生这种情况。但这毕竟不是一个好方法：请注意，对 n -th的每个求值都需要对2 n -1的求值。因此，在导数度上的复杂度为指数。

一种近似不透明函数的第 n 次导数的更好方法是为其拟合第 n 次阶多项式，并对其进行符号/自动区分。或者，如果该功能不是不透明的，则可以通过符号/自动区分其自身。

,

tl; dr：dx分母以指数形式迅速变小，这意味着即使分子中的微小误差也被吹散了。

让我们对第一个“差”近似值（三阶导数）进行一些方程式推理。

dif (dif (dif exp))
= { definition of dif }
dif (dif (\x -> (exp (x+dx) - exp x)/dx))
= { definition of dif }
dif (\y -> ((\x -> (exp (x+dx) - exp x)/dx) (y+dx)
          - (\x -> (exp (x+dx) - exp x)/dx) y
           )/dx)
= { questionable algebra }
dif (\y -> (exp (y + 2*dx) - 2*exp (y + dx) + exp y)/dx^2)
= { alpha }
dif (\x -> (exp (x + 2*dx) - 2*exp (x + dx) + exp x)/dx^2)
= { definition of dif and questionable algebra }
\x -> (exp (x + 3*dx) - 3*exp (x + 2*dx) + 3*exp (x + dx) - exp x)/dx^3

希望现在您可以看到我们正在进入的模式：随着我们使用越来越多的导数，分子中的错误变得越来越严重（因为我们正在计算exp距原始点越来越远，x + 3*dx的距离是例如的三倍，而分母对错误的敏感性变得更高（因为我们正在计算dx^n的第n个导数）。通过三阶导数，这两个因素变得站不住脚了：

> exp (3*dx) - 3*exp (2*dx) + 3*exp (dx) - exp 0
-4.440892098500626e-16
> dx^3
9.999999999999999e-19

因此，您可以看到，尽管分子的误差仅为5e-16，但是分母对误差的敏感度很高，以至于您开始看不到答案。