Menu Close

I1M2018: Cálculo numérico en Haskell

En la tercera parte de la clase de hoy de Informática de 1º del Grado en Matemáticas se han explicado las soluciones de los ejercicios de la relación 19, en la que se definen funciones para resolver los siguientes problemas de cálculo numérico:

  • diferenciación numérica,
  • cálculo de la raíz cuadrada mediante el método de Herón,
  • cálculo de los ceros de una función por el método de Newton y
  • cálculo de funciones inversas.

Un aspecto a destacar desde el punto de vista de la programación es el uso de la abstracción de procedimientos.

Los ejercicios, y sus soluciones, se muestran a continuación.

-- ---------------------------------------------------------------------
-- Introducción                                                       --
-- ---------------------------------------------------------------------
 
-- En esta relación se definen funciones para resolver los siguientes
-- problemas de cálculo numérico:
-- + diferenciación numérica,
-- + cálculo de la raíz cuadrada mediante el método de Herón,
-- + cálculo de los ceros de una función por el método de Newton y
-- + cálculo de funciones inversas.
 
-- ---------------------------------------------------------------------
-- Importación de librerías                                           --
-- ---------------------------------------------------------------------
 
import Test.QuickCheck
 
-- ---------------------------------------------------------------------
-- Diferenciación numérica                                            --
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 1.1. Definir la función
--    derivada :: Double -> (Double -> Double) -> Double -> Double
-- tal que (derivada a f x) es el valor de la derivada de la función f
-- en el punto x con aproximación a. Por ejemplo, 
--    derivada 0.001 sin pi  ==  -0.9999998333332315
--    derivada 0.001 cos pi  ==  4.999999583255033e-4
-- ---------------------------------------------------------------------
 
derivada :: Double -> (Double -> Double) -> Double -> Double
derivada a f x = (f (x+a) - f x)/a  
 
-- ---------------------------------------------------------------------
-- Ejercicio 1.2. Definir las funciones
--    derivadaBurda :: (Double -> Double) -> Double -> Double
--    derivadaFina  :: (Double -> Double) -> Double -> Double
--    derivadaSuper :: (Double -> Double) -> Double -> Double
-- tales que 
--    * (derivadaBurda f x) es el valor de la derivada de la función f 
--      en el punto x con aproximación 0.01,
--    * (derivadaFina f x) es el valor de la derivada de la función f 
--      en el punto x con aproximación 0.0001.
--    * (derivadaSuper f x) es el valor de la derivada de la función f 
--      en el punto x con aproximación 0.000001.
-- Por ejemplo,
--    derivadaBurda cos pi  ==  4.999958333473664e-3
--    derivadaFina  cos pi  ==  4.999999969612645e-5
--    derivadaSuper cos pi  ==  5.000444502911705e-7
-- ---------------------------------------------------------------------
 
derivadaBurda :: (Double -> Double) -> Double -> Double
derivadaBurda = derivada 0.01  
 
derivadaFina :: (Double -> Double) -> Double -> Double
derivadaFina  = derivada 0.0001  
 
derivadaSuper :: (Double -> Double) -> Double -> Double
derivadaSuper = derivada 0.000001  
 
-- ---------------------------------------------------------------------
-- Ejercicio 1.3. Definir la función
--    derivadaFinaDelSeno :: Double -> Double
-- tal que (derivadaFinaDelSeno x) es el valor de la derivada fina del
-- seno en x. Por ejemplo,
--    derivadaFinaDelSeno pi  ==  -0.9999999983354436
-- ---------------------------------------------------------------------
 
derivadaFinaDelSeno :: Double -> Double
derivadaFinaDelSeno = derivadaFina sin  
 
-- ---------------------------------------------------------------------
-- Cálculo de la raíz cuadrada                                        --
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 2.1. En los siguientes apartados de este ejercicio se va a
-- calcular la raíz cuadrada de un número basándose en las siguientes
-- propiedades:
-- + Si y es una aproximación de la raíz cuadrada de x, entonces 
--   (y+x/y)/2 es una aproximación mejor. 
-- + El límite de la sucesión definida por 
--       x_0     = 1 
--       x_{n+1} = (x_n+x/x_n)/2
--   es la raíz cuadrada de x.
-- 
-- Definir, por recursión, la función
--    raiz :: Double -> Double
-- tal que (raiz x) es la raíz cuadrada de x calculada usando la
-- propiedad anterior con una aproximación de 0.00001 y tomando como
-- valor inicial 1. Por ejemplo, 
--    raiz 9  ==  3.000000001396984
-- ---------------------------------------------------------------------
 
