Análisis

Método de bisección para aproximar raíces de funciones

PorJosé A. Alonso 24-junio-202221-junio-2022

El método de bisección para calcular un cero de una función en el intervalo [a,b] se basa en el teorema de Bolzano:

«Si f(x) es una función continua en el intervalo [a, b], y si, además, en los extremos del intervalo la función f(x) toma valores de signo opuesto (f(a) * f(b) < 0), entonces existe al menos un valor c en (a, b) para el que f(c) = 0".

El método para calcular un cero de la función f en el intervalo [a,b] con un error menor que e consiste en tomar el punto medio del intervalo c = (a+b)/2 y considerar los siguientes casos:

Si |f(c)| < e, hemos encontrado una aproximación del punto que anula f en el intervalo con un error aceptable.
Si f(c) tiene signo distinto de f(a), repetir el proceso en el intervalo [a,c].
Si no, repetir el proceso en el intervalo [c,b].

Definir la función

   biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

1	biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double

tal que (biseccion f a b e) es una aproximación del punto del intervalo [a,b] en el que se anula la función f, con un error menor que e, calculada mediante el método de la bisección. Por ejemplo,

   biseccion (\x -> x^2 - 3) 0 5 0.01             ==  1.7333984375
   biseccion (\x -> x^3 - x - 2) 0 4 0.01         ==  1.521484375
   biseccion cos 0 2 0.01                         ==  1.5625
   biseccion (\x -> log (50-x) - 4) (-10) 3 0.01  ==  -5.125

biseccion (\x -> x^2 - 3) 0 5 0.01 == 1.7333984375

biseccion (\x -> x^3 - x - 2) 0 4 0.01 == 1.521484375

biseccion cos 0 2 0.01 == 1.5625

biseccion (\x -> log (50-x) - 4) (-10) 3 0.01 == -5.125

Soluciones

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

-- 1ª solución
-- ===========

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion1 f a b e  
    | abs (f c) < e = c
    | f a * f c < 0 = biseccion1 f a c e
    | otherwise     = biseccion1 f c b e
    where c = (a+b)/2

-- 2ª solución
-- ===========

biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion2 f a b e = aux a b
  where aux a' b' | abs (f c) < e = c
                  | f a' * f c < 0 = aux a' c 
                  | otherwise     = aux c b'
          where c = (a'+b')/2

-- Comprobación de equivalencia
-- ============================

newtype Polinomio = P [Int]
  deriving Show

valorPolinomio :: Polinomio -> Double -> Double
valorPolinomio (P cs) x =
  sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio
polinomioArbitrario 0 = return (P [])
polinomioArbitrario n = do
  c <- choose (-10,10)
  (P xs) <- polinomioArbitrario (n `div` 2)
  return (P (c:xs))

instance Arbitrary Polinomio where
  arbitrary = sized polinomioArbitrario

-- La propiedad es
prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property
prop_biseccion p a b e =
  f a * f b < 0 && e > 0 ==>
  biseccion1 f a b e =~ biseccion2 f a b e
  where
    f = valorPolinomio p 
    x =~ y = abs (x - y) < 0.001

-- La comprobación es
--    λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion
--    +++ OK, passed 100 tests; 2156 discarded.

import Test.QuickCheck (Arbitrary, Gen, Property, (==>), arbitrary, choose, sized, quickCheck)

-- 1ª solución

-- ===========

biseccion1 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion1 f a b e

| abs (f c) < e = c

| f a * f c < 0 = biseccion1 f a c e

| otherwise = biseccion1 f c b e

where c = (a+b)/2

-- 2ª solución

-- ===========

biseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double

biseccion2 f a b e = aux a b

where aux a' b' | abs (f c) < e = c

| f a' * f c < 0 = aux a' c

| otherwise = aux c b'

where c = (a'+b')/2

-- Comprobación de equivalencia

-- ============================

newtype Polinomio = P [Int]

deriving Show

valorPolinomio :: Polinomio -> Double -> Double

valorPolinomio (P cs) x =

sum [fromIntegral c * x^n | (c,n) <- zip cs [0..]]

polinomioArbitrario :: Int -> Gen Polinomio

polinomioArbitrario 0 = return (P [])

polinomioArbitrario n = do

c <- choose (-10,10)

(P xs) <- polinomioArbitrario (n `div` 2)

return (P (c:xs))

instance Arbitrary Polinomio where

arbitrary = sized polinomioArbitrario

-- La propiedad es

