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
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!!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 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] -- ヤコビの反復法 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 x2 = (take (row+1) (reverse xs)) ++ (drop (row+1) x0) in if row >= (n - 1) then reverse xs else (row_loop (row+1) a b x2 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] -- ガウス・ザイデル法 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) -- 前進代入 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:[] -- 後退代入 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 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(\col -> col * 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 -- 前進代入 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:[] -- 後退代入 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 -- L^Tx=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, col) -> col * 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 -- 前進代入 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:[] -- 後退代入 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 -- DL^Tx=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 "" -- 各列で 一番値が大きい行を 探す 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) -- 前進代入 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:[] -- 後退代入 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 y -- Ux=y から x を求める (後退代入) let x = backward_substitution (n-1) a2 y putStrLn "X" disp_vector (reverse x)