backward_subst.hs

-- backwards substitution solving

type LEq = [Double]
type Sol = Double
type XVec = [Double]

bsolve' :: [LEq] -> [Sol] -> XVec -> XVec
bsolve' _ [] _ = []
bsolve' (x:xs) (y:ys) v = let 
                              x' = dropWhile (\x -> x == 0) x
                              v' = [if length x' == 1 then y / (head x') else (y - (sum [i*j | (i,j) <- zip (tail x') (tail v)])) / (head x')]
                          in (bsolve' xs ys ([head v] ++ v' ++ tail v)) ++ v'
-- e.g. given 
-- [1,0.5,-2,4].x = -8
-- [0,3,3,2].x = 3
-- [0,0,1,5].x = -4
-- [0,0,0,2].x = 6
-- and using form x = [x1,x2..xn]
-- [[0,0,0,2], [0,0,1,5], [0,3,3,2], [1,0.5,-2,4]] -- eqs: xn, xn-1+xn, ..
-- [6, 3, -4, -8] -- solutions to xn, xn-1+xn, ..
-- start passing [0] for XVec
-- 
-- ghci> bsolve [[0,0,0,2], [0,0,1,5], [0,3,3,2], [1,0.5,-2,4]] [6, -4, 3, -8] [0]
-- [-67.0,18.0,-19.0,3.0] -- [x1, x2, x3 , x4]
-- ghci> bsolve [[0,0,5], [0,1,2], [2,3,-4]] [15,3,10] [0]
-- [15.5,-3.0,3.0]