Gaussian elimination: Difference between revisions

From scratch
(From scratch)
Line 1,413:
 
=={{header|Haskell}}==
 
===Version 1===
===From scratch===
<lang Haskell>
 
isMatrix xs = null xs || all ((== (length.head $ xs)).length) xs
 
isSquarredMatrix xs = null xs || all ((== (length xs)).length) xs
 
multiply:: Num a => [[a]] -> [[a]] -> [[a]]
multiply us vs = map (mult [] vs) us
where
mult xs [] _ = xs
mult xs _ [] = xs
mult [] (zs:zss) (y:ys) = mult (map (\v -> v*y) zs) zss ys
mult xs (zs:zss) (y:ys) = mult (zipWith (\u v -> u+v*y) xs zs) zss ys
 
gauss::[[Double]] -> [[Double]] -> [[Double]]
gauss xs bs = map (map fromRational) $ solveGauss (toR xs) (toR bs)
where toR = map $ map toRational
 
solveGauss:: (Fractional a, Ord a) => [[a]] -> [[a]] -> [[a]]
solveGauss xs bs | null xs || null bs || length xs /= length bs || (not $ isSquarredMatrix xs) || (not $ isMatrix bs) = []
| otherwise = uncurry solveTriangle $ triangle xs bs
 
solveTriangle::(Fractional a,Eq a) => [[a]] -> [[a]] -> [[a]]
solveTriangle us _ | not.null.dropWhile ((/= 0).head) $ us = []
solveTriangle ([c]:as) (b:bs) = go as bs [map (/c) b]
where
val us vs ws = let u = head us in map (/u) $ zipWith (-) vs (head $ multiply [tail us] ws)
go [] _ zs = zs
go _ [] zs = zs
go (x:xs) (y:ys) zs = go xs ys $ (val x y zs):zs
 
triangle::Num a => [[a]] -> [[a]] -> ([[a]],[[a]])
triangle xs bs = triang ([],[]) (xs,bs)
where
triang ts (_,[]) = ts
triang ts ([],_) = ts
triang (os,ps) zs = triang (us:os,cs:ps).unzip $ [(fun tus vs, fun cs es) | (v:vs,es) <- zip uss css,let fun = zipWith (\x y -> v*x - u*y)]
where ((us@(u:tus)):uss,cs:css) = bubble zs
 
bubble:: (Num a, Ord a) => ([[a]],[[a]]) -> ([[a]],[[a]])
bubble (xs,bs) = (go xs, go bs)
where
idmax = snd.maximum.flip zip [0..].map head $ xs
go ys = let (us,vs) = splitAt idmax ys in vs ++ us
main = do
let a = [[1.00, 0.00, 0.00, 0.00, 0.00, 0.00],
[1.00, 0.63, 0.39, 0.25, 0.16, 0.10],
[1.00, 1.26, 1.58, 1.98, 2.49, 3.13],
[1.00, 1.88, 3.55, 6.70, 12.62, 23.80],
[1.00, 2.51, 6.32, 15.88, 39.90, 100.28],
[1.00, 3.14, 9.87, 31.01, 97.41, 306.02]]
let b = [[-0.01], [0.61], [0.91], [0.99], [0.60], [0.02]]
mapM_ print $ gauss a b
</lang>
{{out}}
<pre>
[-1.0e-2]
[1.6027903945021098]
[-1.6132030599055482]
[1.245494121371424]
[-0.49098971958465265]
[6.576069617523134e-2]
</pre>
 
===Another example===
We use Rational numbers for having more precision. a % b is the rational a / b.
<lang Haskell>
678

edits