ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

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

Clojure ヤコビ法 プログラミング 数値計算 連立一次方程式

例として
9x+2y+z+u=20
2x+8y-2z+u=16
-x-2y+7z-2u=8
x-y-2z+6u=17
を考える.
この方程式を上から順に対角線上の変数について解くと
x=(20-2y-z-u)/9
y=(16-2x+2z-u)/8
z=(8+x+2y+2u)/7
u=(17-x+y+2z)/6
となる.
x,y,z,uに適当な値を入れて右辺を計算し、
得られた値を新たなx,y,z,uとして、計算を繰り返す.
漸化式で書くと
x_{n+1}=(20-2y_n-z_n-u_n)/9
y_{n+1}=(16-2x_n+2z_n-u_n)/8
z_{n+1}=(8+x_n+2y_n+2u_n)/7
u_{n+1}=(17-x_n+y_n+2z_n)/6

(def N 4)

(def a [[9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0] [-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0]])
(def b  [20.0 16.0 8.0 17.0])
(def x0 [0.0 0.0 0.0 0.0])

(defn disp_vector [row]
    (doseq [col row]
        (print (format "%14.10f\t" col)))
    (println))

(defn disp_matrix [matrix]
    (doseq [row matrix]
        (doseq [col row]
            (print (format "%14.10f\t" col)))
        (println)))

(defn row_loop [row a b x0 x1]
    (def a1 (map (fn [a b] [a b]) (nth a row) x0))
    (def a2 (concat (take row a1) (drop (+ row 1) a1)))
    (def a3 (map (fn [x] (* (first x) (second x))) a2))
    (def s  (apply + a3))
    (def x  (/ (- (nth b row) s) (nth (nth a row) row)))
    (def xs (cons x x1))

    (if (>= row (- N 1)) 
        (reverse xs)
        (row_loop (inc row) a b x0 xs)))

(defn jacobi [a b x0 xs]
    (def x1 (row_loop 0 a b x0 []))

    (def cnt (count 
        (filter (fn [x] (>= x 0.0000000001))
            (map (fn [x] (. Math abs (- (first x) (second x)))) 
                (map (fn [a b] [a b]) x0 x1)))))

    (if (< cnt 1)
        (vector (reverse xs) x1)
        (jacobi a b x1 (cons x1 xs))))

(def xs (jacobi a b x0 []))

(disp_matrix (first xs))
(println "X")
(disp_vector (second xs))

  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
参考文献