Haskell で連立一次方程式を解く(ヤコビの反復法)

例として




を考える.
この方程式を上から順に対角線上の変数について解くと




となる.
に適当な値を入れて右辺を計算し、
得られた値を新たなとして、計算を繰り返す.
漸化式で書くと

import Text.Printf
import Debug.Trace
import Control.Monad

n = 4::Int

disp_matrix::[[Double]]->IO()
disp_matrix matrix = do
    forM_ matrix $ \row -> do
        forM_ row $ \elem -> do
            printf "%14.10f\t" elem
        putStrLn ""

disp_vector::[Double]->IO()
disp_vector vector = do
    forM_ vector $ \elem -> do
        printf "%14.10f\t" elem 
    putStrLn ""

row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double]
row_loop row a b x0 x1 =
    let
        a1 = zip (a!!row) x0
        a2 = take row a1 ++ drop (row + 1) a1
        a3 = map (\(x, y) -> x * y) $ a2
        s  = sum a3
        x = ((b!!row) - s) / a!!row!!row
    in
        if row <= 0 
            then
                x:x1
            else
                (row_loop (row-1) a b x0 (x:x1))

jacobi::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double])
jacobi a b x0 xs =
    let
        x1  = (row_loop (n - 1) a b x0 [])
        cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1
    in
        if cnt < 1
            then 
                (reverse xs, x1)
            else 
                (jacobi a b x1 (x1:xs))

main = do
    let a  = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]]
    let b  = [20,16,8,17::Double]
    let x0 = [0,0,0,0::Double]

    putStrLn "A"
    disp_matrix a
    putStrLn "B"
    disp_vector b
    putStrLn ""

    let (xs, x1) = jacobi a b x0 []
    disp_matrix xs

    putStrLn "X"
    disp_vector x1

A
  9.0000000000    2.0000000000    1.0000000000    1.0000000000
  2.0000000000    8.0000000000   -2.0000000000    1.0000000000
 -1.0000000000   -2.0000000000    7.0000000000   -2.0000000000
  1.0000000000   -1.0000000000   -2.0000000000    6.0000000000
B
 20.0000000000   16.0000000000    8.0000000000   17.0000000000

  2.2222222222    2.0000000000    1.1428571429    2.8333333333
  1.3359788360    1.3759920635    2.8412698413    3.1772486772
  1.2477219283    1.9791666667    2.6346371882    3.7870921517
  1.0688819252    1.8733422960    2.9686056521    3.8334531858
  1.0501396189    1.9957492835    2.9260675555    3.9569452792
  1.0139431776    1.9743638243    2.9936469635    3.9662907960
  1.0101482880    1.9991395970    2.9850360597    3.9912857623
  1.0028221093    1.9948112226    2.9987141438    3.9931772381
  1.0020540192    1.9998258539    2.9969712901    3.9982362335
  1.0005711965    1.9989497885    2.9997397420    3.9986190691
  1.0004157346    1.9999647527    2.9993869874    3.9996430127
  1.0001156105    1.9997874366    2.9999473236    3.9997204988
  1.0000841449    1.9999928659    2.9998759259    3.9999277456
  1.0000233996    1.9999569771    2.9999893383    3.9999434288
  1.0000170310    1.9999985561    2.9999748873    3.9999853757
  1.0000047361    1.9999912921    2.9999978421    3.9999885500
  1.0000034471    1.9999997077    2.9999949172    3.9999970400
  1.0000009586    1.9999982375    2.9999995632    3.9999976825
  1.0000006977    1.9999999408    2.9999989712    3.9999994009
  1.0000001940    1.9999996433    2.9999999116    3.9999995309
  1.0000001412    1.9999999880    2.9999997918    3.9999998787
  1.0000000393    1.9999999278    2.9999999821    3.9999999051
  1.0000000286    1.9999999976    2.9999999579    3.9999999755
  1.0000000079    1.9999999854    2.9999999964    3.9999999808
  1.0000000058    1.9999999995    2.9999999915    3.9999999950
  1.0000000016    1.9999999970    2.9999999993    3.9999999961
  1.0000000012    1.9999999999    2.9999999983    3.9999999990
  1.0000000003    1.9999999994    2.9999999999    3.9999999992
  1.0000000002    2.0000000000    2.9999999997    3.9999999998
  1.0000000001    1.9999999999    3.0000000000    3.9999999998
  1.0000000000    2.0000000000    2.9999999999    4.0000000000
X
  1.0000000000    2.0000000000    3.0000000000    4.0000000000

参考文献