(printf "%d\n" (+ 3 5)) (printf "%d\n" (- 3 5)) (printf "%d\n" (* 3 5)) (printf "%f\n" (Math/pow 3 5)) (printf "%f\n" (double (/ 5 3))) (printf "%f\n" (/ 5.0 3)) (printf "%f\n" (/ 5 3.0)) (printf "%d\n" (quot 5 3)) (printf "%d\n" (rem 5 3)) (print (format "%d\n" (* 3 5))) (println (format "%3d" (* 3 5))) (println (format "%23.20f" (/ 5.0 3)))
(def i (* 3 5)) (printf "3 * 5 = %d\n" i)
(loop [n 1] (if (< n 10) (do (printf "%d, " n) (recur (inc n))) (println "")))
(loop [n 1] (if (< n 10) (if (= (mod n 3) 0) (do (printf "%d, " n) (recur (inc n))) (recur (inc n))) (println "")))
(defn total [n] (if (zero? n) 0 (if (= (mod n 3) 0) (+' n (total (dec n))) (total (dec n))))) (println (total 99))
(range 1 10)
(filter #(= (mod % 3) 0) (range 1 10)) (filter #(= (rem % 3) 0) (range 1 10))
(apply + (filter #(= (mod % 3) 0) (range 1 100))) (reduce + (filter #(= (mod % 3) 0) (range 1 100)))
; 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 (defn sn [a, lim] (def n (quot lim a)) ; 項数:n = 上限:lim / 公差:a (def l (* a n)) ; 末項:l = 項数:n * 公差:a (quot (* (+ a l) n) 2)) ; 総和:sn = (初項:a + 末項:l) * 項数:n / 2 ; 3 の倍数の合計 (sn 3 999)
(def n 10000) (quot (* n (+ n 1)) 2)
(def n 5000)
(* n (+ n 1))
(def n 5000) (int (Math/pow n 2))
(def n 1000) (quot (* (* n (+ n 1)) (+ (* 2 n) 1)) 6)
(def n 100) (int (quot (* (Math/pow n 2) (Math/pow (+ n 1) 2)) 4))
(def n 10) (def a 2) (def r 3) (int (quot (* a (- (Math/pow r n) 1)) (- r 1)))
(range 0 10) (map #(+ 5 (* 3 %)) (range 0 10)) (apply * (map #(+ 5 (* 3 %)) (range 0 10)))
; 等差数列の積 (defn product [m d n] (if (zero? n) 1 (* m (product (+ m d) d (- n 1))))) ; 初項 5, 公差 3, 項数 10 の数列の積 (product 5 3 10)
; 階乗を求める関数 (defn Fact [n] (if (zero? n) 1 (* n (Fact (- n 1))))) ; 10の階乗 (Fact 10)
; 下降階乗冪 (defn FallingFact [x n] (if (<= n 1) x (* x (FallingFact (- x 1) (- n 1))))) ; 10 から 6 までの 総乗 (FallingFact 10 5)
; 上昇階乗冪 (defn RisingFact [x n] (if (<= n 1) x (* x (RisingFact (+ x 1) (- n 1))))) ; 10 から 14 までの 総乗 (RisingFact 10 5)
; 階乗 (defn Fact [n] (if (zero? n) 1 (* n (Fact (- n 1))))) ; 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) (def n 10) (def r 5) (quot (Fact n) (Fact (- n r))) ; 下降階乗冪 (defn FallingFact [x n] (if (<= n 1) x (* x (FallingFact (- x 1) (- n 1))))) ; 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) (def n 10) (def r 5) (FallingFact n r)
; 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) (def n 10) (def r 5) (Math/pow n r)
; 組合せ (defn Comb [n r] (cond (or (= r 0) (= r n)) 1 (= r 1) n true (+ (Comb (- n 1) (- r 1)) (Comb (- n 1) r)))) ; 異なる 10 個のものから 5 個取ってできる組合せの総数 (def n 10) (def r 5) (Comb n r)
; 組合せ (defn Comb [n r] (cond (or (= r 0) (= r n)) 1 (= r 1) n true (+ (Comb (- n 1) (- r 1)) (Comb (- n 1) r)))) ; 異なる 10 個のものから重複を許して 5 個取る組合せの総数 (def n 10) (def r 5) (Comb (+ n (- r 1)) r)
;自作の正弦関数 (defn mySin [x n nega numerator denominator y] (def m (* 2 n)) (def denom (* denominator (* (+ m 1) m))) (def nume (* numerator (* x x))) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= a 0.00000000001) y (+ y (mySin x (+ n 1) (not nega) nume denom (if nega a (- a)))))) (doseq [degree (filter #(or (= (mod % 30) 0) (= (mod % 45) 0)) (map #(* % 15) (range 0 25)))] (do (def radian (/ (* degree (. Math PI)) 180.0)) ;自作の正弦関数 (def d1 (mySin radian 1 false radian 1.0 radian)) ;標準の正弦関数 (def d2 (. Math sin radian)) ;標準関数との差異 (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の余弦関数 (defn myCos [x n nega numerator denominator y] (def m (* 2 n)) (def denom (* denominator (* m (- m 1)))) (def nume (* numerator (* x x))) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= a 0.00000000001) y (+ y (myCos x (+ n 1) (not nega) nume denom (if nega a (- a)))))) (doseq [degree (filter #(or (= (mod % 30) 0) (= (mod % 45) 0)) (map #(* % 15) (range 0 25)))] (do (def radian (/ (* degree (. Math PI)) 180.0)) ;自作の余弦関数 (def d1 (myCos radian 1 false 1.0 1.0 1.0)) ;標準の余弦関数 (def d2 (. Math cos radian)) ;標準関数との差異 (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の正接関数 (defn myTan [x x2 n t] (def denom (/ x2 (- n t))) (def nume (- n 2)) (if (<= nume 1) (/ x (- 1 denom)) (myTan x x2 nume denom))) (doseq [degree (map #(- % 90) (filter #(not (zero? (mod % 180))) (map #(* % 15) (range 0 12))))] (do (def radian (/ (* degree (. Math PI)) 180.0)) (def x2 (* radian radian)) ;自作の正接関数 (def d1 (myTan radian x2 15 0)) ;標準の正接関数 (def d2 (. Math tan radian)) ;標準関数との差異 (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2)))))
;自作の指数関数 (defn myExp [x n numerator denominator y] (def denom (* denominator n)) (def nume (* numerator x)) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= (. Math abs a) 0.00000000001) y (+ y (myExp x (+ n 1) nume denom a)))) (doseq [x (map #(/ (- % 10) 4.0) (range 0 21))] (do ;標準の指数関数 (def d1 (. Math exp x)) ;自作の指数関数 (def d2 (myExp x 1 1.0 1.0 1.0)) ;標準関数との差異 (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の指数関数 (defn myExp [x x2 n t] (def denom (/ x2 (+ n t))) (def nume (- n 4)) (if (< nume 6) (+ 1.0 (/ (* 2.0 x) (+ (- 2.0 x) denom))) (myExp x x2 nume denom))) (doseq [x (map #(/ (- % 10) 4.0) (range 0 20))] (do ;標準の指数関数 (def d1 (. Math exp x)) ;自作の指数関数 (def x2 (* x x)) ;標準関数との差異 (def d2 (myExp x x2 30 0.0)) ;30:必要な精度が得られるよう, 6から始めて4ずつ増加させる (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の対数関数 (defn myLog [x numerator denominator y] (def denom (+ denominator 2.0)) (def nume (* numerator (* x2 x2))) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= (. Math abs a) 0.00000000001) y (+ y (myLog x nume denom a)))) (doseq [x (map #(/ % 5.0) (range 1 21))] (do ;標準の対数関数 (def d1 (. Math log x)) ;自作の対数関数 (def x2 (/ (- x 1.0) (+ x 1.0))) ;標準関数との差異 (def d2 (* 2.0 (myLog x2 x2 1.0 x2))) (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の対数関数 (defn myLog [x n t] (do (def n2 (if (zero? (rem n 2)) 2 n)) (def n3 (if (> n 3) n2 n)) (def x2 (if (> n 3) (* x (quot n 2)) x)) (def t2 (/ x2 (+ n3 t))) (if (<= n 2) (/ x (+ 1.0 t2)) (myLog x (- n 1) t2)))) (doseq [x (map #(/ % 5.0) (range 1 21))] (do ;標準の対数関数 (def d1 (. Math log x)) ;自作の対数関数 (def d2 (myLog (- x 1.0) 27 0.0)) ;27:必要な精度が得られる十分大きな奇数 ;標準関数との差異 (println (format "%5.2f : %13.10f - %13.10f = %13.10f" x d1 d2 (- d1 d2)))))
;自作の双曲線正弦関数 (defn mySinh [x n numerator denominator y] (def m (* 2 n)) (def denom (* denominator (* (+ m 1) m))) (def nume (* numerator (* x x))) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= (. Math abs a) 0.00000000001) y (+ y (mySinh x (+ n 1) nume denom a)))) (doseq [i (map #(- % 10) (range 0 21))] (do (def x (double i)) ;自作の双曲線正弦関数 (def d1 (mySinh x 1 x 1.0 x)) ;標準の双曲線正弦関数 (def d2 (. Math sinh x)) ;標準関数との差異 (println (format "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (- d1 d2)))))
;自作の双曲線余弦関数 (defn myCosh [x n numerator denominator y] (def m (* 2 n)) (def denom (* denominator (* m (- m 1)))) (def nume (* numerator (* x x))) (def a (/ nume denom)) ;十分な精度になったら処理を抜ける (if (<= (. Math abs a) 0.00000000001) y (+ y (myCosh x (+ n 1) nume denom a)))) (doseq [i (map #(- % 10) (range 0 21))] (do (def x (double i)) ;自作の双曲線余弦関数 (def d1 (myCosh x 1 1.0 1.0 1.0)) ;標準の双曲線余弦関数 (def d2 (. Math cosh x)) ;標準関数との差異 (println (format "%3d : %17.10f - %17.10f = %13.10f" i d1 d2 (- d1 d2)))))
;自作の逆正接関数 (defn myAtan [x x2 n t] (def m (quot n 2)) (def denom (/ (* (* m m) x2) (+ n t))) (def nume (- n 2)) (if (<= nume 1) (/ x (+ 1 denom)) (myAtan x x2 nume denom))) (doseq [degree (map #(- (* % 15) 45) (range 0 7))] (do (def radian (/ (* degree (. Math PI)) 180.0)) (def x2 (* radian radian)) ;自作の逆正接関数 (def d1 (myAtan radian x2 23 0)) ;標準の逆正接関数 (def d2 (. Math atan radian)) ;標準関数との差異 (println (format "%3d : %13.10f - %13.10f = %13.10f" degree d1 d2 (- d1 d2))))) -45 : -0.6657737500 - -0.6657737500 = 0.0000000000 -30 : -0.4823479071 - -0.4823479071 = 0.0000000000 -15 : -0.2560527700 - -0.2560527700 = 0.0000000000 0 : 0.0000000000 - 0.0000000000 = 0.0000000000 15 : 0.2560527700 - 0.2560527700 = 0.0000000000 30 : 0.4823479071 - 0.4823479071 = -0.0000000000 45 : 0.6657737500 - 0.6657737500 = -0.0000000000 nil (doseq [i (map #(+ (* % 2) 1) (range 5 16))] (do (def radian 1) (def x2 (* radian radian)) ;自作の逆正接関数 (def d1 (myAtan radian x2 i 0)) ;標準関数との差異 (def pi (* d1 4.0)) (println (format "%2d : %13.10f, %13.10f" i pi (- (pi (. Math PI)))))) 11 : 3.1414634146, -0.0001292390 13 : 3.1416149068, 0.0000222532 15 : 3.1415888251, -0.0000038285 17 : 3.1415933119, 0.0000006583 19 : 3.1415925404, -0.0000001131 21 : 3.1415926730, 0.0000000194 23 : 3.1415926503, -0.0000000033 25 : 3.1415926542, 0.0000000006 27 : 3.1415926535, -0.0000000001 29 : 3.1415926536, 0.0000000000 31 : 3.1415926536, -0.0000000000 nil
(defn f[x] (/ 4 (+ 1 (* x x)))) (def a 0) (def b 1) ; 台形則で積分 (doseq [j (range 1 11)] (def n (Math/pow 2 j)) (def h (/ (- b a) n)) (def x a) (def s 0) (doseq [i (range 1 n)] (def x (+ x h)) (def s (+ s (f x)))) (def t1 (double (* h (+ (/ (+ (f a) (f b)) 2) s)))) (def t2 (- t1 (. Math PI))) ; 結果を π と比較 (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x] (/ 4 (+ 1 (* x x)))) (def a 0) (def b 1) ; 中点則で積分 (doseq [j (range 1 11)] (def n (Math/pow 2 j)) (def h (/ (- b a) n)) (def x (+ a (/ h 2))) (def s 0) (doseq [i (range 1 (+ n 1))] (def s (+ s (f x))) (def x (+ x h))) (def t1 (double (* h s))) (def t2 (- t1 (. Math PI))) ; 結果を π と比較 (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x] (/ 4 (+ 1 (* x x)))) (def a 0) (def b 1) ; Simpson則で積分 (doseq [j (range 1 6)] (def n (Math/pow 2 j)) (def h (/ (- b a) n)) (def x (+ a h)) (def s2 0) (def s4 0) (doseq [i (range 1 (+ (/ n 2) 1))] (def s4 (+ s4 (f x))) (def x (+ x h)) (def s2 (+ s2 (f x))) (def x (+ x h)) ) (def w2 (double (+ (* (- s2 (f b)) 2) (f a) (f b)))) (def w4 (* s4 4)) (def t1 (double (/ (* (+ w2 w4) h) 3))) (def t2 (- t1 (. Math PI))) ; 結果を π と比較 (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
(defn f[x] (/ 4 (+ 1 (* x x)))) (def a 0) (def b 1) (def t (list)) ; 台形則で積分 (doseq [j (range 1 11)] (def n (Math/pow 2 j)) (def h (/ (- b a) n)) (def x a) (def s 0) (doseq [i (range 1 n)] (def x (+ x h)) (def s (+ s (f x)))) (def t1 (double (* h (+ (/ (+ (f a) (f b)) 2) s)))) (def t2 (- t1 (. Math PI))) ; 結果を保存 (def t (cons t1 t))) ; Richardsonの補外法 (def t (reverse t)) (doseq [j (range 2 7)] (def n (Math/pow 4 (- j 1))) (def s (list)) (doseq [i (range j 7)] (def t1 (+ (second t) (/ (- (second t) (first t)) (- n 1)))) (def t2 (- t1 (. Math PI))) ; 結果を π と比較 (if (= i j) (println (format "%2d : %13.10f, %13.10f" j t1 t2))) (def t (rest t)) (def s (cons t1 s))) (def t (reverse s)))
(def N 7) ; 元の関数 (defn f[x] (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2)))))) ; Lagrange (ラグランジュ) 補間 (defn lagrange [d x y] (def sum_list (list 0.0)) (doseq [i (range 0 N)] (def prod_list (list (nth y i))) (doseq [j (range 0 N)] (if (= i j) (def prod_list (cons 1 prod_list)) (do (def w1 (/ (- d (nth x j)) (- (nth x i) (nth x j)))) (def prod_list (cons w1 prod_list))))) (def w2 (reduce * prod_list)) (def sum_list (cons w2 sum_list))) (reduce + sum_list)) ; 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (def x (map #(- (* % 1.5) 4.5) (range 0 N))) (def y (map #(f %) x)) ; 0.5刻みで 与えられていない値を補間 (def d1 (map #(- (* % 0.5) 4.5) (range 0 19))) (def d2 (map #(f %) d1)) (def d3 (map #(lagrange % x y) d1)) (doseq [d (map list d1 d2 d3)] (def d1 (nth d 0)) (def d2 (nth d 1)) (def d3 (nth d 2)) (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7) ; 元の関数 (defn f[x] (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2)))))) ; Neville (ネヴィル) 補間 (defn neville [d x y] (def w (cons y nil)) (doseq [j (range 1 N)] (def t_list (list)) (doseq [i (range 0 (- N j))] (def t (nth w 0)) (def t_list (cons (+ (nth t (+ i 1)) (/ (* (- (nth t (+ i 1)) (nth t i)) (- d (nth x (+ i j)))) (- (nth x (+ i j)) (nth x i)))) t_list ) ) ) (def w (cons (reverse t_list) w)) ) (nth (nth w 0) 0) ) ; 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (def x (map #(- (* % 1.5) 4.5) (range 0 N))) (def y (map #(f %) x)) ; 0.5刻みで 与えられていない値を補間 (def d1 (map #(- (* % 0.5) 4.5) (range 0 19))) (def d2 (map #(f %) d1)) (def d3 (map #(neville % x y) d1)) (doseq [d (map list d1 d2 d3)] (def d1 (nth d 0)) (def d2 (nth d 1)) (def d3 (nth d 2)) (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7) ; 元の関数 (defn f[x] (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2)))))) ; Newton (ニュートン) 補間 (defn newton [d x a] (def sum_list (list (nth a 0))) (doseq [i (range 1 N)] (def prod_list (list (nth a i))) (doseq [j (range 0 i)] (def prod_list (cons (- d (nth x j)) prod_list))) (def w (reduce * prod_list)) (def sum_list (cons w sum_list))) (reduce + sum_list)) ; 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (def x (map #(- (* % 1.5) 4.5) (range 0 N))) (def y (map #(f %) x)) ; 差分商の表を作る (def d (cons y nil)) (doseq [i (range 1 N)] (def w (nth d 0)) (def t (list)) (doseq [j (range 0 (- N i))] (def t (cons (/ (- (nth w (+ j 1)) (nth w j)) (- (nth x (+ j i)) (nth x j))) t ) ) ) (def d (cons (reverse t) d)) ) (def d (reverse d)) ; n階差分商 (def a (map #(nth (nth d %) 0) (range 0 N))) ; 0.5刻みで 与えられていない値を補間 (def d1 (map #(- (* % 0.5) 4.5) (range 0 19))) (def d2 (map #(f %) d1)) (def d3 (map #(newton % x a) d1)) (doseq [d (map list d1 d2 d3)] (def d1 (nth d 0)) (def d2 (nth d 1)) (def d3 (nth d 2)) (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7) (def Nx2 14) ; 元の関数 (defn f[x] (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2)))))) ; 導関数 (defn fd[x] (+ (- 1.0 (/ (Math/pow x 2.0) 2)) (/ (Math/pow x 4.0) (* 4 (* 3 2))))) ; Hermite (エルミート) 補間 (defn hermite [d z a] (def sum_list (list (nth a 0))) (doseq [i (range 1 Nx2)] (def prod_list (list (nth a i))) (doseq [j (range 0 i)] (def prod_list (cons (- d (nth z j)) prod_list))) (def w (reduce * prod_list)) (def sum_list (cons w sum_list))) (reduce + sum_list)) ; 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (def x (map #(- (* % 1.5) 4.5) (range 0 N))) (def y (map #(f %) x)) (def yd (map #(fd %) x)) ; 差分商の表を作る (def z (map #(nth x %) (map #(quot % 2) (range 0 Nx2)))) (def w (map #(nth y %) (map #(quot % 2) (range 0 Nx2)))) (def d (cons w nil)) (doseq [i (range 1 Nx2)] (def w (nth d 0)) (def t (list)) (doseq [j (range 0 (- Nx2 i))] (def t (cons (if (and (= i 1) (= (rem j 2) 0)) (nth yd (quot j 2)) (/ (- (nth w (+ j 1)) (nth w j)) (- (nth z (+ j i)) (nth z j))) ) t ) ) ) (def d (cons (reverse t) d)) ) (def d (reverse d)) ; n階差分商 (def a (map #(nth (nth d %) 0) (range 0 Nx2))) ; 0.5刻みで 与えられていない値を補間 (def d1 (map #(- (* % 0.5) 4.5) (range 0 19))) (def d2 (map #(f %) d1)) (def d3 (map #(hermite % z a) d1)) (doseq [d (map list d1 d2 d3)] (def d1 (nth d 0)) (def d2 (nth d 1)) (def d3 (nth d 2)) (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
(def N 7) ; 元の関数 (defn f[x] (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2)))))) ; Spline (スプライン) 補間 (defn spline [d x y z] ; 補間関数値がどの区間にあるか (def k -1) (doseq [i (range 1 N)] (if (and (<= d (nth x i)) (< k 0)) (def k (- i 1)) ) ) (if (< k 0) (def k N)) (def d1 (- (nth x (+ k 1)) d)) (def d2 (- d (nth x k))) (def d3 (- (nth x (+ k 1)) (nth x k))) (def a1 (* (nth z k) (Math/pow d1 3.0))) (def a2 (* (nth z (+ k 1)) (Math/pow d2 3.0))) (def a3 (/ (+ a1 a2) (* 6.0 d3))) (def b1 (/ (nth y k) d3)) (def b2 (* (nth z k) d3)) (def b3 (/ b2 6.0)) (def b4 (* (- b1 b3) d1)) (def c1 (/ (nth y (+ k 1)) d3)) (def c2 (* (nth z (+ k 1)) d3)) (def c3 (/ c2 6.0)) (def c4 (* (- c1 c3) d2)) (+ (+ a3 b4) c4) ) ; 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット (def x (map #(- (* % 1.5) 4.5) (range 0 N))) (def y (map #(f %) x)) ; 3項方程式の係数の表を作る (def a (list 0.0)) (def b (list 0.0)) (def c (list 0.0)) (def d (list 0.0)) (doseq [i (range 1 (- N 1))] (def a (cons (- (nth x i ) (nth x (- i 1))) a)) (def b (cons (* 2.0 (- (nth x (+ i 1)) (nth x (- i 1)))) b)) (def c (cons (- (nth x (+ i 1)) (nth x i )) c)) (def d (cons (* 6.0 (- (/ (- (nth y(+ i 1)) (nth y i)) (- (nth x(+ i 1)) (nth x i))) (/ (- (nth y i) (nth y (- i 1))) (- (nth x i) (nth x (- i 1)))))) d))) (def a (reverse a)) (def b (reverse b)) (def c (reverse c)) (def d (reverse d)) ; 3項方程式を解く (ト−マス法) (def g (list (nth b 1) 0.0)) (def s (list (nth d 1) 0.0)) (doseq [i (range 2 (- N 1))] (def w (nth g 0)) (def g (cons (- (nth b i) (/ (* (nth a i) (nth c (- i 1))) w)) g)) (def s (cons (- (nth d i) (/ (* (nth a i) (nth s 0)) w)) s)) ) (def g (reverse g)) (def s (reverse s)) (def z (list 0.0)) (def z (cons (/ (nth s (- N 2)) (nth g (- N 2))) z)) (doseq [i (reverse (range 1 (- N 2)))] (def z (cons (/ (- (nth s i) (* (nth c i) (nth z 0))) (nth g i)) z)) ) (def z (cons 0.0 z)) ; 0.5刻みで 与えられていない値を補間 (def d1 (map #(- (* % 0.5) 4.5) (range 0 19))) (def d2 (map #(f %) d1)) (def d3 (map #(spline % x y z) d1)) (doseq [d (map list d1 d2 d3)] (def d1 (nth d 0)) (def d2 (nth d 1)) (def d3 (nth d 2)) (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ; Euler法 (defn euler[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置 (def wx (+ x (* h vx))) (def wy (+ y (* h vy))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f" t vx vy wx wy)) ; 速度 (def wvx (+ vx (* h (fx vx vy)))) (def wvy (+ vy (* h (fy vx vy)))) (if (>= wy 0.0) (euler (+ i 1) wvx wvy wx wy) nil)) (euler 1 vx vy x y)
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ; Heun法 (defn heun[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置・速度 (def wx2 (+ x (* h vx ))) (def wy2 (+ y (* h vy ))) (def wvx2 (+ vx (* h (fx vx vy) ))) (def wvy2 (+ vy (* h (fy vx vy) ))) (def wx (+ x (* h (/ (+ vx wvx2 ) 2)))) (def wy (+ y (* h (/ (+ vy wvy2 ) 2)))) (def wvx (+ vx (* h (/ (+ (fx vx vy) (fx wvx2 wvy2)) 2)))) (def wvy (+ vy (* h (/ (+ (fy vx vy) (fy wvx2 wvy2)) 2)))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy)) (if (>= wy 0.0) (heun (+ i 1) wvx wvy wx wy) nil)) (heun 1 vx vy x y)
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ;中点法 (defn midpoint[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置・速度 (def wvx1 (* h (fx vx vy))) (def wvy1 (* h (fy vx vy))) (def wvx2 (+ vx (/ wvx1 2.0))) (def wvy2 (+ vy (/ wvy1 2.0))) (def wvx (+ vx (* h (fx wvx2 wvy2)))) (def wvy (+ vy (* h (fy wvx2 wvy2)))) (def wx (+ x (* h wvx2))) (def wy (+ y (* h wvy2))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy)) (if (>= wy 0.0) (midpoint (+ i 1) wvx wvy wx wy) nil)) (midpoint 1 vx vy x y)
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ; Runge-Kutta法 (defn rungekutta[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置・速度 (def wx1 (* h vx)) (def wy1 (* h vy)) (def wvx1 (* h (fx vx vy))) (def wvy1 (* h (fy vx vy))) (def wvx5 (+ vx (/ wvx1 2.0))) (def wvy5 (+ vy (/ wvy1 2.0))) (def wx2 (* h wvx5)) (def wy2 (* h wvy5)) (def wvx2 (* h (fx wvx5 wvy5))) (def wvy2 (* h (fy wvx5 wvy5))) (def wvx6 (+ vx (/ wvx2 2.0))) (def wvy6 (+ vy (/ wvy2 2.0))) (def wx3 (* h wvx6)) (def wy3 (* h wvy6)) (def wvx3 (* h (fx wvx6 wvy6))) (def wvy3 (* h (fy wvx6 wvy6))) (def wvx7 (+ vx wvx3)) (def wvy7 (+ vy wvy3)) (def wx4 (* h wvx7)) (def wy4 (* h wvy7)) (def wvx4 (* h (fx wvx7 wvy7))) (def wvy4 (* h (fy wvx7 wvy7))) (def wx (+ x (/ (+ wx1 (* wx2 2.0) (* wx3 2.0) wx4) 6.0))) (def wy (+ y (/ (+ wy1 (* wy2 2.0) (* wy3 2.0) wy4) 6.0))) (def wvx (+ vx (/ (+ wvx1 (* wvx2 2.0) (* wvx3 2.0) wvx4) 6.0))) (def wvy (+ vy (/ (+ wvy1 (* wvy2 2.0) (* wvy3 2.0) wvy4) 6.0))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy)) (if (>= wy 0.0) (rungekutta (+ i 1) wvx wvy wx wy) nil)) (rungekutta 1 vx vy x y)
; 重力加速度 (def g -9.8) ; 空気抵抗係数 (def k -0.01) ; 時間間隔(秒) (def h 0.01) ; 角度 (def degree 45.0) (def radian (* degree (/ (. Math PI) 180.0))) ; 初速 250 km/h -> 秒速に変換 (def v (quot (* 250 1000) 3600)) ; 水平方向の速度 (def vx (* v (. Math cos radian))) ; 鉛直方向の速度 (def vy (* v (. Math sin radian))) ; 位置 (def x 0.0) (def y 0.0) ; 空気抵抗による水平方向の減速分 (defn fx[vx vy] (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vx))) ; 重力と空気抵抗による鉛直方向の減速分 (defn fy[vx vy] (+ g (* k (* (. Math sqrt (+ (* vx vx) (* vy vy))) vy)))) ; Runge-Kutta-Gill法 (defn rungekuttagill[i vx vy x y] ; 経過秒数 (def t (* i h)) ; 位置・速度 (def wx1 (* h vx)) (def wy1 (* h vy)) (def wvx1 (* h (fx vx vy))) (def wvy1 (* h (fy vx vy))) (def wvx5 (+ vx (/ wvx1 2.0))) (def wvy5 (+ vy (/ wvy1 2.0))) (def wx2 (* h wvx5)) (def wy2 (* h wvy5)) (def wvx2 (* h (fx wvx5 wvy5))) (def wvy2 (* h (fy wvx5 wvy5))) (def wvx6 (+ vx (* wvx1 (/ (- (. Math sqrt 2.0) 1.0) 2.0)) (* wvx2 (- 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wvy6 (+ vy (* wvy1 (/ (- (. Math sqrt 2.0) 1.0) 2.0)) (* wvy2 (- 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wx3 (* h wvx6)) (def wy3 (* h wvy6)) (def wvx3 (* h (fx wvx6 wvy6))) (def wvy3 (* h (fy wvx6 wvy6))) (def wvx7 (+ (- vx (/ wvx2 (. Math sqrt 2.0))) (* wvx3 (+ 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wvy7 (+ (- vy (/ wvy2 (. Math sqrt 2.0))) (* wvy3 (+ 1.0 (/ 1.0 (. Math sqrt 2.0)))))) (def wx4 (* h wvx7)) (def wy4 (* h wvy7)) (def wvx4 (* h (fx wvx7 wvy7))) (def wvy4 (* h (fy wvx7 wvy7))) (def wx (+ x (/ (+ wx1 (* wx2 (- 2.0 (. Math sqrt 2.0))) (* wx3 (+ 2.0 (. Math sqrt 2.0))) wx4) 6.0))) (def wy (+ y (/ (+ wy1 (* wy2 (- 2.0 (. Math sqrt 2.0))) (* wy3 (+ 2.0 (. Math sqrt 2.0))) wy4) 6.0))) (def wvx (+ vx (/ (+ wvx1 (* wvx2 (- 2.0 (. Math sqrt 2.0))) (* wvx3 (+ 2.0 (. Math sqrt 2.0))) wvx4) 6.0))) (def wvy (+ vy (/ (+ wvy1 (* wvy2 (- 2.0 (. Math sqrt 2.0))) (* wvy3 (+ 2.0 (. Math sqrt 2.0))) wvy4) 6.0))) (println (format "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%9.5f" t wvx wvy wx wy)) (if (>= wy 0.0) (rungekuttagill (+ i 1) wvx wvy wx wy) nil)) (rungekuttagill 1 vx vy x y)
(defn f[x] (- (* x x) 2.0)) (defn bisection [a b] ; 区間 (a, b) の中点 c = (a + b) / 2 (def c (/ (+ a b) 2.0)) (println (format "%12.10f\t%13.10f" c (- c (Math/sqrt 2)))) (def fc (f c)) (if (< (. Math abs fc) 0.00000000001) c (if (< fc 0) (bisection c b) (bisection a c)))) (def a 1.0) (def b 2.0) (println (format "%12.10f" (bisection a b)))
(defn f[x] (- (* x x) 2.0)) (defn falseposition [a b] ; 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 (def c (/ (- (* a (f b)) (* b (f a))) (- (f b) (f a)))) (println (format "%12.10f\t%13.10f" c (- c (Math/sqrt 2)))) (def fc (f c)) (if (< (. Math abs fc) 0.00000000001) c (if (< fc 0) (falseposition c b) (falseposition a c)))) (def a 1.0) (def b 2.0) (println (format "%12.10f" (falseposition a b)))
(defn g[x] (+ (/ x 2.0) (/ 1.0 x))) (defn iterative [x0] (def x1 (g x0)) (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2)))) (if (< (. Math abs (- x1 x0)) 0.00000000001) x1 (iterative x1))) (def x 1.0) (println (format "%12.10f" (iterative x)))
(defn f0[x] (- (* x x) 2.0)) (defn f1[x] (* 2.0 x)) (defn newton [x0] (def x1 (- x0 (/ (f0 x0) (f1 x0)))) (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2)))) (if (< (. Math abs (- x1 x0)) 0.00000000001) x1 (newton x1))) (def x 2.0) (println (format "%12.10f" (newton x)))
(defn f0[x] (- (* x x) 2.0)) (defn f1[x] (* 2.0 x)) (defn f2[x] 2.0) (defn bailey [x0] (def x1 (- x0 (/ (f0 x0) (- (f1 x0) (/ (* (f0 x0) (f2 x0)) (* 2.0 (f1 x0))))))) (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2)))) (if (< (. Math abs (- x1 x0)) 0.00000000001) x1 (bailey x1))) (def x 2.0) (println (format "%12.10f" (bailey x)))
(defn f[x] (- (* x x) 2.0)) (defn secant [x0 x1] (def x2 (- x1 (/ (* (f x1) (- x1 x0)) (- (f x1) (f x0))))) (println (format "%12.10f\t%13.10f" x2 (- x2 (Math/sqrt 2)))) (if (< (. Math abs (- x2 x1)) 0.00000000001) x2 (secant x1 x2))) (def x0 1.0) (def x1 2.0) (println (format "%12.10f" (secant x0 x1)))
(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]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (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))
(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]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (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)) (def x2 (concat (take (inc row) (reverse xs)) (drop (inc row) x0))) (if (>= row (- N 1)) (reverse xs) (row_loop (inc row) a b x2 xs))) (defn gauss [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) (gauss a b x1 (cons x1 xs)))) ;ガウス・ザイデル法 (def xs (gauss a b x0 [])) (disp_matrix (first xs)) (println "X") (disp_vector (second xs))
(def N 4) (def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]]) (def b [8.0 17.0 20.0 16.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;各列で 一番値が大きい行を 探す (defn get_max_row [row col a max_row max_val] ;一番値が大きい行 (def ws1 (if (< max_val (. Math abs (nth (nth a row) col))) (vector row (. Math abs (nth (nth a row) col))) (vector max_row max_val))) (def max_row2 (first ws1)) (def max_val2 (second ws1)) (if (>= row (dec (count a))) (vector max_row2 max_val2) (get_max_row (inc row) col a max_row2 max_val2))) ;ピボット選択 (defn pivoting [pivot a b a2 b2] (def ws1 (get_max_row 0 pivot a 0 0.0)) (def max_row (first ws1)) (def max_val (second ws1)) (def a3 (cons (nth a max_row) a2)) (def b3 (cons (nth b max_row) b2)) (def a4 (concat (take max_row a) (drop (inc max_row) a))) (def b4 (concat (take max_row b) (drop (inc max_row) b))) (if (>= pivot (dec N)) (vector (reverse a3) (reverse b3)) (pivoting (inc pivot) a4 b4 a3 b3))) ;前進消去 (defn forward_elim_loop [pivot row a b] (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot))) (def a_map (map (fn [a_piv a_row] (- a_row (* a_piv s))) (drop pivot (nth a pivot)) (drop pivot (nth a row)))) (def a1 (concat (take pivot (nth a row)) a_map)) (def a2 (concat (take row a) (cons (vec a1) (drop (inc row) a)))) (def x (- (nth b row) (* (nth b pivot) s))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (< row (dec N)) (forward_elim_loop pivot (inc row) a2 b2) (vector a2 b2))) (defn forward_elimination [pivot a b] (def ws1 (if (< pivot (dec N)) (forward_elim_loop pivot (inc pivot) a b) (vector a b))) (def a2 (first ws1)) (def b2 (second ws1)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2 b2) (vector a2 b2))) ;後退代入 (defn backward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (drop (inc row) (nth a row)) (drop (inc row) b))) (def s (apply + a_map)) (def x (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) ;ピボット選択 (def ws1 (pivoting 0 a b [] [])) (def a1 (first ws1)) (def b1 (vec (second ws1))) (println "pivoting") (println "A") (disp_matrix a1) (println "B") (disp_vector b1) (println "") ;前進消去 (def ws2 (forward_elimination 0 a1 b1)) (def a2 (first ws2)) (def b2 (vec (second ws2))) (println "forward elimination") (println "A") (disp_matrix a2) (println "B") (disp_vector b2) (println "") ;後退代入 (def x (backward_substitution (dec N) a2 b2)) (println "X") (disp_vector (reverse x))
(def N 4) (def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]]) (def b [8.0 17.0 20.0 16.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;各列で 一番値が大きい行を 探す (defn get_max_row [row col a max_row max_val] ;一番値が大きい行 (def ws1 (if (< max_val (. Math abs (nth (nth a row) col))) (vector row (. Math abs (nth (nth a row) col))) (vector max_row max_val))) (def max_row2 (first ws1)) (def max_val2 (second ws1)) (if (>= row (dec (count a))) (vector max_row2 max_val2) (get_max_row (inc row) col a max_row2 max_val2))) ;ピボット選択 (defn pivoting [pivot a b a2 b2] (def ws1 (get_max_row 0 pivot a 0 0.0)) (def max_row (first ws1)) (def max_val (second ws1)) (def a3 (cons (nth a max_row) a2)) (def b3 (cons (nth b max_row) b2)) (def a4 (concat (take max_row a) (drop (inc max_row) a))) (def b4 (concat (take max_row b) (drop (inc max_row) b))) (if (>= pivot (dec N)) (vector (reverse a3) (reverse b3)) (pivoting (inc pivot) a4 b4 a3 b3))) ;前進消去 (defn forward_elim_loop [pivot row a b] (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot))) (def a_map (map (fn [a_piv a_row] (- a_row (* a_piv s))) (nth a pivot) (nth a row))) (def a2 (concat (take row a) (cons (vec a_map) (drop (inc row) a)))) (def x (- (nth b row) (* (nth b pivot) s))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (def next_row (if (= (inc row) pivot) (+ row 2) (inc row))) (if (< next_row N) (forward_elim_loop pivot next_row a2 b2) (vector a2 b2))) (defn forward_elimination [pivot a b] (def row (if (= pivot 0) 1 0)) (def ws1 (if (< pivot N) (forward_elim_loop pivot row a b) (vector a b))) (def a2 (first ws1)) (def b2 (second ws1)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2 b2) (vector a2 b2))) ;後退代入 (defn backward_substitution [row a b] (def x (/ (nth b row) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) ;ピボット選択 (def ws1 (pivoting 0 a b [] [])) (def a1 (first ws1)) (def b1 (vec (second ws1))) (println "pivoting") (println "A") (disp_matrix a1) (println "B") (disp_vector b1) (println "") ;前進消去 (def ws2 (forward_elimination 0 a1 b1)) (def a2 (first ws2)) (def b2 (vec (second ws2))) (println "forward elimination") (println "A") (disp_matrix a2) (println "B") (disp_vector b2) (println "") ;後退代入 (def x (backward_substitution (dec N) a2 b2)) (println "X") (disp_vector (reverse x))
(def N 4) (def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]]) (def b [8.0 17.0 20.0 16.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;各列で 一番値が大きい行を 探す (defn get_max_row [row col a max_row max_val] ;一番値が大きい行 (def ws1 (if (< max_val (. Math abs (nth (nth a row) col))) (vector row (. Math abs (nth (nth a row) col))) (vector max_row max_val))) (def max_row2 (first ws1)) (def max_val2 (second ws1)) (if (>= row (dec (count a))) (vector max_row2 max_val2) (get_max_row (inc row) col a max_row2 max_val2))) ;ピボット選択 (defn pivoting [pivot a b a2 b2] (def ws1 (get_max_row 0 pivot a 0 0.0)) (def max_row (first ws1)) (def max_val (second ws1)) (def a3 (cons (nth a max_row) a2)) (def b3 (cons (nth b max_row) b2)) (def a4 (concat (take max_row a) (drop (inc max_row) a))) (def b4 (concat (take max_row b) (drop (inc max_row) b))) (if (>= pivot (dec N)) (vector (reverse a3) (reverse b3)) (pivoting (inc pivot) a4 b4 a3 b3))) ;前進消去 (defn forward_elim_loop [pivot row a b] (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot))) (def a_map (map (fn [a_piv a_row] (- a_row (* a_piv s))) (drop (inc pivot) (nth a pivot)) (drop (inc pivot) (nth a row)))) (def a1 (concat (take pivot (nth a row)) (cons s a_map))) (def a2 (concat (take row a) (cons (vec a1) (drop (inc row) a)))) (def x (- (nth b row) (* (nth b pivot) s))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (< row (dec N)) (forward_elim_loop pivot (inc row) a2 b) (vector a2 b))) (defn forward_elimination [pivot a b] (def ws1 (if (< pivot (dec N)) (forward_elim_loop pivot (inc pivot) a b) (vector a b))) (def a2 (first ws1)) (def b2 (second ws1)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2 b2) (vector a2 b2))) ;前進代入 (defn forward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (take row (nth a row)) b)) (def s (apply + a_map)) (def y (- (nth b row) s)) (def b2 (concat (take row b) (cons y (drop (inc row) b)))) (if (< row (dec N)) (cons y (forward_substitution (inc row) a b2)) (cons y []))) ;後退代入 (defn backward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (drop (inc row) (nth a row)) (drop (inc row) b))) (def s (apply + a_map)) (def x (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) ;ピボット選択 (def ws1 (pivoting 0 a b [] [])) (def a1 (first ws1)) (def b1 (vec (second ws1))) (println "pivoting") (println "A") (disp_matrix a1) (println "B") (disp_vector b1) (println "") ;前進消去 (def ws2 (forward_elimination 0 a1 b1)) (def a2 (first ws2)) (def b2 (vec (second ws2))) (println "LU") (disp_matrix a2) ;Ly=b から y を求める (前進代入) (def y (forward_substitution 0 a2 b2)) (println "Y") (disp_vector (reverse y)) ;Ux=y から x を求める (後退代入) (def x (backward_substitution (dec N) a2 y)) (println "X") (disp_vector (reverse x))
(def N 4) (def a [[5.0 2.0 3.0 4.0] [2.0 10.0 6.0 7.0] [3.0 6.0 15.0 9.0] [4.0 7.0 9.0 20.0]]) (def b [34.0 68.0 96.0 125.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;前進消去 (defn forward_elim_loop [pivot row a a_sqrt] (def a_map (map (fn [a_row a_piv] (* a_row a_piv)) (take pivot (nth a row)) (take pivot (nth a pivot)))) (def s (apply + a_map)) (def a_col (/ (- (nth (nth a row) pivot) s) a_sqrt)) (def a_row (concat (take pivot (nth a row )) (cons a_col (drop (inc pivot) (nth a row))))) (def a_piv (concat (take row (nth a pivot)) (cons a_col (drop (inc row) (nth a pivot))))) (def a1 (concat (take pivot a ) (cons (vec a_piv) (drop (inc pivot) a)))) (def a2 (concat (take row a1) (cons (vec a_row) (drop (inc row) a)))) (if (< row (dec N)) (forward_elim_loop pivot (inc row) a2 a_sqrt) a2)) (defn forward_elimination [pivot a] (def a_map (map (fn [col] (* col col)) (take pivot (nth a pivot)))) (def s (apply + a_map)) ;ここで根号の中が負の値になると計算できない! (def a_sqrt (. Math sqrt (- (nth (nth a pivot) pivot) s))) (def a2 (forward_elim_loop pivot pivot a a_sqrt)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2) a2)) ;前進代入 (defn forward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (take row (nth a row)) b)) (def s (apply + a_map)) (def y (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons y (drop (inc row) b)))) (if (< row (dec N)) (cons y (forward_substitution (inc row) a b2)) (cons y []))) ;後退代入 (defn backward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (drop (inc row) (nth a row)) (drop (inc row) b))) (def s (apply + a_map)) (def x (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) (println "A") (disp_matrix a) (println "B") (disp_vector b) (println "") ;前進消去 (def a2 (forward_elimination 0 a)) (println "LL^T") (disp_matrix a2) ;Ly=b から y を求める (前進代入) (def y (forward_substitution 0 a2 b)) (println "Y") (disp_vector y) ;L^Tx=y から x を求める (後退代入) (def x (backward_substitution (dec N) a2 y)) (println "X") (disp_vector (reverse x))
(def N 4) (def a [[5.0 2.0 3.0 4.0] [2.0 10.0 6.0 7.0] [3.0 6.0 15.0 9.0] [4.0 7.0 9.0 20.0]]) (def b [34.0 68.0 96.0 125.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;前進消去 (defn forward_elim_loop [pivot col a] (def a_map (map (fn [k a_piv a_col] (* (* a_piv a_col) (nth (nth a k) k))) (range 0 col) (take col (nth a pivot)) (take col (nth a col)))) (def s (if (> (count a_map) 0) (apply + a_map) 0)) (def a_col (/ (- (nth (nth a pivot) col) s) (nth (nth a col) col))) (def a_row (concat (take col (nth a pivot)) (cons a_col (drop (inc col) (nth a pivot))))) (def a_piv (concat (take pivot (nth a col )) (cons a_col (drop (inc pivot) (nth a col))))) (def a1 (concat (take col a ) (cons (vec a_piv) (drop (inc col) a)))) (def a2 (concat (take pivot a1) (cons (vec a_row) (drop (inc pivot) a)))) (def next_col (if (= (inc col) pivot) (inc pivot) (inc col))) (if (< next_col pivot) (forward_elim_loop pivot next_col a2) a2)) (defn forward_elimination [pivot a] (def next_col (if (= pivot 0) 1 0)) (def a1 (forward_elim_loop pivot next_col a)) (def a_map (map (fn [k col] (* (* col col) (nth (nth a1 k) k))) (range 0 pivot) (take pivot (nth a1 pivot)))) (def s (apply + a_map)) (def a_col (- (nth (nth a1 pivot) pivot) s)) (def a_row (concat (take pivot (nth a1 pivot)) (cons a_col (drop (inc pivot) (nth a1 pivot))))) (def a2 (concat (take pivot a1) (cons a_row (drop (inc pivot) a1)))) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2) a2)) ;前進代入 (defn forward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (take row (nth a row)) b)) (def s (apply + a_map)) (def y (- (nth b row) s)) (def b2 (concat (take row b) (cons y (drop (inc row) b)))) (if (< row (dec N)) (cons y (forward_substitution (inc row) a b2)) (cons y []))) ;後退代入 (defn backward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* (* a_col b_col) (nth (nth a row) row))) (drop (inc row) (nth a row)) (drop (inc row) b))) (def s (apply + a_map)) (def x (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) (println "A") (disp_matrix a) (println "B") (disp_vector b) (println "") ;前進消去 (def a2 (forward_elimination 0 a)) (println "LDL^T") (disp_matrix a2) ;Ly=b から y を求める (前進代入) (def y (forward_substitution 0 a2 b)) (println "Y") (disp_vector y) ;DL^Tx=y から x を求める (後退代入) (def x (backward_substitution (dec N) a2 y)) (println "X") (disp_vector (reverse x))
(def N 4) (def a [[-1.0 -2.0 7.0 -2.0] [1.0 -1.0 -2.0 6.0] [9.0 2.0 1.0 1.0] [2.0 8.0 -2.0 1.0]]) (def b [8.0 17.0 20.0 16.0]) ;1次元配列を表示 (defn disp_vector [row] (doseq [col row] (print (format "%14.10f\t" col))) (println)) ;2次元配列を表示 (defn disp_matrix [matrix] (doseq [row matrix] (doseq [col row] (print (format "%14.10f\t" col))) (println))) ;各列で 一番値が大きい行を 探す (defn get_max_row [row col a max_row max_val] ;一番値が大きい行 (def ws1 (if (< max_val (. Math abs (nth (nth a row) col))) (vector row (. Math abs (nth (nth a row) col))) (vector max_row max_val))) (def max_row2 (first ws1)) (def max_val2 (second ws1)) (if (>= row (dec (count a))) (vector max_row2 max_val2) (get_max_row (inc row) col a max_row2 max_val2))) ;ピボット選択 (defn pivoting [pivot a b a2 b2] (def ws1 (get_max_row 0 pivot a 0 0.0)) (def max_row (first ws1)) (def max_val (second ws1)) (def a3 (cons (nth a max_row) a2)) (def b3 (cons (nth b max_row) b2)) (def a4 (concat (take max_row a) (drop (inc max_row) a))) (def b4 (concat (take max_row b) (drop (inc max_row) b))) (if (>= pivot (dec N)) (vector (reverse a3) (reverse b3)) (pivoting (inc pivot) a4 b4 a3 b3))) ;前進消去 (defn forward_elim_loop [pivot row a b] (def s (/ (nth (nth a row) pivot) (nth (nth a pivot) pivot))) (def a_map (map (fn [a_piv a_row] (- a_row (* a_piv s))) (drop (inc pivot) (nth a pivot)) (drop (inc pivot) (nth a row)))) (def a1 (concat (take pivot (nth a row)) (cons s a_map))) (def a2 (concat (take row a) (cons (vec a1) (drop (inc row) a)))) (def x (- (nth b row) (* (nth b pivot) s))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (< row (dec N)) (forward_elim_loop pivot (inc row) a2 b) (vector a2 b))) (defn forward_elimination [pivot a b] (def ws1 (if (< pivot (dec N)) (forward_elim_loop pivot (inc pivot) a b) (vector a b))) (def a2 (first ws1)) (def b2 (second ws1)) (if (< pivot (dec N)) (forward_elimination (inc pivot) a2 b2) (vector a2 b2))) ;前進代入 (defn forward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (take row (nth a row)) b)) (def s (apply + a_map)) (def y (- (nth b row) s)) (def b2 (concat (take row b) (cons y (drop (inc row) b)))) (if (< row (dec N)) (cons y (forward_substitution (inc row) a b2)) (cons y []))) ;後退代入 (defn backward_substitution [row a b] (def a_map (map (fn [a_col b_col] (* a_col b_col)) (drop (inc row) (nth a row)) (drop (inc row) b))) (def s (apply + a_map)) (def x (/ (- (nth b row) s) (nth (nth a row) row))) (def b2 (concat (take row b) (cons x (drop (inc row) b)))) (if (> row 0) (cons x (backward_substitution (dec row) a b2)) (cons x []))) ;ピボット選択 (def ws1 (pivoting 0 a b [] [])) (def a1 (first ws1)) (def b1 (vec (second ws1))) (println "pivoting") (println "A") (disp_matrix a1) (println "B") (disp_vector b1) (println "") ;前進消去 (def ws2 (forward_elimination 0 a1 b1)) (def a2 (first ws2)) (def b2 (vec (second ws2))) (println "LU") (disp_matrix a2) ;Ly=b から y を求める (前進代入) (def y (forward_substitution 0 a2 b2)) (println "Y") (disp_vector (reverse y)) ;Ux=y から x を求める (後退代入) (def x (backward_substitution (dec N) a2 y)) (println "X") (disp_vector (reverse x))