Numerical integration: Difference between revisions
Content deleted Content added
SqrtNegInf (talk | contribs) m →{{header|Perl 6}}: cosmetic fix |
→{{header|Haskell}}: Added section labels to output, pruned some redundancy in output code |
||
Line 2,276: | Line 2,276: | ||
The whole program: |
The whole program: |
||
<lang haskell>approx |
<lang haskell>approx |
||
:: Fractional a |
|||
=> (a1 -> a) -> [a1] -> [a] -> a |
|||
approx f xs ws = |
|||
sum |
|||
[ w * f x |
|||
| (x, w) <- zip xs ws ] |
|||
integrateOpen |
|||
integrateOpen :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a |
|||
:: Fractional a |
|||
integrateOpen v vs f a b n = approx f xs ws * h / v where |
|||
=> a -> [a] -> (a -> a) -> a -> a -> Int -> a |
|||
m = fromIntegral (length vs) * n |
|||
integrateOpen v vs f a b n = approx f xs ws * h / v |
|||
h = (b-a) / fromIntegral m |
|||
where |
|||
ws = concat $ replicate n vs |
|||
m = fromIntegral (length vs) * n |
|||
c = a + h/2 |
|||
h = (b - a) / fromIntegral m |
|||
ws = concat $ replicate n vs |
|||
c = a + h / 2 |
|||
xs = |
|||
[ c + h * fromIntegral i |
|||
| i <- [0 .. m - 1] ] |
|||
integrateClosed |
|||
integrateClosed :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a |
|||
:: Fractional a |
|||
integrateClosed v vs f a b n = approx f xs ws * h / v where |
|||
=> a -> [a] -> (a -> a) -> a -> a -> Int -> a |
|||
m = fromIntegral (length vs - 1) * n |
|||
integrateClosed v vs f a b n = approx f xs ws * h / v |
|||
h = (b-a) / fromIntegral m |
|||
where |
|||
ws = overlap n vs |
|||
m = fromIntegral (length vs - 1) * n |
|||
xs = [a + h * fromIntegral i | i <- [0..m]] |
|||
h = (b - a) / fromIntegral m |
|||
ws = overlap n vs |
|||
xs = |
|||
[ a + h * fromIntegral i |
|||
| i <- [0 .. m] ] |
|||
overlap |
|||
overlap :: Num a => Int -> [a] -> [a] |
|||
:: Num a |
|||
overlap n [] = [] |
|||
=> Int -> [a] -> [a] |
|||
overlap n (x:xs) = x : inter n xs where |
|||
overlap n [] = [] |
|||
inter 1 ys = ys |
|||
overlap n (x:xs) = x : inter n xs |
|||
where |
|||
inter n [y] = (x+y) : inter (n-1) xs |
|||
inter |
inter 1 ys = ys |
||
inter n [] = x : inter (n - 1) xs |
|||
inter n [y] = (x + y) : inter (n - 1) xs |
|||
intLeftRect = integrateClosed 1 [1,0] |
|||
inter n (y:ys) = y : inter n ys |
|||
intMidRect = integrateOpen 1 [1] |
|||
intRightRect = integrateClosed 1 [0,1] |
|||
intTrapezium = integrateClosed 2 [1,1] |
|||
intSimpson = integrateClosed 3 [1,4,1] |
|||
uncurry4 :: (t1 -> t2 -> t3 -> t4 -> t) -> (t1, t2, t3, t4) -> t |
|||
uncurry4 f ~(a, b, c, d) = f a b c d |
uncurry4 f ~(a, b, c, d) = f a b c d |
||
-- TEST ---------------------------------------------------------------------- |
|||
main = do |
|||
ms |
|||
let m1 = "rectangular left: " |
|||
:: Fractional a |
|||
let m2 = "rectangular middle: " |
|||
=> [(String, (a -> a) -> a -> a -> Int -> a)] |
|||
let m3 = "rectangular right: " |
|||
ms = |
|||
let m4 = "trapezium: " |
|||
[ ("rectangular left", integrateClosed 1 [1, 0]) |
|||
let m5 = "simpson: " |
|||
, ("rectangular middle", integrateOpen 1 [1]) |
|||
, ("rectangular right", integrateClosed 1 [0, 1]) |
|||
let arg1 = ((\x -> x ^ 3), 0, 1, 100) |
|||
, ("trapezium", integrateClosed 2 [1, 1]) |
|||
putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg1) |
|||
, ("simpson", integrateClosed 3 [1, 4, 1]) |
|||
putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg1) |
|||
] |
|||
putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg1) |
|||
putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg1) |
|||
putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg1) |
|||
putStrLn "" |
|||
let arg2 = ((\x -> 1 / x), 1, 100, 1000) |
|||
putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg2) |
|||
putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg2) |
|||
putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg2) |
|||
putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg2) |
|||
putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg2) |
|||
putStrLn "" |
|||
let arg3 = ((\x -> x), 0, 5000, 5000000) |
|||
putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg3) |
|||
putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg3) |
|||
putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg3) |
|||
putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg3) |
|||
putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg3) |
|||
putStrLn "" |
|||
expressions |
|||
let arg4 = ((\x -> x), 0, 6000, 6000000) |
|||
:: (Fractional a, Num t, Num t1, Num t2) |
|||
putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg4) |
|||
=> [(String, (a -> a, t, t1, t2))] |
|||
putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg4) |
|||
expressions = |
|||
putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg4) |
|||
[ ("x^3", ((^ 3), 0, 1, 100)) |
|||
putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg4) |
|||
, ("1/x", ((1 /), 1, 100, 1000)) |
|||
putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg4)</lang> |
|||
, ("x", (id, 0, 5000, 500000)) |
|||
, ("x", (id, 0, 6000, 600000)) |
|||
] |
|||
main :: IO () |
|||
Output: |
|||
main = |
|||
<pre>rectangular left: 0.24502500000000005 |
|||
mapM_ |
|||
rectangular middle: 0.24998750000000006 |
|||
(\(s, e@(f, a, b, n)) -> do |
|||
rectangular right: 0.25502500000000006 |
|||
putStrLn $ |
|||
trapezium: 0.25002500000000005 |
|||
justifyLeft 20 ' ' ("f(x) = " ++ s) ++ |
|||
unwords (show <$> [a, b]) ++ " " ++ show n |
|||
mapM_ |
|||
(\(s, integration) -> |
|||
putStrLn |
|||
(justifyLeft 20 ' ' (s ++ ":") ++ show (uncurry4 integration e))) |
|||
ms |
|||
putStrLn []) |
|||
expressions |
|||
where |
|||
justifyLeft n c s = take n (s ++ replicate n c)</lang> |
|||
{{Out}} |
|||
<pre>f(x) = x^3 0.0 1.0 100 |
|||
rectangular left: 0.24502500000000005 |
|||
rectangular middle: 0.24998750000000006 |
|||
rectangular right: 0.25502500000000006 |
|||
trapezium: 0.25002500000000005 |
|||
simpson: 0.25000000000000006 |
|||
f(x) = 1/x 1.0 100.0 1000 |
|||
rectangular left: 4.65499105751468 |
|||
rectangular |
rectangular left: 4.65499105751468 |
||
rectangular |
rectangular middle: 4.604762548678376 |
||
rectangular right: 4.55698105751468 |
|||
trapezium: 4.605986057514681 |
|||
trapezium: 4.605986057514681 |
|||
simpson: 4.605170384957135 |
|||
f(x) = x 0.0 5000.0 500000 |
|||
rectangular left: 1.24999975e7 |
|||
rectangular |
rectangular left: 1.2499975000000006e7 |
||
rectangular |
rectangular middle: 1.2499999999999993e7 |
||
rectangular right: 1.2500025000000006e7 |
|||
trapezium: 1.25e7 |
|||
trapezium: 1.2500000000000006e7 |
|||
simpson: 1.2499999999999998e7 |
|||
f(x) = x 0.0 6000.0 600000 |
|||
rectangular left: 1.7999997000000004e7 |
|||
rectangular |
rectangular left: 1.7999970000000004e7 |
||
rectangular |
rectangular middle: 1.7999999999999993e7 |
||
rectangular right: 1.8000030000000004e7 |
|||
trapezium: 1.8000000000000004e7 |
|||
trapezium: 1.8000000000000004e7 |
|||
simpson: 1.7999999999999996e7</pre> |
|||
Runtime: about 7 seconds. |
Runtime: about 7 seconds. |
||