allsol.hs

import Prelude hiding (null)
import Control.Monad (guard)

import Numeric.AD
import Numeric.Interval


allsol :: (RealFloat a, Ord a)
       => (forall b. Floating b => b -> b) -- 根を求める非線形関数
       -> [Interval a]                     -- 探索する区間
       -> [Interval a]                     -- 根が含まれている区間
allsol _ [] = []
allsol f (x:xs) =
  let result = do
        guard $ 0 `member` f x
        let c = midpoint x
            r = 1 / diff f c
            m = 1 - singleton r * diff f x
            k = singleton (c - r * f c) + m * (x - singleton c)
        guard $ not . null $ k `intersection` x
        if k `isSubsetOf` x
           then Just (Right k)
           else Just (Left $ bisect (k `intersection` x))
   in case result of
        Nothing            ->     allsol f xs
        Just (Right k)     -> k : allsol f xs
        Just (Left (y, z)) ->     allsol f (y : z : xs)


inm :: (RealFloat a, Ord a)
    => (forall b. Floating b => b -> b) -- 根を求める非線形関数
    -> Interval a                       -- 区間の初期値
    -> Interval a                       -- 計算結果の区間
inm f x
  | width x < 1e-6 = x
  | otherwise      =
    let c  = midpoint x
        nx = singleton c - (singleton (f c) / diff f x)
     in inm f (x `intersection` nx)


main :: IO ()
main = do
  let f x = (x - 1) * (x - 2) * (x - 3)
  print $ map (inm f) $ allsol f [0...5]