prop_biseccion :: Polinomio -> Double -> Double -> Double -> Property

prop_biseccion p a b e =

f a * f b < 0 && e > 0 ==>

biseccion1 f a b e =~ biseccion2 f a b e

where

f = valorPolinomio p

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheckWith (stdArgs {maxDiscardRatio=30}) prop_biseccion

-- +++ OK, passed 100 tests; 2156 discarded.

El código se encuentra en GitHub.

Cálculo de la suma 11! + 22! + 33! + … + nn!

PorJosé A. Alonso 22-junio-202221-junio-2022

Definir la función

   suma :: Integer -> Integer

1	suma :: Integer -> Integer

tal que (suma n) es la suma 1·1! + 2·2! + 3·3! + ... + n·n!. Por ejemplo,

   suma 1  ==  1
   suma 2  ==  5
   suma 3  ==  23
   suma 4  ==  119
   suma 5  ==  719
   take 9 (show (suma 70000))  ==  "823780458"

suma 1 == 1

suma 2 == 5

suma 3 == 23

suma 4 == 119

suma 5 == 719

take 9 (show (suma 70000)) == "823780458"

Soluciones

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución
-- ===========

suma1 :: Integer -> Integer
suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer
factorial n = product [1..n]

-- 2ª solución
-- ===========

suma2 :: Integer -> Integer
suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]
factoriales = scanl (*) 1 [2..]

-- 3ª solución
-- ===========

