4

我一直在尝试通过构建短程序来学习 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，这是牛顿方法对所调用值的第一次迭代。

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

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

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

4

2 回答 2

11

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

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

于 2013-11-01T04:13:56.223 回答
7

我忍不住用递归的方式重写它并使用自动微分。当然应该真正使用 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]
于 2013-11-01T15:35:49.673 回答