import Text.Printf main = do print (3 + 5) print (3 - 5) print (3 * 5) print (3 ^ 5) print (5 / 3) print (5 `div` 3) print (div 5 3) print (mod 5 3) print (quot 5 3) print (rem 5 3) print (div (-5) 3) print (mod (-5) 3) print (quot (-5) 3) print (rem (-5) 3) putStrLn ((show (3 * 5))) printf "%3d\n" ((3 * 5)::Int) printf "%23.20f\n" ((5 / 3)::Double) import Text.Printf main = do let i = (3 * 5)::Int printf "3 * 5 = %d\n" i putStrLn ("3 * 5 = " ++ (show i)) import Text.Printf import Control.Monad main = do forM [1..9::Int] $ \i -> printf "%d, " i import Text.Printf import Control.Monad main = do forM [1..9::Int] $ \i -> do if (i `mod` 3 == 0) then printf "%d, " i else printf "" total 0 = 0 total i = if i `mod` 3 == 0 then i + total (i - 1) else total (i - 1) main = print (total 99) [1..9] filter (\x -> mod x 3 == 0) [1..9] foldl (+) 0 $ filter (\x -> mod x 3 == 0) [1..99] sum $ filter (\x -> mod x 3 == 0) [1..99] -- 初項:a, 公差:a で, 上限:lim の数列の総和を返す関数 sn::Int->Int->Int sn a lim = let n = lim `div` a -- 項数:n = 上限:lim / 公差:a l = a * n -- 末項:l = 項数:n * 公差:a in (a + l) * n `div` 2 -- 総和:sn = (初項:a + 末項:l) * 項数:n / 2 -- 3 の倍数の合計 main = putStrLn (show (sn 3 999)) let n = 10000 (n * (n + 1)) `div` 2 let n = 5000 n * (n + 1) let n = 5000 n ^ 2 let n = 1000 n * (n + 1) * (2 * n + 1) `div` 6 let n = 100 (n ^ 2) * ((n + 1) ^ 2) `div` 4 let n = 10 let a = 2 let r = 3 a * (r ^ n - 1) `div` (r -1) [0..9] map(\n -> n * 3 + 5)[0..9] product $ map(\n -> n * 3 + 5)[0..9] -- 等差数列の積 prod m d 0 = 1 prod m d n = m * (prod (m + d) d (n - 1)) -- 初項 5, 公差 3, 項数 10 の数列の積 prod 5 3 10 -- 階乗を求める関数 fact 0 = 1 fact n = n * fact (n-1) -- 10の階乗 fact 10 -- 下降階乗冪 fallingFact x 1 = x fallingFact x n = x * (fallingFact (x - 1) (n - 1)) -- 10 から 6 までの 総乗 fallingFact 10 5 -- 上昇階乗冪 risingFact x 1 = x risingFact x n = x * (risingFact (x + 1) (n - 1)) -- 10 から 14 までの 総乗 risingFact 10 5 -- 階乗 fact 0 = 1 fact n = n * fact (n-1) -- 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) let n = 10 let r = 5 (fact n) `div` (fact (n - r)) -- 下降階乗冪 fallingFact x 1 = x fallingFact x n = x * (fallingFact (x - 1) (n - 1)) -- 順列 (異なる 10 個のものから 5 個取ってできる順列の総数) let n = 10 let r = 5 fallingFact n r -- 重複順列 (異なる 10 個のものから重複を許して 5 個取ってできる順列の総数) let n = 10 let r = 5 n ^ r -- 組合せ comb n 0 = 1 comb n 1 = n comb n r | n == r = 1 | otherwise = (comb (n - 1) (r - 1)) + (comb (n - 1) r) -- 異なる 10 個のものから 5 個取ってできる組合せの総数 let n = 10 let r = 5 comb n r -- 組合せ comb n 0 = 1 comb n 1 = n comb n r | n == r = 1 | otherwise = (comb (n - 1) (r - 1)) + (comb (n - 1) r) -- 異なる 10 個のものから重複を許して 5 個取ってできる組合せの総数 let n = 10 let r = 5 comb (n + r - 1) r -- 自作の正弦関数 mySin::Double->Int->Bool->Double->Double->Double->Double mySin x n nega numerator denominator y = let m = 2 * n denom = denominator * (fromIntegral (m + 1)) * (fromIntegral m) num = numerator * x * x a = num / denom in -- 十分な精度になったら処理を抜ける if a <= 0.00000000001 then y else y + (mySin x (n + 1) (not nega) num denom (if nega then a else -a)) import Text.Printf import Control.Monad forM_ ( filter (\n -> (mod n 30 == 0 || mod n 45 == 0)) $ map (* 15) $ [0..24::Int] ) $ \degree -> do let radian = (fromIntegral degree) * pi / 180.0 -- 自作の正弦関数 let d1 = (mySin radian 1 False radian 1.0 radian) -- 標準の正弦関数 let d2 = sin(radian) -- 標準関数との差異 printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2) -- 自作の余弦関数 myCos::Double->Int->Bool->Double->Double->Double->Double myCos x n nega numerator denominator y = let m = 2 * n denom = denominator * (fromIntegral m) * (fromIntegral (m - 1)) num = numerator * x * x a = num / denom in -- 十分な精度になったら処理を抜ける if a <= 0.00000000001 then y else y + (myCos x (n + 1) (not nega) num denom (if nega then a else -a)) import Text.Printf import Control.Monad forM_ ( filter (\n -> (mod n 30 == 0 || mod n 45 == 0)) $ map (* 15) $ [0..24::Int] ) $ \degree -> do let radian = (fromIntegral degree) * pi / 180.0 -- 自作の余弦関数 let d1 = (myCos radian 1 False 1.0 1.0 1.0) -- 標準の余弦関数 let d2 = cos(radian) -- 標準関数との差異 printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2) -- 自作の正接関数 myTan::Double->Double->Int->Double->Double myTan x x2 n t = let denom = x2 / ((fromIntegral n) - t) num = n - 2 in if num <= 1 then x / (1.0 - denom) else myTan x x2 num denom import Text.Printf import Control.Monad forM_ ( map (\n -> n - 90) $ filter (\n -> mod n 180 /= 0) $ map (* 15) $ [0..12::Int] ) $ \degree -> do let radian = (fromIntegral degree) * pi / 180.0 -- 自作の正接関数 let x2 = radian * radian let d1 = (myTan radian x2 15 0.0) -- 15:必要な精度が得られる十分大きな奇数 -- 標準の正接関数 let d2 = tan(radian) -- 標準関数との差異 printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2) -- 自作の指数関数 myExp::Double->Int->Double->Double->Double->Double myExp x n numerator denominator y = let denom = denominator * (fromIntegral n) num = numerator * x a = num / denom in -- 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myExp x (n + 1) num denom a) import Text.Printf import Control.Monad forM_ ( map (\n -> (fromIntegral (n - 10)) / 4.0) $ [0..20::Int] ) $ \x -> do -- 標準の指数関数 let d1 = exp(x) -- 自作の指数関数 let d2 = (myExp x 1 1.0 1.0 1.0) -- 標準関数との差異 printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2) -- 自作の指数関数 myExp::Double->Double->Int->Double->Double; myExp x x2 n t = let denom = x2 / ((fromIntegral n) + t) num = n - 4 in if num < 6 then 1.0 + ((2.0 * x) / (2.0 - x + denom)) else myExp x x2 num denom import Text.Printf import Control.Monad forM_ ( map (\n -> (fromIntegral (n - 10)) / 4.0) $ [0..20::Int] ) $ \x -> do -- 標準の指数関数 let d1 = exp(x) -- 自作の指数関数 let x2 = x * x let d2 = (myExp x x2 30 0.0) -- 30:必要な精度が得られるよう, 6から始めて4ずつ増加させる -- 標準関数との差異 printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2) -- 自作の対数関数 myLog::Double->Double->Double->Double->Double; myLog x2 numerator denominator y = let denom = denominator + 2.0 num = numerator * x2 * x2 a = num / denom in -- 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myLog x2 num denom a) import Text.Printf import Control.Monad forM_ ( map (\n -> (fromIntegral n) / 5.0) $ [1..20::Int] ) $ \x -> do -- 標準の対数関数 let d1 = log(x) -- 自作の対数関数 let x2 = (x - 1.0) / (x + 1.0) let d2 = 2.0 * (myLog x2 x2 1.0 x2) -- 標準関数との差異 printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2) -- 自作の対数関数 myLog::Double->Int->Double->Double myLog x n t = let n2 = if rem n 2 == 0 then 2 else n n3 = if n > 3 then n2 else n x2 = if n > 3 then x * (fromIntegral (div n 2)) else x t2 = x2 / ((fromIntegral n3) + t) in if n <= 2 then x / (1.0 + t2) else myLog x (n - 1) t2 import Text.Printf import Control.Monad forM_ ( map (\n -> (fromIntegral n) / 5.0) $ [1..20::Int] ) $ \x -> do -- 標準の対数関数 let d1 = log(x) -- 自作の対数関数 let d2 = (myLog (x - 1.0) 27 0.0) -- 27:必要な精度が得られる十分大きな奇数 -- 標準関数との差異 printf "%5.2f : %13.10f - %13.10f = %13.10f\n" x d1 d2 (d1- d2) -- 自作の双曲線正弦関数 mySinh::Double->Int->Double->Double->Double->Double mySinh x n numerator denominator y = let m = 2 * n denom = denominator * (fromIntegral (m + 1)) * (fromIntegral m) num = numerator * x * x a = num / denom in -- 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (mySinh x (n + 1) num denom a) import Text.Printf import Control.Monad forM_ ( map ((-) 10) $ [0..20::Int] ) $ \i -> do -- 自作の双曲線正弦関数 let x = (fromIntegral i) let d1 = (mySinh x 1 x 1.0 x) -- 標準の双曲線正弦関数 let d2 = sinh(x) -- 標準関数との差異 printf "%3d : %17.10f - %17.10f = %13.10f\n" i d1 d2 (d1- d2) -- 自作の双曲線余弦関数 myCosh::Double->Int->Double->Double->Double->Double; myCosh x n numerator denominator y = let m = 2 * n denom = denominator * (fromIntegral m) * (fromIntegral (m - 1)) num = numerator * x * x a = num / denom in -- 十分な精度になったら処理を抜ける if abs(a) <= 0.00000000001 then y else y + (myCosh x (n + 1) num denom a) import Text.Printf import Control.Monad forM_ ( map ((-) 10) $ [0..20::Int] ) $ \i -> do -- 自作の双曲線余弦関数 let x = (fromIntegral i) let d1 = (myCosh x 1 1.0 1.0 1.0) -- 標準の双曲線余弦関数 let d2 = cosh(x) -- 標準関数との差異 printf "%3d : %17.10f - %17.10f = %13.10f\n" i d1 d2 (d1- d2) -- 自作の逆正接関数 myAtan::Double->Double->Int->Double->Double myAtan x x2 n t = let m = (fromIntegral (n `div` 2)) denom = (m * m * x2) / ((fromIntegral n) + t) num = n - 2 in if num <= 1 then x / (1.0 + denom) else myAtan x x2 num denom import Text.Printf import Control.Monad forM_ ( map (\n -> n - 45) $ map (* 15) $ [0..6::Int] ) $ \degree -> do let radian = (fromIntegral degree) * pi / 180.0 -- 自作の逆正接関数 let x2 = radian * radian let d1 = (myAtan radian x2 23 0.0) -- 23:必要な精度が得られる十分大きな奇数 -- 標準の逆正接関数 let d2 = atan(radian) -- 標準関数との差異 printf "%3d : %13.10f - %13.10f = %13.10f\n" degree d1 d2 (d1- d2) -- 自作の逆正接関数 myAtan::Double->Double->Int->Double->Double myAtan x x2 n t = let m = (fromIntegral (n `div` 2)) denom = (m * m * x2) / ((fromIntegral n) + t) num = n - 2 in if num <= 1 then x / (1.0 + denom) else myAtan x x2 num denom import Text.Printf import Control.Monad :{ forM_ ( map (\n -> n * 2 + 1) $ [5..15::Int] ) $ \i -> do -- 11..31 let radian = 1 -- 自作の逆正接関数 let x2 = radian * radian let d1 = (myAtan radian x2 i 0.0) -- i:必要な精度が得られる十分大きな奇数 -- 標準の逆正接関数 let d2 = atan(radian) -- 標準関数との差異 printf "%2d : %13.10f, %13.10f\n" i (d1 * 4) (d1 * 4 - pi) :} printf "%33.30f\n" pi 3.14159 26535 89793 3.14159 26535 89793 23846 26433 83279 50288 … import Text.Printf import Control.Monad f::Double->Double f x = 4 / (1 + (x * x)) main = do 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 import Text.Printf import Control.Monad f::Double->Double f x = 4 / (1 + (x * x)) main = do 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 a2 = a + (h / 2) let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a2) $ [0..(n-1)] let t1 = h * w1 -- 結果を π と比較 let t2 = t1 - pi printf "%3d : %13.10f, %13.10f\n" j t1 t2 import Text.Printf import Control.Monad f::Double->Double f x = 4 / (1 + (x * x)) main = do forM_ ([1..5::Integer]) $ \j -> do let n = 2 ^ j let a = 0.0 let b = 1.0 let h = (b - a) / (fromIntegral n) -- シンプソン則で積分 let w4 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2 - 1) $ [1..(div n 2)] let w2 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2) $ [1..(div n 2)] let t1 = (w2 - (f b)) * 2 + (f a) + (f b) let t2 = w4 * 4 let t3 = (t1 + t2) * h / 3 -- 結果を π と比較 let t4 = t3 - pi printf "%3d : %13.10f, %13.10f\n" j t3 t4 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 Control.Monad -- 台形則で積分して, 結果を保存 t <- forM ([1..6::Integer]) $ \j -> do let n = 2 ^ j let a = 0 let b = 1 let h = (fromIntegral (b - a)) / (fromIntegral n) let x = 0::Double let s = 0::Double let w1 = trapezoidal 1 n h x 0::Double let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2 let t1 = h * (w1 + w2) let t2 = t1 - pi return t1 -- Richardsonの補外法 rich_sub n [] ss = ss rich_sub n (x:[]) ss = ss rich_sub n (x:xs) ss = let x2 = (head xs) s = x2 + (x2 - x) / (fromIntegral (n - 1)) in rich_sub n xs (s:ss) richardson n ([]) = 0 richardson n (x:[]) = x richardson n (x:xs) = let s = rich_sub n (x:xs) [] in richardson (n * 4) (reverse s) -- Richardsonの補外法で積分して, 結果を π と比較 import Text.Printf forM_ ([2..6::Int]) $ \i -> do let s = richardson 4 (take i t) printf "%3d : %13.10f, %13.10f\n" i s (s - pi) :{ let{ 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 } :} :{ t <- forM ([1..6::Integer]) $ \j -> do let n = 2 ^ j let a = 0 let b = 1 let h = (fromIntegral (b - a)) / (fromIntegral n) let x = 0::Double let s = 0::Double let w1 = trapezoidal 1 n h x 0::Double let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2 let t1 = h * (w1 + w2) let t2 = t1 - pi return t1 :} Control.Monad> t [3.1,3.1311764705882354,3.138988494491089,3.140941612041389,3.141429893174975,3.141551963485654] :{ let { rich_sub n [] ss = ss; rich_sub n (x:[]) ss = ss; rich_sub n (x:xs) ss = let x2 = (head xs) s = x2 + (x2 - x) / (fromIntegral (n - 1)) in rich_sub n xs (s:ss) } :} :{ let { richardson n ([]) = 0; richardson n (x:[]) = x; richardson n (x:xs) = let s = rich_sub n (x:xs) [] in richardson (n * 4) (reverse s) } :} :{ forM_ ([2..6::Int]) $ \i -> do let s = richardson 4 (take i t) printf "%3d : %13.10f, %13.10f\n" i s (s - pi) :} /* 2 : 3.1415686275, -0.0000240261 3 : 3.1415940941, 0.0000014405 4 : 3.1415926384, -0.0000000152 5 : 3.1415926536, 0.0000000001 6 : 3.1415926536, -0.0000000000 */ midpoint C:\Users\Ifthen>ghci GHCi, version 7.6.3: http://www.haskell.org/ghc/ :? for help Loading package ghc-prim ... linking ... done. Loading package integer-gmp ... linking ... done. Loading package base ... linking ... done. Prelude> :{ Prelude| let { Prelude| f::Double->Double; Prelude| f x = 4 / (1 + (x * x)) Prelude| } Prelude| :} Prelude> :{ Prelude| let { Prelude| trape::Integer->Integer->Double->Double->Double->Double; Prelude| trape i n h x s Prelude| | i <= (n - 1) = trape (i + 1) n h (x + h) (s + (f (x + h))) Prelude| | otherwise = s Prelude| } Prelude| :} Prelude> import Text.Printf Prelude Text.Printf> import Control.Monad Prelude Text.Printf Control.Monad> :{ Prelude Text.Printf Control.Monad| forM_ ([1..10::Integer]) $ \j -> do Prelude Text.Printf Control.Monad| let n = 2 ^ j Prelude Text.Printf Control.Monad| let a = 0 Prelude Text.Printf Control.Monad| let b = 1 Prelude Text.Printf Control.Monad| let h = (fromIntegral (b - a)) / (fromIntegral n) Prelude Text.Printf Control.Monad| let x = 0::Double Prelude Text.Printf Control.Monad| let s = 0::Double Prelude Text.Printf Control.Monad| Prelude Text.Printf Control.Monad| let w1 = trape 1 n h x 0::Double Prelude Text.Printf Control.Monad| Prelude Text.Printf Control.Monad| let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2::Double Prelude Text.Printf Control.Monad| let t1 = (realToFrac h) * (w1 + w2) Prelude Text.Printf Control.Monad| let t2 = t1 - pi Prelude Text.Printf Control.Monad| printf "%3d : %23.20f, %23.20f\n" j t1 t2 Prelude Text.Printf Control.Monad| :} 1 : 3.10000000000000000000, -0.04159265358979303000 2 : 3.13117647058823540000, -0.01041618300155766600 3 : 3.13898849449108900000, -0.00260415909870426180 4 : 3.14094161204138900000, -0.00065104154840422980 5 : 3.14142989317497500000, -0.00016276041481821935 6 : 3.14155196348565400000, -0.00004069010413898511 7 : 3.14158248106375200000, -0.00001017252604107455 8 : 3.14159011045828150000, -0.00000254313151160090 9 : 3.14159201780692140000, -0.00000063578287168298 10 : 3.14159249464407250000, -0.00000015894572058528 Prelude Text.Printf Control.Monad> :{ Prelude Text.Printf Control.Monad| t <- forM ([1..10::Integer]) $ \j -> do Prelude Text.Printf Control.Monad| let n = 2 ^ j Prelude Text.Printf Control.Monad| let a = 0 Prelude Text.Printf Control.Monad| let b = 1 Prelude Text.Printf Control.Monad| let h = (fromIntegral (b - a)) / (fromIntegral n) Prelude Text.Printf Control.Monad| let x = 0::Double Prelude Text.Printf Control.Monad| let s = 0::Double Prelude Text.Printf Control.Monad| Prelude Text.Printf Control.Monad| let w1 = trape 1 n h x 0::Double Prelude Text.Printf Control.Monad| Prelude Text.Printf Control.Monad| let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2 Prelude Text.Printf Control.Monad| let t1 = h * (w1 + w2) Prelude Text.Printf Control.Monad| let t2 = t1 - pi Prelude Text.Printf Control.Monad| return t1 Prelude Text.Printf Control.Monad| :} Prelude Text.Printf Control.Monad> mapM_ print t 3.1 3.1311764705882354 3.138988494491089 3.140941612041389 3.141429893174975 3.141551963485654 3.141582481063752 3.1415901104582815 3.1415920178069214 3.1415924946440725 Prelude Text.Printf Control.Monad> forM_ t print 3.1 3.1311764705882354 3.138988494491089 3.140941612041389 3.141429893174975 3.141551963485654 3.141582481063752 3.1415901104582815 3.1415920178069214 3.1415924946440725 Prelude Text.Printf Control.Monad> :{ Prelude Text.Printf Control.Monad| let{ Prelude Text.Printf Control.Monad| rich n [] ss = ss; Prelude Text.Printf Control.Monad| rich n (x:[]) ss = ss; Prelude Text.Printf Control.Monad| rich n (x:xs) ss = Prelude Text.Printf Control.Monad| let Prelude Text.Printf Control.Monad| x2 = (head xs) Prelude Text.Printf Control.Monad| s = x2 + (x2- x) / (fromIntegral (n - 1)) Prelude Text.Printf Control.Monad| in Prelude Text.Printf Control.Monad| rich n xs (s:ss) Prelude Text.Printf Control.Monad| } Prelude Text.Printf Control.Monad| :} Prelude Text.Printf Control.Monad> :{ Prelude Text.Printf Control.Monad| let{ Prelude Text.Printf Control.Monad| call_rich n ([]) = 0; Prelude Text.Printf Control.Monad| call_rich n (x:[]) = x; Prelude Text.Printf Control.Monad| call_rich n (x:xs) = Prelude Text.Printf Control.Monad| let Prelude Text.Printf Control.Monad| s = rich n (x:xs) [] Prelude Text.Printf Control.Monad| in Prelude Text.Printf Control.Monad| call_rich (n * 4) (reverse s) Prelude Text.Printf Control.Monad| } Prelude Text.Printf Control.Monad| :} Prelude Text.Printf Control.Monad> :{ Prelude Text.Printf Control.Monad| forM_ ([2..6::Int]) $ \i -> do Prelude Text.Printf Control.Monad| let s = call_rich 4 (take i t) Prelude Text.Printf Control.Monad| printf "%3d : %13.10f, %13.10f\n" i s (s - pi) Prelude Text.Printf Control.Monad| :} 2 : 3.1415686275, -0.0000240261 3 : 3.1415940941, 0.0000014405 4 : 3.1415926384, -0.0000000152 5 : 3.1415926536, 0.0000000001 6 : 3.1415926536, -0.0000000000 Prelude Text.Printf Control.Monad> :q Leaving GHCi. C:\Users\Ifthen> */ import Text.Printf import Control.Monad :{ let { f::Double->Double; f x = 4 / (1 + (x * x)) } :} :{ let { trape::Integer->Integer->Double->Double->Double->Double; trape i n h x s | i <= (n - 1) = trape (i + 1) n h (x + h) (s + (f (x + h))) | otherwise = s } :} :{ forM_ ([1..6::Integer]) $ \j -> do let n = 2 ^ j let a = 0 let b = 1 let h = (fromIntegral (b - a)) / (fromIntegral n) let x = 0::Double let s = 0::Double let w1 = trape 1 n h x 0::Double let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2::Double let t1 = (realToFrac h) * (w1 + w2) let t2 = t1 - pi printf "%3d : %23.20f, %23.20f\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 :{ t <- 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::Double let s = 0::Double let w1 = trape 1 n h x 0::Double let w2 = ((f (fromIntegral a)) + (f (fromIntegral b))) / 2 let t1 = h * (w1 + w2) let t2 = t1 - pi return t1 :} mapM_ print t forM_ t print :{ let{ rich n [] ss = ss; rich n (x:[]) ss = ss; rich n (x:xs) ss = let x2 = (head xs) s = x2 + (x2- x) / (fromIntegral (n - 1)) in rich n xs (s:ss) } :} let s = test1 4 t [] mapM_ print s let s = test1 4 [3.1,3.13] [] mapM_ print s let s = test1 4 [3.1] [] mapM_ print s let s = test1 4 [] [] mapM_ print s :{ let{ call_rich n ([]) = 0; call_rich n (x:[]) = x; call_rich n (x:xs) = let s = rich n (x:xs) [] in call_rich (n * 4) (reverse s) } :} let s = test2 4 [] s let s = test2 4 [3.1] s let s = test2 4 [3.1,3.13] s let s = test2 4 [3.1,3.13,3.1311] s let s = test2 4 t s :{ forM_ ([2..6::Int]) $ \i -> do let s = call_rich 4 (take i t) printf "%3d : %13.10f, %13.10f\n" i s (s - pi) :} let x1 = (head t) let x2 = (head (tail t)) let t2_2 = x2 + ((x2- x1) / (fromIntegral (4 - 1))) printf "%13.10f\n" t2_2 let x1 = (head (tail t)) let x2 = (head (tail (tail t))) let t2_3 = x2 + ((x2- x1) / (fromIntegral (4 - 1))) printf "%13.10f\n" t2_3 let x1 = (head (tail (tail t))) let x2 = (head (tail (tail (tail t)))) let t2_4 = x2 + ((x2- x1) / (fromIntegral (4 - 1))) printf "%13.10f\n" t2_4 let x1 = (head (tail (tail (tail t)))) let x2 = (head (tail (tail (tail (tail t))))) let t2_5 = x2 + ((x2- x1) / (fromIntegral (4 - 1))) printf "%13.10f\n" t2_5 let x1 = (head (tail (tail (tail (tail t))))) let x2 = (head (tail (tail (tail (tail (tail t)))))) let t2_6 = x2 + ((x2- x1) / (fromIntegral (4 - 1))) printf "%13.10f\n" t2_6 let t3_3 = t2_3 + ((t2_3- t2_2) / (fromIntegral (16 - 1))) printf "%13.10f\n" t3_3 let s = test1 4 t [] mapM_ print s in test1 n xs (s:ss) :{ forM_ ([t_list]) $ \t -> do print (show t) :} :{ forM_ ([1..10::Integer]) $ \j -> do n <- return (2 ^ j) a <- return 0 b <- return 1 h <- return ((fromIntegral (b - a)) / (fromIntegral n)) x <- return (0::Double) s <- return (0::Double) w1 <- return (trape_sub 1 n h x 0::Double) w2 <- return (((f (fromIntegral a)) + (f (fromIntegral b))) / 2) t1 <- return (h * (w1 + w2)) t2 <- return (t1 - pi) printf "%3d : %13.10f, %13.10f\n" j t1 t2 :} :{ forM_ ( [1..10::Integer] ) $ \j -> do let n = 2 ^ j let a = 0 let b = 1 let h = (b - a) / n let w = trape_sub (toInteger n) (fromIntegral h) 1 (fromIntegral a) -- printf "%3d : %13.10f\n" j t1 printf "%3d\n" j :} let { trape::Double->Double->Integer->Double; trape a b j = let n = 2 ^ j h = (b - a) / n in if j <= 10 then (trape_sub (toInteger n) (fromIntegral h) 1 (fromIntegral a) 0) then trape a b (j + 1) else 0 } :} import Text.Printf import Control.Monad f::Double->Double f x = 4 / (1 + (x * x)) -- Richardsonの補外法 rich_sub n [] ss = ss rich_sub n (x:[]) ss = ss rich_sub n (x:xs) ss = let x2 = (head xs) s = x2 + (x2 - x) / (fromIntegral (n - 1)) in rich_sub n xs (s:ss) richardson n ([]) = 0 richardson n (x:[]) = x richardson n (x:xs) = let s = rich_sub n (x:xs) [] in richardson (n * 4) (reverse s) main = do let a = 0.0 let b = 1.0 t <- forM ([1..6::Integer]) $ \j -> do let n = 2 ^ j 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) return t1 -- Richardsonの補外法で積分して, 結果を π と比較 forM_ ([2..6::Int]) $ \i -> do let s = richardson 4 (take i t) printf "%3d : %13.10f, %13.10f\n" i s (s - pi) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- Lagrange (ラグランジュ) 補間 lagrange::Double->[Double]->[Double]->Double lagrange d x y = let sum_list = map(\i -> do let prod_list = map(\j -> do (d - x!!(fromIntegral j)) / (x!!(fromIntegral i) - x!!(fromIntegral j))) $ filter (\j -> j /= i) [0..n] product $ y!!(fromIntegral i) : prod_list ) [0..n] in sum $ sum_list main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (lagrange i x y)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- Lagrange (ラグランジュ) 補間 lagrange::Double->[Double]->[Double]->Double lagrange d x y = let sum_list = map(\i -> do let prod_list = map(\j -> do (d - x!!j) / (x!!i - x!!j) ) $ filter (\j -> j /= i) [0..n::Int] product $ y!!i : prod_list ) [0..n::Int] in sum $ sum_list main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (lagrange i x y)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- Neville (ネヴィル) 補間 neville::Double->[Double]->[Double]->Int->Double neville d x w j = let t = map(\i -> do (w!!(i+1) + (w!!(i+1) - w!!i) * (d - x!!(i+j)) / (x!!(i+j) - x!!i)) ) [0..(n-j)::Int] in if j == n then t!!0 else (neville d x t (j + 1)) main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (neville i x y 1)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- Newton (ニュートン) 補間 newton::Double->[Double]->[Double]->Double newton d x a = let sum_list = map(\i -> do let prod_list = map(\j -> do d - x!!j ) $ [0..(i-1)::Int] product $ a!!i : prod_list ) [1..n::Int] in sum $ a!!0 : sum_list -- 差分商の表を作る make_table::[Double]->[Double]->[Double]->Int->[Double] make_table x d a i = let w = map(\j -> do ((d!!(j+1) - d!!j) / (x!!(j+i) - x!!j)) ) $ [0..(n-i)::Int] t = w!!0:a in -- n階差分商 if i == n then t else (make_table x w t (i + 1)) main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x -- 差分商の表を作る let a = reverse (make_table x y [y!!0] 1) -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (newton i x a)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int nx2 = 13 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- 導関数 fd::Double->Double fd x = 1.0 - (x ^ 2) / (fromIntegral 2) + (x ^ 4) / (fromIntegral (4 * 3 * 2)) -- Hermite (エルミート) 補間 hermite::Double->[Double]->[Double]->Double hermite d z a = let sum_list = map(\i -> do let prod_list = map(\j -> do d - z!!j ) $ [0..(i-1)::Int] product $ a!!i : prod_list ) [1..nx2::Int] in sum $ a!!0 : sum_list -- 差分商の表を作る make_table::[Double]->[Double]->[Double]->[Double]->Int->[Double] make_table yd z d a i = let w = map(\j -> do if (i == 1 && (rem j 2) == 0) then yd!!(div j 2) else ((d!!(j+1) - d!!j) / (z!!(j+i) - z!!j)) ) $ [0..(nx2-i)::Int] t = w!!0:a in -- n階差分商 if i == nx2 then t else (make_table yd z w t (i + 1)) main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x let yd = map(\i -> fd(i)) x -- 差分商の表を作る let z = map(\i -> x!!(div i 2)) [0..nx2] let w = map(\i -> y!!(div i 2)) [0..nx2] let a = reverse (make_table yd z w [w!!0] 1) -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (hermite i z a)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf import Control.Monad -- データ点の数 - 1 n = 6 :: Int -- 元の関数 f::Double->Double f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2)) -- 3項方程式を解く (ト−マス法) thomas::Int->[Double]->[Double]->[Double]->[Double]->[Double]->[Double]->([Double], [Double]) thomas i a b c d gs ss = let g = b!!i - a!!i * c!!(i-1) / gs!!0 s = d!!i - a!!i * ss!!0 / gs!!0 in if i == (n - 1) then (reverse (g:gs), reverse (s:ss)) else (thomas (i + 1) a b c d (g:gs) (s:ss)) thomas2::Int->[Double]->[Double]->[Double]->[Double]->[Double] thomas2 i g s c zs = let z = (s!!i - c!!i * zs!!0) / g!!i in if i == 1 then 0.0 : z : zs else thomas2 (i - 1) g s c (z:zs) -- Spline (スプライン) 補間 spline::Double->[Double]->[Double]->[Double]->Double spline d x y z = let -- 補間関数値がどの区間にあるか k = if d <= x!!0 then 0 else (length $ filter (\i -> i < d) $ x) - 1 d1 = x!!(k+1) - d d2 = d - x!!k d3 = x!!(k+1) - x!!k in (z!!k * (d1 ^ 3) + z!!(k+1) * (d2 ^ 3)) / (6.0 * d3) + (y!!k / d3 - z!!k * d3 / 6.0) * d1 + (y!!(k+1) / d3 - z!!(k+1) * d3 / 6.0) * d2 main = do -- 1.5刻みで -4.5〜4.5 まで, 7点だけ値をセット let x = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n] let y = map(\i -> f(i)) x -- 3項方程式の係数の表を作る let a = 0.0 : map(\i -> x!!i - x!!(i-1)) [1..n - 1] let b = 0.0 : map(\i -> 2 * (x!!(i+1) - x!!(i-1))) [1..n - 1] let c = 0.0 : map(\i -> x!!(i+1) - x!!i) [1..n - 1] let d = 0.0 : map(\i -> 6 * ((y!!(i+1) - y!!i) / (x!!(i+1) - x!!i) - (y!!i - y!!(i-1)) / (x!!i - x!!(i-1)))) [1..n - 1] -- 3項方程式を解く (ト−マス法) let (g, s) = (thomas 2 a b c d [b!!1, 0.0] [d!!1, 0.0]) let z = (thomas2 (n - 2) g s c [s!!(n-1) / g!!(n-1), 0.0]) -- 0.5刻みで 与えられていない値を補間 let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18] let d2 = map(\i -> (f i)) d1 let d3 = map(\i -> (spline i x y z)) d1 forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3) import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Euler法 euler::Integer->Double->Double->Double->Double->IO () euler i vx vy x y = let t = (fromIntegral i) * h wx = x + h * vx wy = y + h * vy wvx = vx + h * (fx vx vy) wvy = vy + h * (fy vx vy) in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t vx vy wx wy euler (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Euler法 euler 1 vx vy x y import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Heun法 heun::Integer->Double->Double->Double->Double->IO () heun i vx vy x y = let t = (fromIntegral i) * h wx2 = x + h * vx wy2 = y + h * vy wvx2 = vx + h * (fx vx vy) wvy2 = vy + h * (fy vx vy) wx = x + h * ( vx + wvx2 ) / 2.0 wy = y + h * ( vy + wvy2 ) / 2.0 wvx = vx + h * ((fx vx vy) + (fx wvx2 wvy2)) / 2.0 wvy = vy + h * ((fy vx vy) + (fy wvx2 wvy2)) / 2.0 in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy heun (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Heun法 heun 1 vx vy x y import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- 中点法 midpoint::Integer->Double->Double->Double->Double->IO () midpoint i vx vy x y = let -- 経過秒数 t = (fromIntegral i) * h -- 位置・速度 wvx1 = h * (fx vx vy) wvy1 = h * (fy vx vy) wvx2 = vx + wvx1 / 2.0 wvy2 = vy + wvy1 / 2.0 wvx = vx + h * (fx wvx2 wvy2) wvy = vy + h * (fy wvx2 wvy2) wx = x + h * wvx2 wy = y + h * wvy2 in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy midpoint (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- 中点法 midpoint 1 vx vy x y import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Runge-Kutta法 rungekutta::Integer->Double->Double->Double->Double->IO () rungekutta i vx vy x y = let -- 経過秒数 t = (fromIntegral i) * h -- 位置・速度 wx1 = h * vx wy1 = h * vy wvx1 = h * (fx vx vy) wvy1 = h * (fy vx vy) wvx5 = vx + wvx1 / 2.0 wvy5 = vy + wvy1 / 2.0 wx2 = h * wvx5 wy2 = h * wvy5 wvx2 = h * (fx wvx5 wvy5) wvy2 = h * (fy wvx5 wvy5) wvx6 = vx + wvx2 / 2.0 wvy6 = vy + wvy2 / 2.0 wx3 = h * wvx6 wy3 = h * wvy6 wvx3 = h * (fx wvx6 wvy6) wvy3 = h * (fy wvx6 wvy6) wvx7 = vx + wvx3 wvy7 = vy + wvy3 wx4 = h * wvx7 wy4 = h * wvy7 wvx4 = h * (fx wvx7 wvy7) wvy4 = h * (fy wvx7 wvy7) wx = x + ( wx1 + wx2 * 2.0 + wx3 * 2.0 + wx4) / 6.0 wy = y + ( wy1 + wy2 * 2.0 + wy3 * 2.0 + wy4) / 6.0 wvx = vx + (wvx1 + wvx2 * 2.0 + wvx3 * 2.0 + wvx4) / 6.0 wvy = vy + (wvy1 + wvy2 * 2.0 + wvy3 * 2.0 + wvy4) / 6.0 in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy rungekutta (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Runge-Kutta法 rungekutta 1 vx vy x y import Text.Printf -- 重力加速度 g = -9.8 :: Double -- 空気抵抗係数 k = -0.01 :: Double -- 時間間隔(秒) h = 0.01 :: Double -- 空気抵抗による水平方向の減速分 fx::Double->Double->Double fx vx vy = let v = sqrt(vx * vx + vy * vy) in k * v * vx -- 空気抵抗による鉛直方向の減速分 fy::Double->Double->Double fy vx vy = let v = sqrt(vx * vx + vy * vy) in g + (k * v * vy) -- Runge-Kutta-Gill法 rungekuttagill::Integer->Double->Double->Double->Double->IO () rungekuttagill i vx vy x y = let -- 経過秒数 t = (fromIntegral i) * h -- 位置・速度 wx1 = h * vx wy1 = h * vy wvx1 = h * (fx vx vy) wvy1 = h * (fy vx vy) wvx5 = vx + wvx1 / 2.0 wvy5 = vy + wvy1 / 2.0 wx2 = h * wvx5 wy2 = h * wvy5 wvx2 = h * (fx wvx5 wvy5) wvy2 = h * (fy wvx5 wvy5) wvx6 = vx + wvx1 * ((sqrt(2.0) - 1.0) / 2.0) + wvx2 * (1.0 - 1.0 / sqrt(2.0)) wvy6 = vy + wvy1 * ((sqrt(2.0) - 1.0) / 2.0) + wvy2 * (1.0 - 1.0 / sqrt(2.0)) wx3 = h * wvx6 wy3 = h * wvy6 wvx3 = h * (fx wvx6 wvy6) wvy3 = h * (fy wvx6 wvy6) wvx7 = vx - wvx2 / sqrt(2.0) + wvx3 * (1.0 + 1.0 / sqrt(2.0)) wvy7 = vy - wvy2 / sqrt(2.0) + wvy3 * (1.0 + 1.0 / sqrt(2.0)) wx4 = h * wvx7 wy4 = h * wvy7 wvx4 = h * (fx wvx7 wvy7) wvy4 = h * (fy wvx7 wvy7) wx = x + ( wx1 + wx2 * (2.0 - sqrt(2.0)) + wx3 * (2.0 + sqrt(2.0)) + wx4) / 6.0 wy = y + ( wy1 + wy2 * (2.0 - sqrt(2.0)) + wy3 * (2.0 + sqrt(2.0)) + wy4) / 6.0 wvx = vx + (wvx1 + wvx2 * (2.0 - sqrt(2.0)) + wvx3 * (2.0 + sqrt(2.0)) + wvx4) / 6.0 wvy = vy + (wvy1 + wvy2 * (2.0 - sqrt(2.0)) + wvy3 * (2.0 + sqrt(2.0)) + wvy4) / 6.0 in if y >= 0.0 then do printf "%4.2f\t%8.5f\t%9.5f\t%9.5f\t%8.5f\n" t wvx wvy wx wy rungekuttagill (i+1) wvx wvy wx wy else return () main = do -- 角度 let degree = 45.0 :: Double let radian = degree * pi / 180.0 -- 初速 250 km/h -> 秒速に変換 let v = (fromIntegral (250 * 1000 `div` 3600)) -- 水平方向の速度 let vx = v * cos(radian) -- 鉛直方向の速度 let vy = v * sin(radian) -- 位置 let x = 0.0 let y = 0.0 -- Runge-Kutta-Gill法 rungekuttagill 1 vx vy x y import Text.Printf import Debug.Trace f::Double->Double f x = x * x - 2.0 bisection::Double->Double->Double bisection a b = let -- 区間 (a, b) の中点 c = (a + b) / 2 c = (a + b) / 2.0 fc = (f c) :: Double in if abs(fc) < 0.0000000001 then c else trace (printf "%12.10f\t%13.10f" c (c - ((sqrt 2.0)::Double))) $ if fc < 0.0 then -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在 bisection c b else -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在 bisection a c main = do let a = 1.0 let b = 2.0 printf "%12.10f\n" (bisection a b) import Text.Printf import Debug.Trace f::Double->Double f x = x * x - 2.0 falseposition::Double->Double->Double falseposition a b = let -- 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c = (a * (f b) - b * (f a)) / ((f b) - (f a)) fc = (f c) :: Double in if abs(fc) < 0.0000000001 then c else trace (printf "%12.10f\t%13.10f" c (c - ((sqrt 2.0)::Double))) $ if fc < 0.0 then -- f(c) < 0 であれば, 解は区間 (c, b) の中に存在 falseposition c b else -- f(c) > 0 であれば, 解は区間 (a, c) の中に存在 falseposition a c main = do let a = 1.0 let b = 2.0 printf "%12.10f\n" (falseposition a b) import Text.Printf import Debug.Trace g::Double->Double g x = (x / 2.0) + (1.0 / x) iterative::Double->Double iterative x0 = let x1 = (g x0) in if abs(x1 - x0) < 0.0000000001 then x1 else trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double))) iterative x1 main = do let x = 1.0 printf "%12.10f\n" (iterative x) import Text.Printf import Debug.Trace f0::Double->Double f0 x = x * x - 2.0 f1::Double->Double f1 x = 2.0 * x newton::Double->Double newton x0 = let x1 = x0 - ((f0 x0) / (f1 x0)) in if abs(x1 - x0) < 0.0000000001 then x1 else trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double))) newton x1 main = do let x = 2.0 printf "%12.10f\n" (newton x) import Text.Printf import Debug.Trace f0::Double->Double f0 x = x * x - 2.0 f1::Double->Double f1 x = 2.0 * x f2::Double->Double f2 x = 2.0 bailey::Double->Double bailey x0 = let x1 = x0 - ((f0 x0) / ((f1 x0) - ((f0 x0) * (f2 x0) / (2.0 * (f1 x0))))) in if abs(x1 - x0) < 0.0000000001 then x1 else trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double))) bailey x1 main = do let x = 2.0 printf "%12.10f\n" (bailey x) import Text.Printf import Debug.Trace f::Double->Double f x = x * x - 2.0 secant::Double->Double->Double secant x0 x1 = let x2 = x1 - (f x1) * (x1 - x0) / ((f x1) - (f x0)) in if abs(x2 - x1) < 0.0000000001 then x2 else trace (printf "%12.10f\t%13.10f" x2 (x2 - ((sqrt 2.0)::Double))) secant x1 x2 main = do let x0 = 1.0 let x1 = 2.0 printf "%12.10f\n" (secant x0 x1) import Data.Array import Text.Printf import Debug.Trace import Control.Monad n = 3 :: Int jacobi::Bool->[[Double]]->[Double]->[Double]->IO() jacobi finish a b x0 = do -- let -- x1 = array (0,3) [(0,0.0),(1,0.0),(2,0.0),(3,0.0)] -- in forM_ [0..3] $ \row -> do let a1 = zip (a!!row) x0 let a2 = take row a1 ++ drop (row+1) a1 let a3 = map (\(x, y) -> x * y) $ a2 let s = sum a3 printf "%f\t" s let x1 = ((b!!row) - s) / a!!row!!row printf "%f\n" x1 --x1 // [((fromIntegral row), ((b!!row) - s) / a!!row!!row)] -- printf "%d" (length $ filter (< 0) $ map (\(x,y) -> x - y) $ zip x0 (elems x1)) --putStrLn "" if finish then putStrLn "" else jacobi True a b x0 disp_progress::[[Double]]->[Double]->IO() disp_progress a b = do forM ([0..n]) $ \i -> do forM ([0..n]) $ \j -> do printf "%14.10f\t" (a!!i!!j) printf "%14.10f\t" (b!!i) putStrLn "" putStrLn "" disp_result::[Double]->IO() disp_result a = do forM ([0..n]) $ \i -> do printf "%14.10f\t" (a!!i) putStrLn "" main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17::Double] let x0 = [0,0,0,0::Double] disp_progress a b putStrLn "" jacobi False a b x0 disp_progress a b putStrLn "" -- let a1 = map (\pivot -> map(\x -> x / (a!!pivot!!pivot)) (a!!pivot)) [0..n] -- let b1 = map (\pivot -> b!!pivot / (a!!pivot!!pivot) ) [0..n] --disp_progress a b1 --disp_result b1 putStrLn "" -- while (true) -- { -- double[] x1 = new double[N]; -- bool finish = true; -- for (int i = 0; i < N; i++) -- { -- x1[i] = 0; -- for (int j = 0; j < N; j++) -- if (j != i) -- map (\x -> filter(/= x)[0..3])[0..3] -- x1[i] += a[i,j] * x0[j]; -- -- x1[i] = (b[i] - x1[i]) / a[i,i]; -- if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false; -- } -- -- for (int i = 0; i < N; i++) -- { -- x0[i] = x1[i]; -- Console.Write(string.Format("{0,14:F10}\t", x0[i])); -- } -- Console.WriteLine(); -- if (finish) return; -- } -- while (true) -- { -- double[] x1 = new double[N]; -- bool finish = true; -- for (int i = 0; i < N; i++) -- { -- x1[i] = 0; -- for (int j = 0; j < N; j++) -- if (j != i) -- x1[i] += a[i,j] * x0[j]; -- -- x1[i] = (b[i] - x1[i]) / a[i,i]; -- if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false; -- } -- -- for (int i = 0; i < N; i++) -- { -- x0[i] = x1[i]; -- Console.Write(string.Format("{0,14:F10}\t", x0[i])); -- } -- Console.WriteLine(); -- if (finish) return; -- } import Data.Array import Text.Printf import Debug.Trace import Control.Monad n = 3 :: Int disp_progress::[[Double]]->[Double]->IO() disp_progress a b = do forM ([0..n]) $ \i -> do forM ([0..n]) $ \j -> do printf "%14.10f\t" (a!!i!!j) printf "%14.10f\t" (b!!i) putStrLn "" putStrLn "" disp_result::[Double]->IO() disp_result a = do forM ([0..n]) $ \i -> do printf "%14.10f\t" (a!!i) putStrLn "" row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double] row_loop row a b x0 x1 = let a1 = zip (a!!row) x0 a2 = take row a1 ++ drop (row+1) a1 a3 = map (\(x, y) -> x * y) $ a2 s = sum a3 x = ((b!!row) - s) / a!!row!!row in if row >= n then x:x1 else (row_loop (row+1) a b x0 (x:x1)) main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17::Double] let x0 = [0,0,0,0::Double] disp_progress a b putStrLn "" let x1 = row_loop 0 a b x0 [] disp_result x1 putStrLn "" -- while (true) -- { -- double[] x1 = new double[N]; -- bool finish = true; -- for (int i = 0; i < N; i++) -- { -- x1[i] = 0; -- for (int j = 0; j < N; j++) -- if (j != i) -- x1[i] += a[i,j] * x0[j]; -- -- x1[i] = (b[i] - x1[i]) / a[i,i]; -- if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false; -- } -- -- for (int i = 0; i < N; i++) -- { -- x0[i] = x1[i]; -- Console.Write(string.Format("{0,14:F10}\t", x0[i])); -- } -- Console.WriteLine(); -- if (finish) return; -- } import Text.Printf import Debug.Trace import Control.Monad n = 3::Int disp_progress::[[Double]]->[Double]->IO() disp_progress a b = do forM_ ([0..n]) $ \i -> do forM_ ([0..n]) $ \j -> do printf "%14.10f\t" (a!!i!!j) printf "%14.10f\t" (b!!i) putStrLn "" putStrLn "" disp_result::[Double]->IO() disp_result a = do forM_ ([0..n]) $ \i -> do printf "%14.10f\t" (a!!i) putStrLn "" row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double] row_loop row a b x0 x1 = let a1 = zip (a!!row) x0 a2 = take row a1 ++ drop (row+1) a1 a3 = map (\(x, y) -> x * y) $ a2 s = sum a3 x = ((b!!row) - s) / a!!row!!row in if row <= 0 then x:x1 else (row_loop (row-1) a b x0 (x:x1)) jacobi::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double]) jacobi a b x0 xs = let x1 = (row_loop n a b x0 []) cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1 in if cnt < 1 then (reverse xs, x1) else (jacobi a b x1 (x1:xs)) main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17::Double] let x0 = [0,0,0,0::Double] disp_progress a b putStrLn "" let (xs, x1) = jacobi a b x0 [] putStrLn "mapM_" mapM_ (\x -> (printf "%14.10f\t" x)) $ x1 putStrLn "" putStrLn "forM_1" forM_ x1 $ \x -> (printf "%14.10f\t" x) putStrLn "" putStrLn "forM_2" forM_ x1 $ \x -> do (printf "%14.10f\t" x) putStrLn "" putStrLn "forM_3" forM_ xs $ \x -> do forM_ x $ \x1 -> do (printf "%14.10f\t" x1) putStrLn "" putStrLn "" mapM_ (\x -> (printf "%14.10f\t" x)) $ x1 import Text.Printf import Debug.Trace import Control.Monad n = 3 :: Int disp_progress::[[Double]]->IO() disp_progress a = do forM ([0..n]) $ \i -> do forM ([0..n]) $ \j -> do printf "%14.10f\t" (a!!i!!j) putStrLn "" putStrLn "" --return () jacobi::Bool->[[Double]]->[Double]->([[Double]], [Double]) jacobi finish a b = -- do if finish then (a, b) else --trace (printf "%14.10f", 1.0 ) $ jacobi True a b jacobi2::Bool->[[Double]]->[Double]->IO() jacobi2 finish a b = -- do if finish then putStrLn "" else jacobi2 True a b test::(Show a)=>a->a test x = trace (show x) $ x main_loop::Bool->[[Double]]->[Double]->IO() main_loop finish a b = -- do if finish then putStrLn "" else -- double[] x1 = new double[N]; -- bool finish = true; -- for (int i = 0; i < N; i++) -- { -- x1[i] = 0; -- for (int j = 0; j < N; j++) -- if (j != i) -- x1[i] += a[i,j] * x0[j]; -- -- x1[i] = (b[i] - x1[i]) / a[i,i]; -- if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false; -- } -- -- for (int i = 0; i < N; i++) -- { -- x0[i] = x1[i]; -- Console.Write(string.Format("{0,14:F10}\t", x0[i])); -- } main_loop True a b --main :: IO () main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17] disp_progress a let (x, y) = jacobi False a b printf "%14.10f\n" (x!!0!!0) main_loop False a b -- while (true) -- { -- double[] x1 = new double[N]; -- bool finish = true; -- for (int i = 0; i < N; i++) -- { -- x1[i] = 0; -- for (int j = 0; j < N; j++) -- if (j != i) -- x1[i] += a[i,j] * x0[j]; -- -- x1[i] = (b[i] - x1[i]) / a[i,i]; -- if (Math.Abs(x1[i] - x0[i]) > 0.0000000001) finish = false; -- } -- -- for (int i = 0; i < N; i++) -- { -- x0[i] = x1[i]; -- Console.Write(string.Format("{0,14:F10}\t", x0[i])); -- } -- Console.WriteLine(); -- if (finish) return; -- } forM ([0..n]) $ \i -> do forM ([0..n]) $ \j -> do printf "%14.10f\t" (test (a!!i!!j)) putStrLn "" putStrLn "" --print $ id' (8::Int) --print $ sum' ([1,2,3]::[Int]) --print $ append "abc" "def" --print $ log2 10 --print $ getHeight tree1 --return () import Text.Printf import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- ヤコビの反復法 row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double] row_loop row a b x0 x1 = let a1 = zip (a!!row) x0 a2 = take row a1 ++ drop (row + 1) a1 a3 = map (\(x, y) -> x * y) $ a2 s = sum a3 x = ((b!!row) - s) / a!!row!!row xs = x:x1 in if row >= (n - 1) then reverse xs else (row_loop (row+1) a b x0 xs) jacobi::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double]) jacobi a b x0 xs = let x1 = (row_loop 0 a b x0 []) cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1 in if cnt < 1 then (reverse xs, x1) else (jacobi a b x1 (x1:xs)) main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17::Double] let x0 = [0,0,0,0::Double] putStrLn "A" disp_matrix a putStrLn "B" disp_vector b putStrLn "" -- ヤコビの反復法 let (xs, x1) = jacobi a b x0 [] disp_matrix xs putStrLn "X" disp_vector x1 import Text.Printf import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- ガウス・ザイデル法 row_loop::Int->[[Double]]->[Double]->[Double]->[Double]->[Double] row_loop row a b x0 x1 = let a1 = zip (a!!row) x0 a2 = take row a1 ++ drop (row + 1) a1 a3 = map (\(x, y) -> x * y) $ a2 s = sum a3 x = ((b!!row) - s) / a!!row!!row xs = x:x1 in if row >= (n - 1) then reverse xs else (row_loop (row+1) a b ((take (row+1) (reverse xs)) ++ (drop (row+1) x0)) xs) gauss::[[Double]]->[Double]->[Double]->[[Double]]->([[Double]],[Double]) gauss a b x0 xs = let x1 = (row_loop 0 a b x0 []) cnt = length $ filter (>= 0.0000000001) $ map (\(x,y) -> abs(x - y)) $ zip x0 x1 in if cnt < 1 then (reverse xs, x1) else (gauss a b x1 (x1:xs)) main = do let a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6::Double]] let b = [20,16,8,17::Double] let x0 = [0,0,0,0::Double] putStrLn "A" disp_matrix a putStrLn "B" disp_vector b putStrLn "" -- ガウス・ザイデル法 let (xs, x1) = gauss a b x0 [] disp_matrix xs putStrLn "X" disp_vector x1 import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 各列で 一番値が大きい行を 探す get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double) get_max_row row col a max_row max_val = let -- 一番値が大きい行 (max_row2, max_val2) = if max_val < abs(a!!row!!col) then (row, abs(a!!row!!col)) else (max_row, max_val) in if row >= length(a) - 1 then (max_row2, max_val2) else get_max_row (row+1) col a max_row2 max_val2 -- ピボット選択 pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double]) pivoting pivot a b a2 b2 = let (max_row, max_val) = get_max_row 0 pivot a 0 0.0 a3 = (a!!max_row):a2 b3 = (b!!max_row):b2 a4 = (take (max_row) a) ++ (drop (max_row+1) a) b4 = (take (max_row) b) ++ (drop (max_row+1) b) in if pivot >= (n - 1) then (reverse a3, reverse b3) else pivoting (pivot+1) a4 b4 a3 b3 -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elim_loop pivot row a b = let s = a!!row!!pivot / a!!pivot!!pivot a_zip = zip (drop pivot (a!!pivot)) (drop pivot (a!!row)) a_map = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip new_row = (take pivot (a!!row)) ++ a_map a2 = (take row a ) ++ (new_row:(drop (row+1) a)) x = b!!row - b!!pivot * s b2 = (take row b) ++ (x:(drop (row+1) b)) in if row < (n - 1) then forward_elim_loop pivot (row+1) a2 b2 else (a2, b2) forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elimination pivot a b = let (a2, b2) = if pivot < (n - 1) then forward_elim_loop pivot (pivot+1) a b else (a, b) in if pivot < (n - 1) then forward_elimination (pivot+1) a2 b2 else (a2, b2) -- 後退代入 backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b) s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip x = (b!!row - s) / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]] let b = [8,17,20,16::Double] -- ピボット選択 let (a1, b1) = pivoting 0 a b [] [] putStrLn "Pivoting" putStrLn "A" disp_matrix a1 putStrLn "B" disp_vector b1 putStrLn "" -- 前進消去 let (a2, b2) = forward_elimination 0 a1 b1 putStrLn "Forward Elimination" putStrLn "A" disp_matrix a2 putStrLn "B" disp_vector b2 putStrLn "" -- 後退代入 let x = backward_substitution (n-1) a2 b2 putStrLn "X" disp_vector (reverse x) import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 各列で 一番値が大きい行を 探す get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double) get_max_row row col a max_row max_val = let -- 一番値が大きい行 (max_row2, max_val2) = if max_val < abs(a!!row!!col) then (row, abs(a!!row!!col)) else (max_row, max_val) in if row >= length(a) - 1 then (max_row2, max_val2) else get_max_row (row+1) col a max_row2 max_val2 -- ピボット選択 pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double]) pivoting pivot a b a2 b2 = let (max_row, max_val) = get_max_row 0 pivot a 0 0.0 a3 = (a!!max_row):a2 b3 = (b!!max_row):b2 a4 = (take (max_row) a) ++ (drop (max_row+1) a) b4 = (take (max_row) b) ++ (drop (max_row+1) b) in if pivot >= (n - 1) then (reverse a3, reverse b3) else pivoting (pivot+1) a4 b4 a3 b3 -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elim_loop pivot row a b = let s = a!!row!!pivot / a!!pivot!!pivot a_zip = zip (a!!pivot) (a!!row) new_row = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip a2 = (take row a) ++ (new_row:(drop (row+1) a)) x = b!!row - b!!pivot * s b2 = (take row b) ++ (x:(drop (row+1) b)) next_row = if row + 1 == pivot then row + 2 else row + 1 in if next_row < n then forward_elim_loop pivot next_row a2 b2 else (a2, b2) forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elimination pivot a b = let row = if pivot == 0 then 1 else 0 (a2, b2) = if pivot < n then forward_elim_loop pivot row a b else (a, b) in if pivot < (n - 1) then forward_elimination (pivot+1) a2 b2 else (a2, b2) -- 後退代入 backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let x = b!!row / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]] let b = [8,17,20,16::Double] -- ピボット選択 let (a1, b1) = pivoting 0 a b [] [] putStrLn "Pivoting" putStrLn "A" disp_matrix a1 putStrLn "B" disp_vector b1 putStrLn "" -- 前進消去 let (a2, b2) = forward_elimination 0 a1 b1 putStrLn "Forward Elimination" putStrLn "A" disp_matrix a2 putStrLn "B" disp_vector b2 putStrLn "" -- 後退代入 let x = backward_substitution (n-1) a2 b2 putStrLn "X" disp_vector (reverse x) import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 各列で 一番値が大きい行を 探す get_max_row::Int->Int->[[Double]]->Int->Double->(Int, Double) get_max_row row col a max_row max_val = let -- 一番値が大きい行 (max_row2, max_val2) = if max_val < abs(a!!row!!col) then (row, abs(a!!row!!col)) else (max_row, max_val) in if row >= length(a) - 1 then (max_row2, max_val2) else get_max_row (row+1) col a max_row2 max_val2 -- ピボット選択 pivoting::Int->[[Double]]->[Double]->[[Double]]->[Double]->([[Double]],[Double]) pivoting pivot a b a2 b2 = let (max_row, max_val) = get_max_row 0 pivot a 0 0.0 a3 = (a!!max_row):a2 b3 = (b!!max_row):b2 a4 = (take (max_row) a) ++ (drop (max_row+1) a) b4 = (take (max_row) b) ++ (drop (max_row+1) b) in if pivot >= (n - 1) then (reverse a3, reverse b3) else pivoting (pivot+1) a4 b4 a3 b3 -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elim_loop pivot row a b = let s = a!!row!!pivot / a!!pivot!!pivot a_zip = zip (drop (pivot+1) (a!!pivot)) (drop (pivot+1) (a!!row)) a_map = map (\(a_pivot, a_row) -> a_row - a_pivot * s) a_zip new_row = (take pivot (a!!row)) ++ (s:a_map ) a2 = (take row a ) ++ (new_row:(drop (row+1) a)) in if row < (n - 1) then forward_elim_loop pivot (row+1) a2 b else (a2, b) forward_elimination::Int->[[Double]]->[Double]->([[Double]],[Double]) forward_elimination pivot a b = let (a2, b2) = if pivot < (n - 1) then forward_elim_loop pivot (pivot+1) a b else (a, b) in if pivot < (n - 1) then forward_elimination (pivot+1) a2 b2 else (a2, b2) -- Ly=b から y を求める (前進代入) forward_substitution::Int->[[Double]]->[Double]->[Double] forward_substitution row a b = let a_zip = zip (take row (a!!row)) b s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip y = b!!row - s b2 = (take row b) ++ (y:(drop (row+1) b)) in if row < (n - 1) then y:(forward_substitution (row+1) a b2) else y:[] -- Ux=y から x を求める (後退代入) backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b) s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip x = (b!!row - s) / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[-1,-2,7,-2],[1,-1,-2,6],[9,2,1,1],[2,8,-2,1::Double]] let b = [8,17,20,16::Double] -- ピボット選択 let (a1, b1) = pivoting 0 a b [] [] putStrLn "Pivoting" putStrLn "A" disp_matrix a1 putStrLn "B" disp_vector b1 putStrLn "" -- 前進消去 let (a2, b2) = forward_elimination 0 a1 b1 putStrLn "LU" disp_matrix a2 -- Ly=b から y を求める (前進代入) let y = forward_substitution 0 a2 b2 putStrLn "Y" disp_vector (reverse y) -- Ux=y から x を求める (後退代入) let x = backward_substitution (n-1) a2 y putStrLn "X" disp_vector (reverse x) import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->Double->[[Double]] forward_elim_loop pivot row a a_sqrt = let a_zip = zip (take pivot (a!!row)) (take pivot (a!!pivot)) s = sum $ map(\(a_row, a_pivot) -> a_row * a_pivot) a_zip a_col = (a!!row!!pivot - s) / a_sqrt a_row = (take pivot (a!!row) ) ++ (a_col:(drop (pivot+1) (a!!row))) a_piv = (take row (a!!pivot)) ++ (a_col:(drop (row+1) (a!!pivot))) a1 = (take pivot a ) ++ (a_piv:(drop (pivot+1) a)) a2 = (take row a1) ++ (a_row:(drop (row+1) a)) in if row < (n - 1) then forward_elim_loop pivot (row+1) a2 a_sqrt else a2 forward_elimination::Int->[[Double]]->[[Double]] forward_elimination pivot a = let s = sum $ map(\a_col -> a_col * a_col) $ take pivot (a!!pivot) -- ここで根号の中が負の値になると計算できない! a_sqrt = (sqrt (a!!pivot!!pivot - s)) a2 = forward_elim_loop pivot pivot a a_sqrt in if pivot < (n - 1) then forward_elimination (pivot+1) a2 else a2 -- Ly=b から y を求める (前進代入) forward_substitution::Int->[[Double]]->[Double]->[Double] forward_substitution row a b = let a_zip = zip (take row (a!!row)) b s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip y = (b!!row - s) / a!!row!!row b2 = (take row b) ++ (y:(drop (row+1) b)) in if row < (n - 1) then y:(forward_substitution (row+1) a b2) else y:[] -- Ux=y から x を求める (後退代入) backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b) s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip x = (b!!row - s) / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20::Double]] let b = [34,68,96,125::Double] putStrLn "A" disp_matrix a putStrLn "B" disp_vector b putStrLn "" -- 前進消去 let a2 = forward_elimination 0 a putStrLn "LL^T" disp_matrix a2 -- Ly=b から y を求める (前進代入) let y = forward_substitution 0 a2 b putStrLn "Y" disp_vector y -- Ux=y から x を求める (後退代入) let x = backward_substitution (n-1) a2 y putStrLn "X" disp_vector (reverse x) import Text.Printf import Debug.Trace import Control.Monad n = 4::Int -- 2次元配列を表示 disp_matrix::[[Double]]->IO() disp_matrix matrix = do forM_ matrix $ \row -> do forM_ row $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 1次元配列を表示 disp_vector::[Double]->IO() disp_vector vector = do forM_ vector $ \elem -> do printf "%14.10f\t" elem putStrLn "" -- 前進消去 forward_elim_loop::Int->Int->[[Double]]->[[Double]] forward_elim_loop pivot col a = let a_zip = zip [0..col] $ zip (take col (a!!pivot)) (take col (a!!col)) s = sum $ map(\(k, (a_pivot, a_col)) -> a_pivot * a_col * (a!!k!!k)) a_zip a_col = (a!!pivot!!col - s) / (a!!col!!col) a_row = (take col (a!!pivot)) ++ (a_col:(drop (col+1) (a!!pivot))) a_piv = (take pivot (a!!col )) ++ (a_col:(drop (pivot+1) (a!!col))) a1 = (take col a ) ++ (a_piv:(drop (col+1) a)) a2 = (take pivot a1) ++ (a_row:(drop (pivot+1) a)) next_col = if col + 1 == pivot then pivot + 1 else col + 1 in if next_col < pivot then forward_elim_loop pivot next_col a2 else a2 forward_elimination::Int->[[Double]]->[[Double]] forward_elimination pivot a = let next_col = if pivot == 0 then 1 else 0 a1 = forward_elim_loop pivot next_col a s = sum $ map(\(k, a_col) -> a_col * a_col * a1!!k!!k) $ zip [0..pivot] $ take pivot (a1!!pivot) a_col = a1!!pivot!!pivot - s a_row = (take pivot (a1!!pivot)) ++ (a_col:(drop (pivot+1) (a1!!pivot))) a2 = (take pivot a1 ) ++ (a_row:(drop (pivot+1) a1 )) in if pivot < (n - 1) then forward_elimination (pivot+1) a2 else a2 -- Ly=b から y を求める (前進代入) forward_substitution::Int->[[Double]]->[Double]->[Double] forward_substitution row a b = let a_zip = zip (take row (a!!row)) b s = sum $ map(\(a_col, b_col) -> a_col * b_col) a_zip y = (b!!row - s) b2 = (take row b) ++ (y:(drop (row+1) b)) in if row < (n - 1) then y:(forward_substitution (row+1) a b2) else y:[] -- Ux=y から x を求める (後退代入) backward_substitution::Int->[[Double]]->[Double]->[Double] backward_substitution row a b = let a_zip = zip (drop (row+1) (a!!row)) (drop (row+1) b) s = sum $ map(\(a_col, b_col) -> a_col * a!!row!!row * b_col) a_zip x = (b!!row - s) / (a!!row!!row) b2 = (take row b) ++ (x:(drop (row+1) b)) in if row > 0 then x:(backward_substitution (row-1) a b2) else x:[] main = do let a = [[5,2,3,4],[2,10,6,7],[3,6,15,9],[4,7,9,20::Double]] let b = [34,68,96,125::Double] putStrLn "A" disp_matrix a putStrLn "B" disp_vector b putStrLn "" -- 前進消去 let a2 = forward_elimination 0 a putStrLn "LDL^T" disp_matrix a2 -- Ly=b から y を求める (前進代入) let y = forward_substitution 0 a2 b putStrLn "Y" disp_vector y -- Ux=y から x を求める (後退代入) let x = backward_substitution (n-1) a2 y putStrLn "X" disp_vector (reverse x) |