递归 Haskell；牛顿法：为什么不收敛？

Question

我一直在尝试通过构建短程序来学习 Haskell。我对函数式编程领域有些陌生，但已经阅读了大量内容。

我在 Haskell 中有一个相对较短的递归函数，用于使用牛顿法找到函数的根，达到浮点数允许的精度：

newtonsMethod :: (Ord a, Num a, Fractional a)  => (a -> a) -> (a -> a) -> a -> a
newtonsMethod f f' x
    | f x < epsilon = x
    | otherwise = 
        newtonsMethod f f' (x - (f x / f' x))
    where
        epsilon = last . map (subtract 1) . takeWhile (/= 1) 
            . map (+ 1) . iterate (/2) $ 1

当我在 GHCi 中解释并插入 newtonsMethod (\ x -> cos x + 0.2) (\ x -> -1 * sin x) (-1) 时，我得到 -1.8797716370899549，这是牛顿方法对所调用值的第一次迭代。

我的第一个问题很简单：为什么它只递归一次？如果您发现此代码的结构或明显错误有任何潜在的改进，也请告诉我。

我的第二个问题，涉及更多一点，是这样的：是否有一些干净的方法来测试这个函数的父调用，看看它是否无法收敛，并相应地退出？

提前感谢您提供的任何答案！

Answer 1

它只运行一次，因为-1.8... 是小于epsilon，一个严格的正数。您想检查差异的绝对值是否在公差范围内。

为此类代码获取收敛诊断的一种方法是将结果生成为惰性列表，这与使用 iterate 找到 epsilon 的方式不同。这意味着您可以通过遍历列表来获得最终结果，但您也可以在导致它的结果的上下文中看到它。

Answer 2

我忍不住以递归方式重新编写它并使用自动微分。当然应该真正使用 AD 包：http://hackage.haskell.org/package/ad。这样就不用自己计算导数了，就可以看到方法收敛了。

data Dual = Dual Double Double
  deriving (Eq, Ord, Show)

constD :: Double -> Dual
constD x = Dual x 0

idD :: Double -> Dual
idD x = Dual x 1.0

instance Num Dual where
  fromInteger n             = constD $ fromInteger n
  (Dual x x') + (Dual y y') = Dual (x + y) (x' + y')
  (Dual x x') * (Dual y y') = Dual (x * y) (x * y' + y * x')
  negate (Dual x x')        = Dual (negate x) (negate x')
  signum _                  = undefined
  abs _                     = undefined

instance Fractional Dual where
  fromRational p = constD $ fromRational p
  recip (Dual x x') = Dual (1.0 / x) (-x' / (x * x))

instance Floating Dual where
  pi = constD pi
  exp   (Dual x x') = Dual (exp x)   (x' * exp x)
  log   (Dual x x') = Dual (log x)   (x' / x)
  sqrt  (Dual x x') = Dual (sqrt x)  (x' / (2 * sqrt x))
  sin   (Dual x x') = Dual (sin x)   (x' * cos x)
  cos   (Dual x x') = Dual (cos x)   (x' * (- sin x))
  sinh  (Dual x x') = Dual (sinh x)  (x' * cosh x)
  cosh  (Dual x x') = Dual (cosh x)  (x' * sinh x)
  asin  (Dual x x') = Dual (asin x)  (x' / sqrt (1 - x*x))
  acos  (Dual x x') = Dual (acos x)  (x' / (-sqrt (1 - x*x)))
  atan  (Dual x x') = Dual (atan x)  (x' / (1 + x*x))
  asinh (Dual x x') = Dual (asinh x) (x' / sqrt (1 + x*x))
  acosh (Dual x x') = Dual (acosh x) (x' / (sqrt (x*x - 1)))
  atanh (Dual x x') = Dual (atanh x) (x' / (1 - x*x))

newtonsMethod' :: (Dual -> Dual) -> Double -> [Double]
newtonsMethod' f x = zs
  where
    zs  = x : map g zs
    g y = y - a / b
      where
        Dual a b = f $ idD y

epsilon :: (Eq a, Fractional a) => a
epsilon = last . map (subtract 1) . takeWhile (/= 1) 
          . map (+ 1) . iterate (/2) $ 1

这给出了以下内容

*Main> take 10 $ newtonsMethod' (\x -> cos x + 0.2) (-1)
[-1.0,
-1.8797716370899549,
-1.770515242616871,
-1.7721539749707398,
-1.7721542475852199,
-1.7721542475852274,
-1.7721542475852274,
-1.7721542475852274,
-1.7721542475852274,
-1.7721542475852274]