-- Basada en los siguientes cálculos
--    λ> [suma1 n | n <- [0..10]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]
--    λ> [factorial n | n <- [0..10]]
--    [1,1,2,6,24,120,720,5040,40320,362880,3628800]
--    λ> [factorial n | n <- [1..11]]
--    [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]
--    λ> [factorial n - 1 | n <- [1..11]]
--    [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer
suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_suma :: Positive Integer -> Bool
prop_suma (Positive n) =
  all (== suma1 n)
      [suma2 n,
       suma3 n]
  
-- La comprobación es
--    λ> quickCheck prop_suma
--    +++ OK, passed 100 tests.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> take 5 (show (suma1 4000))
--    "73170"
--    (5.04 secs, 16,225,195,448 bytes)
--    λ> take 5 (show (suma2 4000))
--    "73170"
--    (0.08 secs, 35,862,152 bytes)
--    λ> take 5 (show (suma3 4000))
--    "73170"
--    (0.01 secs, 12,896,968 bytes)
--    
--    
--    λ> take 5 (show (suma2 40000))
--    "83669"
--    (1.82 secs, 4,549,612,264 bytes)
--    λ> take 5 (show (suma3 40000))
--    "83669"
--    (0.24 secs, 1,620,976,984 bytes)

import Test.QuickCheck (Positive (Positive), quickCheck)

-- 1ª solución

-- ===========

suma1 :: Integer -> Integer

suma1 n = sum [k * factorial k | k <- [1..n]]

factorial :: Integer -> Integer

factorial n = product [1..n]

-- 2ª solución

-- ===========

suma2 :: Integer -> Integer

suma2 n = sum (zipWith (*) [1..n] factoriales)

factoriales :: [Integer]

factoriales = scanl (*) 1 [2..]

-- 3ª solución

-- ===========

-- Basada en los siguientes cálculos

-- λ> [suma1 n | n <- [0..10]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

-- λ> [factorial n | n <- [0..10]]

-- [1,1,2,6,24,120,720,5040,40320,362880,3628800]

-- λ> [factorial n | n <- [1..11]]

-- [1,2,6,24,120,720,5040,40320,362880,3628800,39916800]

-- λ> [factorial n - 1 | n <- [1..11]]

-- [0,1,5,23,119,719,5039,40319,362879,3628799,39916799]

suma3 :: Integer -> Integer

suma3 n = factorial (n+1) - 1

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_suma :: Positive Integer -> Bool

prop_suma (Positive n) =

all (== suma1 n)

[suma2 n,

suma3 n]

-- La comprobación es

-- λ> quickCheck prop_suma

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> take 5 (show (suma1 4000))

-- "73170"

-- (5.04 secs, 16,225,195,448 bytes)

-- λ> take 5 (show (suma2 4000))

-- "73170"

-- (0.08 secs, 35,862,152 bytes)

-- λ> take 5 (show (suma3 4000))

-- "73170"

-- (0.01 secs, 12,896,968 bytes)

-- λ> take 5 (show (suma2 40000))

-- "83669"

-- (1.82 secs, 4,549,612,264 bytes)

-- λ> take 5 (show (suma3 40000))

-- "83669"

-- (0.24 secs, 1,620,976,984 bytes)

El código se encuentra en GitHub.

Cálculo aproximado de integrales definidas

PorJosé A. Alonso 21-junio-202221-junio-2022

La integral definida de una función f entre los límites a y b puede calcularse mediante la regla del rectángulo usando la fórmula

   h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

1	h * (f(a+h/2) + f(a+h+h/2) + f(a+2h+h/2) + ... + f(a+n*h+h/2))

con a+n*h+h/2 <= b < a+(n+1)*h+h/2 y usando valores pequeños para h.

Definir la función

   integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

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

   integral 0 1 (^3) 0.01  ==  0.24998750000000042

1	integral 0 1 (^3) 0.01 == 0.24998750000000042

Otros ejemplos son

   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

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

import Test.QuickCheck.HigherOrder (quickCheck')
import Test.QuickCheck (Property, (==>), quickCheck)

-- 1ª solución
-- ===========

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral1 a b f h
  | a+h/2 > b = 0
  | otherwise = h * f (a+h/2) + integral1 (a+h) b f h

-- 2ª solución
-- ===========

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral2 a b f h = aux a where
  aux x | x+h/2 > b = 0
        | otherwise = h * f (x+h/2) + aux (x+h)

-- 3ª solución
-- ===========

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a
integral3 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)

-- Comprobación de equivalencia
-- ============================

-- La propiedad es
prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property
prop_integral a b f h =
  a < b && h > 0 ==> 
  all (=~ integral1 a' b' f' h')
      [integral2 a' b' f' h',
       integral3 a' b' f' h']
  where
    a' = fromIntegral a
    b' = fromIntegral b
    h' = fromIntegral h
    f' = fromIntegral . f. round
    x =~ y = abs (x - y) < 0.001
        
-- La comprobación es
--    λ> quickCheck' prop_integral
--    +++ OK, passed 100 tests; 385 discarded.

-- Comparación de eficiencia
-- =========================

-- La comparación es
--    λ> integral1 0 10 (^3) 0.00001
--    2499.999999881125
--    (2.63 secs, 1,491,006,744 bytes)
--    λ> integral2 0 10 (^3) 0.00001
--    2499.999999881125
--    (1.93 secs, 1,419,006,696 bytes)
--    λ> integral3 0 10 (^3) 0.00001
--    2499.9999998811422
--    (1.28 secs, 817,772,216 bytes)

import Test.QuickCheck.HigherOrder (quickCheck')

import Test.QuickCheck (Property, (==>), quickCheck)

-- 1ª solución

-- ===========

integral1 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral1 a b f h

| a+h/2 > b = 0

| otherwise = h * f (a+h/2) + integral1 (a+h) b f h

-- 2ª solución

-- ===========

integral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral2 a b f h = aux a where

aux x | x+h/2 > b = 0

| otherwise = h * f (x+h/2) + aux (x+h)

-- 3ª solución

-- ===========

integral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a

integral3 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)

-- Comprobación de equivalencia

-- ============================

-- La propiedad es

prop_integral :: Int -> Int -> (Int -> Int) -> Int -> Property

prop_integral a b f h =

a < b && h > 0 ==>

all (=~ integral1 a' b' f' h')

[integral2 a' b' f' h',

integral3 a' b' f' h']

where

a' = fromIntegral a

b' = fromIntegral b

h' = fromIntegral h

f' = fromIntegral . f. round

x =~ y = abs (x - y) < 0.001

-- La comprobación es

-- λ> quickCheck' prop_integral

-- +++ OK, passed 100 tests; 385 discarded.

-- Comparación de eficiencia

-- =========================

-- La comparación es

-- λ> integral1 0 10 (^3) 0.00001

-- 2499.999999881125

-- (2.63 secs, 1,491,006,744 bytes)

-- λ> integral2 0 10 (^3) 0.00001

-- 2499.999999881125

-- (1.93 secs, 1,419,006,696 bytes)

-- λ> integral3 0 10 (^3) 0.00001

-- 2499.9999998811422

-- (1.28 secs, 817,772,216 bytes)

El código se encuentra en GitHub.

Método de bisección para aproximar raíces de funciones

Soluciones

Cálculo de la suma 11! + 22! + 33! + … + nn!

Soluciones

Cálculo aproximado de integrales definidas

Soluciones

Entradas recientes

Buscar

Correo electrónico