{"id":5959,"date":"2018-03-07T17:31:17","date_gmt":"2018-03-07T16:31:17","guid":{"rendered":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/?p=5959"},"modified":"2018-03-11T08:38:35","modified_gmt":"2018-03-11T07:38:35","slug":"i1m2017-calculo-numerico-en-haskell-2o-parte","status":"publish","type":"post","link":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/i1m2017-calculo-numerico-en-haskell-2o-parte\/","title":{"rendered":"I1M2017: C\u00e1lculo num\u00e9rico en Haskell (2\u00ba parte)"},"content":{"rendered":"<p>En la segunda parte de la clase de hoy de <a href=\"http:\/\/www.cs.us.es\/~jalonso\/cursos\/i1m-17\">Inform\u00e1tica de 1\u00ba del Grado en Matem\u00e1ticas<\/a> se han explicado las soluciones de los ejercicios de la relaci\u00f3n 27, en la que se definen funciones para resolver los siguientes problemas de c\u00e1lculo num\u00e9rico:<\/p>\n<ul>\n<li>C\u00e1lculo de l\u00edmites.<\/li>\n<li>C\u00e1lculo de los ceros de una funci\u00f3n por el m\u00e9todo de la bisecci\u00f3n.<\/li>\n<li>C\u00e1lculo de ra\u00edces enteras.<\/li>\n<li>C\u00e1lculo de integrales por el m\u00e9todo de los rect\u00e1ngulos.<\/li>\n<li>Algoritmo de bajada para resolver un sistema triangular inferior.<\/li>\n<\/ul>\n<p>Los ejercicios, y sus soluciones, se muestran a continuaci\u00f3n.<br \/>\n<!--more--><\/p>\n<pre lang=\"haskell\">\n-- ---------------------------------------------------------------------\n-- \u00a7 Librer\u00edas auxiliares                                             --\n-- ---------------------------------------------------------------------\n\nimport Test.QuickCheck\nimport Data.Matrix\n\n-- ---------------------------------------------------------------------\n-- \u00a7 C\u00e1lculo de l\u00edmites                                               --\n-- ---------------------------------------------------------------------\n\n-- ---------------------------------------------------------------------\n-- Ejercicio 1. Definir la funci\u00f3n  \n--    limite :: (Double -> Double) -> Double -> Double\n-- tal que (limite f a) es el valor de f en el primer t\u00e9rmino x tal que, \n-- para todo y entre x+1 y x+100, el valor absoluto de la diferencia\n-- entre f(y) y f(x) es menor que a. Por ejemplo,\n--    limite (\\n -> (2*n+1)\/(n+5)) 0.001  ==  1.9900110987791344\n--    limite (\\n -> (1+1\/n)**n) 0.001     ==  2.714072874546881\n-- ---------------------------------------------------------------------\n\nlimite :: (Double -> Double) -> Double -> Double\nlimite f a = \n  head [f x | x <- [1..],\n              maximum [abs (f y - f x) | y <- [x+1..x+100]] < a]\n\n-- ---------------------------------------------------------------------\n-- \u00a7 Ceros de una funci\u00f3n por el m\u00e9todo de la bisecci\u00f3n               --\n-- ---------------------------------------------------------------------\n\n-- ---------------------------------------------------------------------\n-- Ejercicio 2. El m\u00e9todo de bisecci\u00f3n para calcular un cero de una\n-- funci\u00f3n en el intervalo [a,b] se basa en el teorema de Bolzano: \n--    \"Si f(x) es una funci\u00f3n continua en el intervalo [a, b], y si,\n--    adem\u00e1s, en los extremos del intervalo la funci\u00f3n f(x) toma valores\n--    de signo opuesto (f(a) * f(b) < 0), entonces existe al menos un\n--    valor c en (a, b) para el que f(c) = 0\".\n--\n-- El m\u00e9todo para calcular un cero de la funci\u00f3n f en el intervalo [a,b]\n-- con un error menor que e consiste en tomar el punto medio del\n-- intervalo c = (a+b)\/2 y considerar los siguientes casos:\n-- (*) Si |f(c)| < e, hemos encontrado una aproximaci\u00f3n del punto que\n--     anula f en el intervalo con un error aceptable.\n-- (*) Si f(c) tiene signo distinto de f(a), repetir el proceso en el\n--     intervalo [a,c].\n-- (*) Si no, repetir el proceso en el intervalo [c,b].\n-- \n-- Definir la funci\u00f3n\n--    biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double\n-- tal que (biseccion f a b e) es una aproximaci\u00f3n del punto del\n-- intervalo [a,b] en el que se anula la funci\u00f3n f, con un error menor\n-- que e, calculada mediante el m\u00e9todo de la bisecci\u00f3n. Por ejemplo,\n--    biseccion (\\x -> x^2 - 3) 0 5 0.01             ==  1.7333984375\n--    biseccion (\\x -> x^3 - x - 2) 0 4 0.01         ==  1.521484375\n--    biseccion cos 0 2 0.01                         ==  1.5625\n--    biseccion (\\x -> log (50-x) - 4) (-10) 3 0.01  ==  -5.125\n-- ---------------------------------------------------------------------\n\n-- 1\u00aa soluci\u00f3n\nbiseccion :: (Double -> Double) -> Double -> Double -> Double -> Double\nbiseccion f a b e  \n  | abs (f c) < e   = c\n  | (f a)*(f c) < 0 = biseccion f a c e\n  | otherwise       = biseccion f c b e\n  where c = (a+b)\/2\n\n-- 2\u00aa soluci\u00f3n\nbiseccion2 :: (Double -> Double) -> Double -> Double -> Double -> Double\nbiseccion2 f a b e = aux a b\n  where aux a b | abs (f c) < e   = c\n                | (f a)*(f c) < 0 = aux a c \n                | otherwise       = aux c b\n          where c = (a+b)\/2\n\n-- ---------------------------------------------------------------------\n-- \u00a7 C\u00e1lculo de ra\u00edces enteras                                        --\n-- ---------------------------------------------------------------------\n\n-- ---------------------------------------------------------------------\n-- Ejercicio 3. Definir la funci\u00f3n \n--    raizEnt :: Integer -> Integer -> Integer\n-- tal que (raizEnt x n) es la ra\u00edz entera n-\u00e9sima de x; es decir, el\n-- mayor n\u00famero entero y tal que y^n <= x. Por ejemplo,\n--    raizEnt  8 3      ==  2\n--    raizEnt  9 3      ==  2\n--    raizEnt 26 3      ==  2\n--    raizEnt 27 3      ==  3\n--    raizEnt (10^50) 2 ==  10000000000000000000000000\n--\n-- Comprobar con QuickCheck que para todo n\u00famero natural n, \n--     raizEnt (10^(2*n)) 2 == 10^n\n-- ---------------------------------------------------------------------\n\n-- 1\u00aa definici\u00f3n\nraizEnt1 :: Integer -> Integer -> Integer\nraizEnt1 x n =\n    last (takeWhile (\\y -> y^n <= x) [0..])\n\n-- 2\u00aa definici\u00f3n         \nraizEnt2 :: Integer -> Integer -> Integer\nraizEnt2 x n =\n    floor ((fromIntegral x)**(1 \/ fromIntegral n))\n\n-- Nota. La definici\u00f3n anterior falla para n\u00fameros grandes. Por ejemplo,\n--    \u03bb> raizEnt2 (10^50) 2 == 10^25\n--    False\n          \n-- 3\u00aa definici\u00f3n          \nraizEnt3 :: Integer -> Integer -> Integer\nraizEnt3 x n = aux (1,x)\n  where aux (a,b) | d == x    = c\n                  | c == a    = c\n                  | d < x     = aux (c,b)\n                  | otherwise = aux (a,c) \n          where c = (a+b) `div` 2\n                d = c^n\n\n-- Comparaci\u00f3n de eficiencia\n--    \u03bb> raizEnt1 (10^14) 2\n--    10000000\n--    (6.15 secs, 6,539,367,976 bytes)\n--    \u03bb> raizEnt2 (10^14) 2\n--    10000000\n--    (0.00 secs, 0 bytes)\n--    \u03bb> raizEnt3 (10^14) 2\n--    10000000\n--    (0.00 secs, 25,871,944 bytes)\n--    \n--    \u03bb> raizEnt2 (10^50) 2\n--    9999999999999998758486016\n--    (0.00 secs, 0 bytes)\n--    \u03bb> raizEnt3 (10^50) 2\n--    10000000000000000000000000\n--    (0.00 secs, 0 bytes)\n                        \n-- La propiedad es                        \nprop_raizEnt :: Integer -> Bool\nprop_raizEnt n =\n    raizEnt3 (10^(2*m)) 2 == 10^m\n    where m = abs n\n\n-- La comprobaci\u00f3n es              \n--    \u03bb> quickCheck prop_raizEnt\n--    +++ OK, passed 100 tests.\n\n-- ---------------------------------------------------------------------\n-- \u00a7 Integraci\u00f3n por el m\u00e9todo de los rect\u00e1ngulos                     --\n-- ---------------------------------------------------------------------\n\n-- ---------------------------------------------------------------------\n-- Ejercicio 4. La integral definida de una funci\u00f3n f entre los l\u00edmites\n-- a y b puede calcularse mediante la regla del rect\u00e1ngulo\n-- (ver en http:\/\/bit.ly\/1FDhZ1z) usando la f\u00f3rmula  \n--    h * (f(a+h\/2) + f(a+h+h\/2) + f(a+2h+h\/2) + ... + f(a+n*h+h\/2))\n-- con a+n*h+h\/2 <= b < a+(n+1)*h+h\/2 y usando valores peque\u00f1os para h.\n-- \n-- Definir la funci\u00f3n\n--    integral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a\n-- tal que (integral a b f h) es el valor de dicha expresi\u00f3n. Por\n-- ejemplo, el c\u00e1lculo de la integral de f(x) = x^3 entre 0 y 1, con\n-- paso 0.01, es \n--    integral 0 1 (^3) 0.01  ==  0.24998750000000042\n-- Otros ejemplos son\n--    integral 0 1 (^4) 0.01                        ==  0.19998333362500048\n--    integral 0 1 (\\x -> 3*x^2 + 4*x^3) 0.01       ==  1.9999250000000026\n--    log 2 - integral 1 2 (\\x -> 1\/x) 0.01         ==  3.124931644782336e-6\n--    pi - 4 * integral 0 1 (\\x -> 1\/(x^2+1)) 0.01  ==  -8.333333331389525e-6\n-- ---------------------------------------------------------------------\n\n-- 1\u00aa soluci\u00f3n\n-- ===========\n\nintegral :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a\nintegral a b f h = h * suma (a+h\/2) b (+h) f\n\n-- (suma a b s f) es l valor de\n--    f(a) + f(s(a)) + f(s(s(a)) + ... + f(s(...(s(a))...))\n-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,\n--    suma 2 5 (1+) (^3)  ==  224\nsuma :: (Ord t, Num a) => t -> t -> (t -> t) -> (t -> a) -> a\nsuma a b s f = sum [f x | x <- sucesion a b s]\n\n-- (sucesion x y s) es la lista\n--    [a, s(a), s(s(a), ..., s(...(s(a))...)]\n-- hasta que s(s(...(s(a))...)) > b. Por ejemplo,\n--    sucesion 3 20 (+2)  ==  [3,5,7,9,11,13,15,17,19]\nsucesion :: Ord a => a -> a -> (a -> a) -> [a]\nsucesion a b s = takeWhile (<=b) (iterate s a)\n\n-- 2\u00aa soluci\u00f3n\n-- ===========\n\nintegral2 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a\nintegral2 a b f h\n    | a+h\/2 > b = 0\n    | otherwise = h * f (a+h\/2) + integral2 (a+h) b f h\n\n-- 3\u00aa soluci\u00f3n\n-- ===========\n\nintegral3 :: (Fractional a, Ord a) => a -> a -> (a -> a) -> a -> a\nintegral3 a b f h = aux a where\n    aux x | x+h\/2 > b = 0\n          | otherwise = h * f (x+h\/2) + aux (x+h)\n\n-- Comparaci\u00f3n de eficiencia\n--    ghci> integral 0 10 (^3) 0.00001\n--    2499.9999998811422\n--    (4.62 secs, 1084774336 bytes)\n--    ghci> integral2 0 10 (^3) 0.00001\n--    2499.999999881125\n--    (7.90 secs, 1833360768 bytes)\n--    ghci> integral3 0 10 (^3) 0.00001\n--    2499.999999881125\n--    (7.27 secs, 1686056080 bytes)\n\n-- ---------------------------------------------------------------------\n-- \u00a7 Algoritmo de bajada para resolver un sistema triangular inferior --\n-- ---------------------------------------------------------------------\n\n-- ---------------------------------------------------------------------\n-- Ejercicio 5. Un sistema de ecuaciones lineales Ax = b es triangular\n-- inferior si todos los elementos de la matriz A que est\u00e1n por encima\n-- de la diagonal principal son nulos; es decir, es de la forma\n--    a(1,1)*x(1)                                               = b(1)\n--    a(2,1)*x(1) + a(2,2)*x(2)                                 = b(2)\n--    a(3,1)*x(1) + a(3,2)*x(2) + a(3,3)*x(3)                   = b(3)\n--    ...\n--    a(n,1)*x(1) + a(n,2)*x(2) + a(n,3)*x(3) +...+ a(x,x)*x(n) = b(n)\n--\n-- El sistema es compatible si, y s\u00f3lo si, el producto de los elementos\n-- de la diagonal principal es distinto de cero. En este caso, la\n-- soluci\u00f3n se puede calcular mediante el algoritmo de bajada: \n--    x(1) = b(1) \/ a(1,1)\n--    x(2) = (b(2) - a(2,1)*x(1)) \/ a(2,2)\n--    x(3) = (b(3) - a(3,1)*x(1) - a(3,2)*x(2)) \/ a(3,3)\n--    ...\n--    x(n) = (b(n) - a(n,1)*x(1) - a(n,2)*x(2) -...- a(n,n-1)*x(n-1)) \/ a(n,n)\n-- \n-- Definir la funci\u00f3n \n--    bajada :: Matrix Double -> Matrix Double -> Matrix Double\n-- tal que (bajada a b) es la soluci\u00f3n, mediante el algoritmo de bajada,\n-- del sistema compatible triangular superior ax = b. Por ejemplo, \n--    ghci> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]\n--    ghci> let b = fromLists [[3],[6.5],[10]]\n--    ghci> bajada a b\n--    ( 1.5 )\n--    ( 2.0 )\n--    ( 0.0 )\n-- Es decir, la soluci\u00f3n del sistema\n--    2x            = 3\n--    3x + y        = 6.5\n--    4x + 2y + 5 z = 10\n-- es x=1.5, y=2 y z=0.\n-- ---------------------------------------------------------------------\n\nbajada :: Matrix Double -> Matrix Double -> Matrix Double\nbajada a b = fromLists [[x i] | i <- [1..m]]\n    where m = nrows a\n          x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) \/ a!(k,k)\n<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>En la segunda parte de la clase de hoy de Inform\u00e1tica de 1\u00ba del Grado en Matem\u00e1ticas se han explicado las soluciones de los ejercicios de la relaci\u00f3n 27, en la que se definen funciones para resolver los siguientes problemas de c\u00e1lculo num\u00e9rico: C\u00e1lculo de l\u00edmites. C\u00e1lculo de los ceros de una funci\u00f3n por el&#8230;<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"_kad_post_transparent":"","_kad_post_title":"","_kad_post_layout":"","_kad_post_sidebar_id":"","_kad_post_content_style":"","_kad_post_vertical_padding":"","_kad_post_feature":"","_kad_post_feature_position":"","_kad_post_header":false,"_kad_post_footer":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"footnotes":"","_jetpack_memberships_contains_paid_content":false},"categories":[265],"tags":[270,316],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_likes_enabled":false,"_links":{"self":[{"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/posts\/5959"}],"collection":[{"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/comments?post=5959"}],"version-history":[{"count":2,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/posts\/5959\/revisions"}],"predecessor-version":[{"id":5961,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/posts\/5959\/revisions\/5961"}],"wp:attachment":[{"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/media?parent=5959"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/categories?post=5959"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.glc.us.es\/~jalonso\/vestigium\/wp-json\/wp\/v2\/tags?post=5959"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}