Haskell で 積分(台形則)して π を求める

本日発売！

さまざまな言語で数値計算 第２巻 数値積分

πの求め方

積分(台形則)

f::Double->Double
f x = 4 / (1 + (x * x))

-- 台形則
trapezoidal::Integer->Integer->Double->Double->Double->Double
trapezoidal i n h x s
    | i <= (n - 1) = trapezoidal (i + 1) n h (x + h) (s + (f (x + h)))
    | otherwise    = s

import Text.Printf
import Control.Monad

forM_ ([1..10::Integer]) $ \j -> do
    let n = 2 ^ j
    let a = 0
    let b = 1
    let h = (fromIntegral (b - a)) / (fromIntegral n)
    let x = 0
    -- 台形則で積分
    let w1 = trapezoidal 1 n h x 0
    let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2
    let t1 = h * (w1 + w2)
    -- 結果を π と比較
    let t2 = t1 - pi
    printf "%3d : %13.10f, %13.10f\n" j t1 t2

  1 :  3.1000000000, -0.0415926536
  2 :  3.1311764706, -0.0104161830
  3 :  3.1389884945, -0.0026041591
  4 :  3.1409416120, -0.0006510415
  5 :  3.1414298932, -0.0001627604
  6 :  3.1415519635, -0.0000406901
  7 :  3.1415824811, -0.0000101725
  8 :  3.1415901105, -0.0000025431
  9 :  3.1415920178, -0.0000006358
 10 :  3.1415924946, -0.0000001589

こうしてやれば、再帰関数「trapezoidal」は不要だと、後で気づいた.

forM_ ([1..10::Integer]) $ \j -> do
    let n = 2 ^ j 
    let a = 0.0
    let b = 1.0
    let h = (b - a) / (fromIntegral n)
    -- 台形則で積分
    let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a) $ [1..(n - 1)]
    let w2 = ((f a) + (f b)) / 2
    let t1 = h * (w1 + w2)
    -- 結果を π と比較
    let t2 = t1 - pi
    printf "%3d : %13.10f, %13.10f\n" j t1 t2

自分で直接再帰を書くより fold や reduce を使え, と言うのは、こうゆうことか

参考文献