数値積分

ずいぶん空いてしまいました。でも頑張るぞ。
積分というものは殆どの場合、解析的には不可能なので実際には数値的に行うことが多いです。今回はそんな数値積分の代表的(ていうか簡単な)方法として2種類、台形近似とシンプソン近似です。考え方は両方とも同じで、積分区間を適当に分割、それぞれの分割された領域でグラフを解析的に積分可能な関数で近似、それを計算して最後に足し合わせるわけです。積分可能な関数は、台形近似が1次関数で、シンプソンが２次関数です。まあ詳しいことは例によってwikipediaにでも聞いてください（＾＾；
まあそんなわけでちゃっちゃと実装しちゃいました。

main = do print( trapezoid 10 f 0.0 1.0 )
          print( simpson   10 f 0.0 1.0 )

-- 被積分関数
f x = 4.0 / ( 1.0 + x ** 2 )

-- 台形近似
trapezoid 1 f a b = ((f a) + (f b)) * (b-a) / 2
trapezoid n f a b = let m = (b+a)/2
                    in ( trapezoid (n-1) f a m ) + ( trapezoid (n-1) f m b )

-- シンプソン
simpson 1 f a b = ((f a) + 4*(f ((a+b)/2)) + (f b)) * (b-a) / 6
simpson n f a b = let m = (b+a)/2
                  in ( simpson (n-1) f a m ) + ( simpson (n-1) f m b )

特筆することは無いですね。定義どおり書いただけです。ただこの２つ、最初にも書いたように近似関数が違うだけで考え方は一緒です。なのでその共通する部分をくくり出せないかな？と考えたらくくり出せました（＾＾；

main = do print( divsum trapezoid 10 f 0.0 1.0 )
          print( divsum simpson   10 f 0.0 1.0 )

-- 分割・合計
divsum s 1 f a b = s f a b
divsum s n f a b = let m = (b+a)/2
                   in ( divsum s (n-1) f a m ) + ( divsum s (n-1) f m b )

-- 台形近似
trapezoid f a b = ((f a) + (f b)) * (b-a) / 2

-- シンプソン
simpson f a b = ((f a) + 4*(f ((a+b)/2)) + (f b)) * (b-a) / 6

というわけで、ちょうどいいリハビリでした。