raiz :: Double -> Double
raiz x = raizAux 1
  where raizAux y | aceptable y = y
                  | otherwise   = raizAux (mejora y)
        aceptable y = abs(y*y-x) < 0.00001
        mejora y    = 0.5*(y+x/y)
 
-- ---------------------------------------------------------------------
-- Ejercicio 3.2. Definir el operador 
--    (~=) :: Double -> Double -> Bool
-- tal que (x ~= y) si |x-y| < 0.001. Por ejemplo,
--    3.05 ~= 3.07        ==  False
--    3.00005 ~= 3.00007  == True
-- ---------------------------------------------------------------------
 
infix 5 ~=
(~=) :: Double -> Double -> Bool
x ~= y = abs (x-y) < 0.001
 
-- ---------------------------------------------------------------------
-- Ejercicio 3.3. Comprobar con QuickCheck que si x es positivo,
-- entonces 
--    (raiz x)^2 ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_raiz :: Double -> Bool
prop_raiz x =
  raiz x' ^ 2 ~= x'
  where x' = abs x
 
-- La comprobación es
--    ghci> quickCheck prop_raiz
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ejercicio 3.4. Definir por recursión la función
--    until' :: (a -> Bool) -> (a -> a) -> a -> a
-- tal que (until' p f x) es el resultado de aplicar la función f a x el
-- menor número posible de veces, hasta alcanzar un valor que satisface
-- el predicado p. Por ejemplo, 
--    until' (>1000) (2*) 1  ==  1024
--
-- Nota: until' es equivalente a la predefinida until.
-- ---------------------------------------------------------------------
 
until' :: (a -> Bool) -> (a -> a) -> a -> a
until' p f x | p x       = x
             | otherwise = until' p f (f x)
 
-- ---------------------------------------------------------------------
-- Ejercicio 3.5. Definir, por iteración con until, la función
--    raizI :: Double -> Double
-- tal que (raizI x) es la raíz cuadrada de x calculada usando la
-- propiedad anterior. Por ejemplo, 
--    raizI 9  ==  3.000000001396984
-- ---------------------------------------------------------------------
 
raizI :: Double -> Double
raizI x = until aceptable mejora 1
  where aceptable y = abs(y*y-x) < 0.00001
        mejora y    = 0.5*(y+x/y)
 
-- ---------------------------------------------------------------------
-- Ejercicio 3.6. Comprobar con QuickCheck que si x es positivo,
-- entonces 
--    (raizI x)^2 ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_raizI :: Double -> Bool
prop_raizI x =
   raizI x' ^ 2 ~= x'
  where x' = abs x           
 
-- La comprobación es
--    ghci> quickCheck prop_raizI
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ceros de una función                                               --
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 4. Los ceros de una función pueden calcularse mediante el
-- método de Newton basándose en las siguientes propiedades:
-- + Si b es una aproximación para el punto cero de f, entonces 
--   b-f(b)/f'(b) es una mejor aproximación.
-- + el límite de la sucesión x_n definida por
--      x_0     = 1 
--      x_{n+1} = x_n-f(x_n)/f'(x_n)
--   es un cero de f.
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 4.1. Definir, por recursión, la función
--    puntoCero :: (Double -> Double) -> Double
-- tal que (puntoCero f) es un cero de la función f calculado usando la
-- propiedad anterior. Por ejemplo, 
--    puntoCero cos  ==  1.5707963267949576
-- ---------------------------------------------------------------------
 
puntoCero :: (Double -> Double) -> Double
puntoCero f = puntoCeroAux f 1
  where puntoCeroAux f' x | aceptable x = x
                          | otherwise   = puntoCeroAux f' (mejora x)
        aceptable b = abs (f b) < 0.00001
        mejora b    = b - f b / derivadaFina f b
 
-- ---------------------------------------------------------------------
-- Ejercicio 4.2. Definir, por iteración con until, la función
--    puntoCeroI :: (Double -> Double) -> Double
-- tal que (puntoCeroI f) es un cero de la función f calculado usando la
-- propiedad anterior. Por ejemplo, 
--    puntoCeroI cos  ==  1.5707963267949576
-- ---------------------------------------------------------------------
 
puntoCeroI :: (Double -> Double) -> Double
puntoCeroI f = until aceptable mejora 1
  where aceptable b = abs (f b) < 0.00001
        mejora b    = b - f b / derivadaFina f b
 
-- ---------------------------------------------------------------------
-- Funciones inversas                                                 --
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 5. En este ejercicio se usará la función puntoCero para
-- definir la inversa de distintas funciones.
-- ---------------------------------------------------------------------
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.1. Definir, usando puntoCero, la función
--    raizCuadrada :: Double -> Double
-- tal que (raizCuadrada x) es la raíz cuadrada de x. Por ejemplo,
--    raizCuadrada 9  ==  3.000000002941184
-- ---------------------------------------------------------------------
 
raizCuadrada :: Double -> Double
raizCuadrada a = puntoCero f
  where f x = x*x-a  
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.2. Comprobar con QuickCheck que si x es positivo,
-- entonces 
--    (raizCuadrada x)^2 ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_raizCuadrada :: Double -> Bool
prop_raizCuadrada x =
  raizCuadrada x' ^ 2 ~= x'
  where x' = abs x
 
-- La comprobación es
--    ghci> quickCheck prop_raizCuadrada
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.3. Definir, usando puntoCero, la función
--    raizCubica :: Double -> Double
-- tal que (raizCubica x) es la raíz cúbica de x. Por ejemplo,
--    raizCubica 27  ==  3.0000000000196048
-- ---------------------------------------------------------------------
 
raizCubica :: Double -> Double
raizCubica a = puntoCero f
  where f x = x*x*x-a 
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.4. Comprobar con QuickCheck que si x es positivo,
-- entonces 
--    (raizCubica x)^3 ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_raizCubica :: Double -> Bool
prop_raizCubica x =
  raizCubica x ^ 3 ~= x
 
-- La comprobación es
--    ghci> quickCheck prop_raizCubica
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.5. Definir, usando puntoCero, la función
--    arcoseno :: Double -> Double
-- tal que (arcoseno x) es el arcoseno de x. Por ejemplo,
--    arcoseno 1  == 1.5665489428306574
-- ---------------------------------------------------------------------
 
arcoseno :: Double -> Double
arcoseno a = puntoCero f
  where f x = sin x - a  
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.6. Comprobar con QuickCheck que si x está entre 0 y 1,
-- entonces 
--    sin (arcoseno x) ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_arcoseno :: Double -> Bool
prop_arcoseno x =
  sin (arcoseno x') ~= x'
  where x' = abs (x - fromIntegral (truncate x))
 
-- La comprobación es
--    ghci> quickCheck prop_arcoseno
--    OK, passed 100 tests.
 
-- Otra forma de expresar la propiedad es
prop_arcoseno2 :: Property
prop_arcoseno2 =
  forAll (choose (0,1)) $ \x -> sin (arcoseno x) ~= x
 
-- La comprobación es
--    ghci> quickCheck prop_arcoseno2
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.7. Definir, usando puntoCero, la función
--    arcocoseno :: Double -> Double
-- tal que (arcoseno x) es el arcoseno de x. Por ejemplo,
--    arcocoseno 0  == 1.5707963267949576
-- ---------------------------------------------------------------------
 
arcocoseno :: Double -> Double
arcocoseno a = puntoCero f
  where f x = cos x - a  
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.8. Comprobar con QuickCheck que si x está entre 0 y 1,
-- entonces 
--    cos (arcocoseno x) ~= x
-- ---------------------------------------------------------------------
 
-- La propiedad es
prop_arcocoseno :: Double -> Bool
prop_arcocoseno x =
  cos (arcocoseno x') ~= x'
  where x' = abs (x - fromIntegral (truncate x))
 
-- La comprobación es
--    ghci> quickCheck prop_arcocoseno
--    OK, passed 100 tests.
 
-- Otra forma de expresar la propiedad es
prop_arcocoseno2 :: Property
prop_arcocoseno2 =
  forAll (choose (0,1)) $ \x -> cos (arcocoseno x) ~= x
 
-- La comprobación es
--    ghci> quickCheck prop_arcocoseno2
--    OK, passed 100 tests.
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.7. Definir, usando puntoCero, la función
--    inversa :: (Double -> Double) -> Double -> Double
-- tal que (inversa g x) es el valor de la inversa de g en x. Por
-- ejemplo, 
--    inversa (^2) 9  ==  3.000000002941184
-- ---------------------------------------------------------------------
 
inversa :: (Double -> Double) -> Double -> Double
inversa g a = puntoCero f
  where f x = g x - a  
 
-- ---------------------------------------------------------------------
-- Ejercicio 5.8. Redefinir, usando inversa, las funciones raizCuadrada,
-- raizCubica, arcoseno y arcocoseno.
-- ---------------------------------------------------------------------
 
raizCuadrada', raizCubica', arcoseno', arcocoseno' :: Double -> Double
raizCuadrada' = inversa (^2)
raizCubica'   = inversa (^3)
arcoseno'     = inversa sin  
arcocoseno'   = inversa cos

Finalmente, se ha propuesto la resolución de la relación 20 (para ampliar los problemas de cálculo numérico) y la relación 21 (para programar las soluciones de los problemas anteriores en Maxima).

I1M2018