(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))) |