Hatena::ブログ(Diary)

ONLY DO WHAT ONLY YOU CAN DO このページをアンテナに追加 RSSフィード

2013.12.29 Clojureでラグランジュ補間

Clojureで関数の近似(ラグランジュ補間)

f:id:fornext1119:20131213103126p:image
f:id:fornext1119:20131213103127p:image
をラグランジュ補間で近似する

n+1個の点 (x0, y0), (x1, y1) … (xn, yn)
が与えられているとき, これらすべての点を通る n次式は次のように表すことができる.
f:id:fornext1119:20131207120958p:image
この式を使って, 与えられた点以外の点の値を求める.

(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))))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0701.clj
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667    -0.00000
-3.50   -0.73099    -0.73099    -0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964    -0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167    -0.00000
-0.50   -0.47943    -0.47943    -0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000
参考文献



スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/fornext1119/20131229/p3