Cálculo numérico

Aproximación del número pi

PorJosé A. Alonso 30-junio-20222-julio-2022

Una forma de aproximar el número π es usando la siguiente igualdad:

    π         1     1·2     1·2·3     1·2·3·4     
   --- = 1 + --- + ----- + ------- + --------- + ....
    2         3     3·5     3·5·7     3·5·7·9

π 1 1·2 1·2·3 1·2·3·4

--- = 1 + --- + ----- + ------- + --------- + ....

2 3 3·5 3·5·7 3·5·7·9

Es decir, la serie cuyo término general n-ésimo es el cociente entre el producto de los primeros n números y los primeros n números impares:

               Π i   
   s(n) =  -----------
            Π (2*i+1)

Π i

s(n) = -----------

Π (2*i+1)

Definir la función

   aproximaPi :: Integer -> Double

1	aproximaPi :: Integer -> Double

tal que (aproximaPi n) es la aproximación del número π calculada con la serie anterior hasta el término n-ésimo. Por ejemplo,

   aproximaPi 10     == 3.1411060206
   aproximaPi 20     == 3.1415922987403397
   aproximaPi 30     == 3.1415926533011596
   aproximaPi 40     == 3.1415926535895466
   aproximaPi 50     == 3.141592653589793
   aproximaPi (10^4) == 3.141592653589793
   pi                == 3.141592653589793

aproximaPi 10 == 3.1411060206

aproximaPi 20 == 3.1415922987403397

aproximaPi 30 == 3.1415926533011596

aproximaPi 40 == 3.1415926535895466

aproximaPi 50 == 3.141592653589793

aproximaPi (10^4) == 3.141592653589793

pi == 3.141592653589793

Soluciones

import Data.Ratio ((%))
import Data.List (genericTake)
import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

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

aproximaPi1 :: Integer -> Double
aproximaPi1 n = 
  fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

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

aproximaPi2 :: Integer -> Double
aproximaPi2 0 = 2
aproximaPi2 n = 
  aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

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

aproximaPi3 :: Integer -> Double
aproximaPi3 n = 
  fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,
--    λ> take 10 numeradores
--    [1,2,6,24,120,720,5040,40320,362880,3628800]
numeradores :: [Integer]
numeradores = scanl (*) 1 [2..]

-- denominadores es la sucesión de los denominadores. Por ejemplo,
--    λ> take 10 denominadores
--    [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]
denominadores :: [Integer]
denominadores = scanl (*) 3 [5, 7..]
 
-- 4ª solución
-- ===========

aproximaPi4 :: Integer -> Double
aproximaPi4 n = 
  read (x : "." ++ xs)
  where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer
aproximaPi4' n = 
  2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores))) 
  where p = 10^n

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

-- La propiedad es
prop_aproximaPi :: Property
prop_aproximaPi =
  forAll (arbitrary `suchThat` (> 3)) $ \n ->
  all (=~ aproximaPi1 n)
      [aproximaPi2 n,
       aproximaPi3 n,
       aproximaPi4 n]

(=~) :: Double -> Double -> Bool
x =~ y = abs (x - y) < 0.001

-- La comprobación es
--    λ> quickCheck prop_aproximaPi
--    +++ OK, passed 100 tests.

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

-- La comparación es
--    λ> aproximaPi1 3000 
--    3.141592653589793
--    (4.96 secs, 27,681,824,408 bytes)
--    λ> aproximaPi2 3000 
--    3.1415926535897922
--    (3.00 secs, 20,496,194,496 bytes)
--    λ> aproximaPi3 3000 
--    3.141592653589793
--    (3.13 secs, 13,439,528,432 bytes)
--    λ> aproximaPi4 3000 
--    3.141592653589793
--    (0.09 secs, 23,142,144 bytes)

import Data.Ratio ((%))

import Data.List (genericTake)

import Test.QuickCheck (Property, arbitrary, forAll, suchThat, quickCheck)

-- 1ª solución

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

aproximaPi1 :: Integer -> Double

aproximaPi1 n =

fromRational (2 * sum [product [1..i] % product [1,3..2*i+1] | i <- [0..n]])

-- 2ª solución

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

aproximaPi2 :: Integer -> Double

aproximaPi2 0 = 2

aproximaPi2 n =

aproximaPi2 (n-1) + fromRational (2 * product [1..n] % product [3,5..2*n+1])

-- 3ª solución

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

aproximaPi3 :: Integer -> Double

aproximaPi3 n =

fromRational (2 * (1 + sum (zipWith (%) numeradores (genericTake n denominadores))))

-- numeradores es la sucesión de los numeradores. Por ejemplo,

-- λ> take 10 numeradores

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

numeradores :: [Integer]

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

-- denominadores es la sucesión de los denominadores. Por ejemplo,

-- λ> take 10 denominadores

-- [3,15,105,945,10395,135135,2027025,34459425,654729075,13749310575]

denominadores :: [Integer]

denominadores = scanl (*) 3 [5, 7..]

-- 4ª solución

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

aproximaPi4 :: Integer -> Double

aproximaPi4 n =

read (x : "." ++ xs)

where (x:xs) = show (aproximaPi4' n)

aproximaPi4' :: Integer -> Integer

aproximaPi4' n =

2 * (p + sum (zipWith div (map (*p) numeradores) (genericTake n denominadores)))

where p = 10^n

-- Comprobación de equivalencia

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

-- La propiedad es

prop_aproximaPi :: Property

prop_aproximaPi =

forAll (arbitrary `suchThat` (> 3)) $ \n ->

all (=~ aproximaPi1 n)

[aproximaPi2 n,

aproximaPi3 n,

aproximaPi4 n]

(=~) :: Double -> Double -> Bool

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

-- La comprobación es

-- λ> quickCheck prop_aproximaPi

-- +++ OK, passed 100 tests.

-- Comparación de eficiencia

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

-- La comparación es

-- λ> aproximaPi1 3000

-- 3.141592653589793

-- (4.96 secs, 27,681,824,408 bytes)

-- λ> aproximaPi2 3000

-- 3.1415926535897922

-- (3.00 secs, 20,496,194,496 bytes)

-- λ> aproximaPi3 3000

-- 3.141592653589793

-- (3.13 secs, 13,439,528,432 bytes)

-- λ> aproximaPi4 3000

-- 3.141592653589793

-- (0.09 secs, 23,142,144 bytes)

El código se encuentra en GitHub.

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 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.

Aproximación del número pi

Soluciones

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

Soluciones

Cálculo aproximado de integrales definidas

Soluciones

Entradas recientes

Buscar

Correo electrónico