Permalink
Browse files

try implementing a Durand-Kerner polynomial solver; it seems to have …

…fiddly convergence issues to deal with too, plus you have to deal with complex numbers and it's difficult if you only want real solutions...
  • Loading branch information...
1 parent c2ee898 commit eb67bd34562de55b8f201b5bc44fcf4e38f63057 Brent Yorgey committed May 5, 2011
Showing with 51 additions and 0 deletions.
  1. +51 −0 misc/DKSolve.hs
View
@@ -0,0 +1,51 @@
+------------------------------------------------------------
+-- Durand-Kerner method
+------------------------------------------------------------
+
+-- See http://en.wikipedia.org/wiki/Durand–Kerner_method
+
+
+import Data.Complex
+import Data.List (inits, tails)
+
+eps :: Double
+eps = 1e-14
+
+-- | Given as input a list of polynomial coefficients (least
+-- significant first), return a list of the /real/ roots.
+durandKerner :: [Double] -> [Double]
+durandKerner as = map realPart . filter ((<(sqrt eps)) . abs . imagPart) . fixedPt eps (dkIter as) $ initial
+ where initial = take (length as - 1) $ iterate (*(0.4 :+ 0.9)) 1
+
+-- | Given the polynomial coefficients, perform one iteration of the
+-- D-K method.
+dkIter :: [Double] -> [Complex Double] -> [Complex Double]
+dkIter as rs = zipWith (-) rs (zipWith (/) evals denoms)
+ where evals = map (eval as) rs
+ denoms = zipWith mkDenom rs (drops rs) -- (skipZip rs' rs)
+ mkDenom r = product . map ((-) r)
+
+drops :: [a] -> [[a]]
+drops [x] = [[]]
+drops (x:xs) = xs : map (x:) (drops xs)
+
+{-
+skipZip :: [a] -> [a] -> [[a]]
+skipZip xs ys = zipWith (++) (initsL xs) (tail (tails ys))
+
+initsL :: [a] -> [[a]]
+initsL xs = [] : initsL' xs
+ where initsL' [] = []
+ initsL' (x:xs) = map (x:) (initsL xs)
+-}
+
+-- | Evaluate a polynomial for a complex input.
+eval :: [Double] -> Complex Double -> Complex Double
+eval as x = foldr (\a v -> (a :+ 0) + x*v) 0 as
+
+type C = Complex Double
+
+fixedPt :: Double -> ([C] -> [C]) -> [C] -> [C]
+fixedPt eps f as | all ((<eps) . realPart . abs) $ zipWith (-) as as' = as
+ | otherwise = fixedPt eps f as'
+ where as' = f as

0 comments on commit eb67bd3

Please sign in to comment.