Integración por el método de los rectángulos
La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo (ver en http://bit.ly/1FDhZ1z) usando la fórmula
\[ h(f(a+\frac{h}{2}) + f(a+h+\frac{h}{2}) + f(a+2h+\frac{h}{2}) + … + f(a+nh+\frac{h}{2}))\]
con \(a+nh+\dfrac{h}{2} \leq b < a+(n+1)h+\dfrac{h}{2}\) y usando valores pequeños para \(h\).
Definir la función
1 |
integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a |
tal que integral a b f h
es el valor de dicha expresión. Por ejemplo, el cálculo de la integral de (f(x) = x^3) entre 0 y 1, con paso 0.01, es
1 |
integral 0 1 (^3) 0.01 == 0.24998750000000042 |
Otros ejemplos son
1 2 3 4 |
integral 0 1 (^4) 0.01 == 0.19998333362500048 integral 0 1 (\x -> 3*x^2 + 4*x^3) 0.01 == 1.9999250000000026 log 2 - integral 1 2 (\x -> 1/x) 0.01 == 3.124931644782336e-6 pi - 4 * integral 0 1 (\x -> 1/(x^2+1)) 0.01 == -8.333333331389525e-6 |
Soluciones
A continuación se muestran las soluciones en Haskell y las soluciones en Python.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 |
module Integracion_por_rectangulos where import Test.Hspec (Spec, hspec, it, shouldBe, shouldSatisfy) -- 1ª solución -- =========== integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral a b f h = h * suma (a+h/2) b (+h) f -- (suma a b s f) es l valor de -- f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...)) -- hasta que s(s(...(s(a))...)) > b. Por ejemplo, -- suma 2 5 (1+) (^3) == 224 suma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a suma a b s f = sum [f x | x <- sucesion a b s] -- (sucesion x y s) es la lista -- [a, s(a), s(s(a), ..., s(...(s(a))...)] -- hasta que s(s(...(s(a))...)) > b. Por ejemplo, -- sucesion 3 20 (+2) == [3,5,7,9,11,13,15,17,19] sucesion :: Ord a => a -> a -> (a -> a) -> [a] sucesion a b s = takeWhile (<=b) (iterate s a) -- 2ª solución -- =========== integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral2 a b f h | a+h/2 > b = 0 | otherwise = h * f (a+h/2) + integral2 (a+h) b f h -- 3ª solución -- =========== integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a integral3 a b f h = aux a where aux x | x+h/2 > b = 0 | otherwise = h * f (x+h/2) + aux (x+h) -- Comparación de eficiencia -- λ> integral 0 10 (^3) 0.00001 -- 2499.9999998811422 -- (4.62 secs, 1084774336 bytes) -- λ> integral2 0 10 (^3) 0.00001 -- 2499.999999881125 -- (7.90 secs, 1833360768 bytes) -- λ> integral3 0 10 (^3) 0.00001 -- 2499.999999881125 -- (7.27 secs, 1686056080 bytes) -- Verificación -- ============ verifica :: IO () verifica = hspec spec spec :: Spec spec = do it "e1" $ integral' 0 1 (^(3::Int)) 0.01 `shouldBe` 0.24998750000000042 it "e2" $ integral' 0 1 (^(4::Int)) 0.01 `shouldBe` 0.19998333362500048 it "e3" $ integral' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldBe` 1.9999250000000026 it "e4" $ log 2 - integral' 1 2 (1 /) 0.01 `shouldBe` 3.124931644782336e-6 it "e5" $ pi - 4 * integral' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldBe` -8.333333331389525e-6 it "e1b" $ integral2' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042) it "e2b" $ integral2' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048) it "e3b" $ integral2' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026) it "e4b" $ log 2 - integral2' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6) it "e5b" $ pi - 4 * integral2' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6)) it "e1c" $ integral3' 0 1 (^(3::Int)) 0.01 `shouldSatisfy` (~= 0.24998750000000042) it "e2c" $ integral3' 0 1 (^(4::Int)) 0.01 `shouldSatisfy` (~= 0.19998333362500048) it "e3c" $ integral3' 0 1 (\x -> 3*x^(2::Int) + 4*x^(3::Int)) 0.01 `shouldSatisfy` (~= 1.9999250000000026) it "e4c" $ log 2 - integral3' 1 2 (1 /) 0.01 `shouldSatisfy` (~= 3.124931644782336e-6) it "e5c" $ pi - 4 * integral3' 0 1 (\x -> 1/(x^(2::Int)+1)) 0.01 `shouldSatisfy` (~= (-8.333333331389525e-6)) where integral', integral2', integral3' :: Double -> Double -> (Double -> Double) -> Double -> Double integral' = integral integral2' = integral2 integral3' = integral3 a ~= b = abs (a - b) < 0.00001 -- La verificación es -- λ> verifica -- -- Finished in 0.0058 seconds -- 15 examples, 0 failures |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 |
from math import log, pi from typing import Callable # 1ª solución # =========== # sucesion(x, y, s) es la lista # [a, s(a), s(s(a), ..., s(...(s(a))...)] # hasta que s(s(...(s(a))...)) > b. Por ejemplo, # sucesion(3, 20, lambda x : x+2) == [3,5,7,9,11,13,15,17,19] def sucesion(a: float, b: float, s: Callable[[float], float]) -> list[float]: xs = [] while a <= b: xs.append(a) a = s(a) return xs # suma(a, b, s, f) es el valor de # f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...)) # hasta que s(s(...(s(a))...)) > b. Por ejemplo, # suma(2, 5, lambda x: x+1, lambda x: x**3) == 224 def suma(a: float, b: float, s: Callable[[float], float], f: Callable[[float], float]) -> float: return sum(f(x) for x in sucesion(a, b, s)) def integral(a: float, b: float, f: Callable[[float], float], h: float) -> float: return h * suma(a+h/2, b, lambda x: x+h, f) # 2ª solución # =========== def integral2(a: float, b: float, f: Callable[[float], float], h: float) -> float: if a+h/2 > b: return 0 return h * f(a+h/2) + integral2(a+h, b, f, h) # 3ª solución # =========== def integral3(a: float, b: float, f: Callable[[float], float], h: float) -> float: def aux(x: float) -> float: if x+h/2 > b: return 0 return h * f(x+h/2) + aux(x+h) return aux(a) # Verificación # ============ def test_integral() -> None: def aproximado(a: float, b: float) -> bool: return abs(a - b) < 0.00001 assert integral(0, 1, lambda x : x**3, 0.01) == 0.24998750000000042 assert integral(0, 1, lambda x : x**4, 0.01) == 0.19998333362500054 assert integral(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01) == 1.9999250000000026 assert log(2) - integral(1, 2, lambda x : 1/x, 0.01) == 3.124931644782336e-6 assert pi - 4 * integral(0, 1, lambda x : 1/(x**2+1), 0.01) == -8.333333331389525e-6 assert aproximado(integral2(0, 1, lambda x : x**3, 0.01), 0.24998750000000042) assert aproximado(integral2(0, 1, lambda x : x**4, 0.01), 0.19998333362500054) assert aproximado(integral2(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01), 1.9999250000000026) assert aproximado(log(2) - integral2(1, 2, lambda x : 1/x, 0.01), 3.124931644782336e-6) assert aproximado(pi - 4 * integral2(0, 1, lambda x : 1/(x**2+1), 0.01), -8.333333331389525e-6) assert aproximado(integral3(0, 1, lambda x : x**3, 0.01), 0.24998750000000042) assert aproximado(integral3(0, 1, lambda x : x**4, 0.01), 0.19998333362500054) assert aproximado(integral3(0, 1, lambda x : 3*x**2 + 4*x**3, 0.01), 1.9999250000000026) assert aproximado(log(2) - integral3(1, 2, lambda x : 1/x, 0.01), 3.124931644782336e-6) assert aproximado(pi - 4 * integral3(0, 1, lambda x : 1/(x**2+1), 0.01), -8.333333331389525e-6) print("Verificado") # La verificación es # >>> test_integral() # Verificado |