You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
61 lines
1.9 KiB
Haskell
61 lines
1.9 KiB
Haskell
module Main where
|
|
|
|
import Numeric.LinearAlgebra as LA
|
|
import Graphics.EasyPlot
|
|
|
|
data TriangType = Up | Low | None deriving Show
|
|
|
|
a n = build (n, n) $ gen n :: Matrix Double
|
|
|
|
gen n x y | x == (fromIntegral $ n - 1) = 1
|
|
| x == y = 1
|
|
| x < y = -1
|
|
| x > y = 0
|
|
|
|
luGet a = case lu a of
|
|
(l, u , p, s) -> (p, ((l, Low), (u, Up)))
|
|
|
|
solveTriang (a, t) f = case t of
|
|
Up -> let rs = toRows a in
|
|
calcu (cols a) (reverse rs) f
|
|
Low -> let rs = toRows a in
|
|
calcl (cols a) rs f
|
|
None -> error "Not triang"
|
|
|
|
calcu nrs rs f = calc' rs 0 [] where
|
|
calc' [] n us = us
|
|
calc' (r:rs) n us = calc' rs (n + 1) (ui : us) where
|
|
ui = (fi - (subVector 0 n $ vector us) <.> (subVector (nrs - n) n r)) / ri where
|
|
fi = f `atIndex` (nrs - n - 1)
|
|
ri = (r `atIndex` (nrs - n - 1))
|
|
|
|
calcl nrs rs f = calc' rs 0 [] where
|
|
calc' [] n us = reverse us
|
|
calc' (r:rs) n us = calc' rs (n + 1) (ui : us) where
|
|
ui = (fi - ((subVector 0 n $ vector $ reverse us) <.> (subVector 0 n r))) / ri where
|
|
fi = f `atIndex` n
|
|
ri = r `atIndex` n
|
|
|
|
-- Для вычисления числа обусловенности
|
|
norm1 = maximum . map (sum . map abs) . toLists
|
|
norm2 = norm1 . tr
|
|
norm3 m = sqrt $ maximum $ map magnitude $ toList $ eigenvalues $ m LA.<> (tr m)
|
|
|
|
mu norm m = norm m * (norm $ inv m)
|
|
|
|
collectDots norm 1 = [(0, mu norm (a 1))]
|
|
collectDots norm n = (fromIntegral n, mu norm (a n)) : collectDots norm (n - 1)
|
|
|
|
main :: IO ()
|
|
main = do
|
|
-- Пример решения системы n = 5
|
|
let m = a 5
|
|
let f = vector $ take 5 $ repeat 1
|
|
let v = vector $ solveTriang (fst $ snd $ luGet m) f
|
|
let u = solveTriang (snd $ snd $ luGet m) v
|
|
--print u :: Vector R
|
|
|
|
let dots = collectDots norm1 100
|
|
plot X11 $ Data2D [Title "norm1", Style Lines] [] dots
|
|
print dots
|