Exercitium

Algoritmo de bajada para resolver un sistema triangular inferior

PorJosé A. Alonso 4-diciembre-202321-noviembre-2023

Un sistema de ecuaciones lineales $Ax = b$ es triangular inferior si todos los elementos de la matriz $A$ que están por encima de la diagonal principal son nulos; es decir, es de la forma
\begin{align}
&a_{1 1}x_1 &= b_1 \\
&a_{2 1}x_1 + a_{2 2}x_2 &= b_2 \\
&a_{3 1}x_1 + a_{3 2}x_2 + a_{3 3}x_3 &= b_3 \\
&… & \\
&a_{n 1}x_1 + a_{n 2}x_2 + a_{n 3}x_3 +…+ a_{n n}x_n &= b_n
\end{align}

El sistema es compatible si, y sólo si, el producto de los elementos de la diagonal principal es distinto de cero. En este caso, la solución se puede calcular mediante el algoritmo de bajada:
\begin{align}
x_1 &= \frac{b_1}{a_{1 1}} \\
x_2 &= \frac{b_2 – a_{2 1}x_1}{a_{2 2}} \\
x_3 &= \frac{b_3 – a_{3 1}x_1 – a_{3 2}x_2}{a_{3 3}} \\
… \\
x_n &= \frac{b_n – a_{n 1}x_1 – a_{n 2}x_2 – … – a_{n,n-1}x_{n-1}}{a_{n n}}
\end{align}

Definir la función

   bajada :: Matrix Double -> Matrix Double -> Matrix Double

1	bajada :: Matrix Double -> Matrix Double -> Matrix Double

tal que bajada a b es la solución, mediante el algoritmo de bajada, del sistema compatible triangular superior ax = b. Por ejemplo,

   λ> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]
   λ> let b = fromLists [[3],[6.5],[10]]
   λ> bajada a b
   ( 1.5 )
   ( 2.0 )
   ( 0.0 )

λ> let a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]

λ> let b = fromLists [[3],[6.5],[10]]

λ> bajada a b

( 1.5 )

( 2.0 )

( 0.0 )

Es decir, la solución del sistema
\begin{align}
2x &= 3 \\
3x + y &= 6.5 \\
4x + 2y + 5z &= 10
\end{align}
es $x=1.5$, $y=2$ y $z=0$.

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

module Algoritmo_de_bajada where

import Data.Matrix (Matrix, (!), fromLists, nrows, toLists)
import Test.Hspec (Spec, hspec, it, shouldBe)

bajada :: Matrix Double -> Matrix Double -> Matrix Double
bajada a b = fromLists [[x i] | i <- [1..m]]
  where m   = nrows a
        x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "ej1" $
    toLists (bajada a b) `shouldBe` [[1.5],[2.0],[0.0]]
    where
      a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]
      b = fromLists [[3],[6.5],[10]]

-- La verificación es
--    λ> verifica
--
--    Finished in 0.0007 seconds
--    1 example, 0 failures

module Algoritmo_de_bajada where

import Data.Matrix (Matrix, (!), fromLists, nrows, toLists)

import Test.Hspec (Spec, hspec, it, shouldBe)

bajada :: Matrix Double -> Matrix Double -> Matrix Double

bajada a b = fromLists [[x i] | i <- [1..m]]

where m = nrows a

x k = (b!(k,1) - sum [a!(k,j) * x j | j <- [1..k-1]]) / a!(k,k)

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "ej1" $

toLists (bajada a b) `shouldBe` [[1.5],[2.0],[0.0]]

where

a = fromLists [[2,0,0],[3,1,0],[4,2,5.0]]

b = fromLists [[3],[6.5],[10]]

-- La verificación es

-- λ> verifica

-- Finished in 0.0007 seconds

-- 1 example, 0 failures

Soluciones en Python

def bajada(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:
    n = len(a)
    def x(k: int) -> float:
        return (b[k][0] - sum((a[k][j] * x(j) for j in range(0, k)))) / a[k][k]
    return [[x(i)] for i in range(0, n)]

# Verificación
# ============

def test_bajada() -> None:
    assert bajada([[2,0,0],[3,1,0],[4,2,5.0]], [[3],[6.5],[10]]) == \
        [[1.5], [2.0], [0.0]]
    print("Verificado")

# La verificación es
#    >>> test_bajada()
#    Verificado

def bajada(a: list[list[float]], b: list[list[float]]) -> list[list[float]]:

n = len(a)

def x(k: int) -> float:

return (b[k][0] - sum((a[k][j] * x(j) for j in range(0, k)))) / a[k][k]

return [[x(i)] for i in range(0, n)]

# Verificación

# ============

def test_bajada() -> None:

assert bajada([[2,0,0],[3,1,0],[4,2,5.0]], [[3],[6.5],[10]]) == \

[[1.5], [2.0], [0.0]]

print("Verificado")

# La verificación es

# >>> test_bajada()

# Verificado

Integración por el método de los rectángulos

PorJosé A. Alonso 29-noviembre-202320-noviembre-2023

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

   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

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

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

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

Soluciones en Python

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

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

Raíces enteras

PorJosé A. Alonso 24-noviembre-202313-noviembre-2023

Definir la función

   raizEnt :: Integer -> Integer -> Integer

1	raizEnt :: Integer -> Integer -> Integer

tal que raizEnt x n es la raíz entera n-ésima de x; es decir, el mayor número entero y tal que $y^n \leq x$. Por ejemplo,

   raizEnt  8 3      ==  2
   raizEnt  9 3      ==  2
   raizEnt 26 3      ==  2
   raizEnt 27 3      ==  3
   raizEnt (10^50) 2 ==  10000000000000000000000000

raizEnt 8 3 == 2

raizEnt 9 3 == 2

raizEnt 26 3 == 2

raizEnt 27 3 == 3

raizEnt (10^50) 2 == 10000000000000000000000000

Comprobar con QuickCheck que para todo número natural n,

    raizEnt (10^(2*n)) 2 == 10^n

1	raizEnt (10^(2*n)) 2 == 10^n

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

module Raices_enteras where

import Test.Hspec (Spec, hspec, it, shouldBe)
import Test.QuickCheck (quickCheck)

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

raizEnt1 :: Integer -> Integer -> Integer
raizEnt1 x n =
  last (takeWhile (\y -> y^n <= x) [0..])

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

raizEnt2 :: Integer -> Integer -> Integer
raizEnt2 x n =
  floor ((fromIntegral x)**(1 / fromIntegral n))

-- Nota. La solución anterior falla para números grandes. Por ejemplo,
--    λ> raizEnt2 (10^50) 2 == 10^25
--    False

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

raizEnt3 :: Integer -> Integer -> Integer
raizEnt3 x n = aux (1,x)
  where aux (a,b) | d == x    = c
                  | c == a    = c
                  | d < x     = aux (c,b)
                  | otherwise = aux (a,c)
          where c = (a+b) `div` 2
                d = c^n

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

--    λ> raizEnt1 (10^14) 2
--    10000000
--    (6.15 secs, 6,539,367,976 bytes)
--    λ> raizEnt2 (10^14) 2
--    10000000
--    (0.00 secs, 0 bytes)
--    λ> raizEnt3 (10^14) 2
--    10000000
--    (0.00 secs, 25,871,944 bytes)
--
--    λ> raizEnt2 (10^50) 2
--    9999999999999998758486016
--    (0.00 secs, 0 bytes)
--    λ> raizEnt3 (10^50) 2
--    10000000000000000000000000
--    (0.00 secs, 0 bytes)

-- Comprobación de la propiedad
-- ============================

-- La propiedad es
prop_raizEnt :: Integer -> Bool
prop_raizEnt n =
  raizEnt3 (10^(2*m)) 2 == 10^m
  where m = abs n

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

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    raizEnt1  8 3 `shouldBe` 2
  it "e2" $
    raizEnt1  9 3 `shouldBe` 2
  it "e3" $
    raizEnt1 26 3 `shouldBe` 2
  it "e4" $
    raizEnt1 27 3 `shouldBe` 3
  it "e5" $
    raizEnt2  8 3 `shouldBe` 2
  it "e6" $
    raizEnt2  9 3 `shouldBe` 2
  it "e7" $
    raizEnt2 26 3 `shouldBe` 2
  it "e8" $
    raizEnt2 27 3 `shouldBe` 3
  it "e9" $
    raizEnt3  8 3 `shouldBe` 2
  it "e10" $
    raizEnt3  9 3 `shouldBe` 2
  it "e11" $
    raizEnt3 26 3 `shouldBe` 2
  it "e12" $
    raizEnt3 27 3 `shouldBe` 3

-- La verificación es
--    λ> verifica
--
--    e1
--    e2
--    e3
--    e4
--    e5
--    e6
--    e7
--    e8
--    e9
--    e10
--    e11
--    e12
--
--    Finished in 0.0007 seconds
--    12 examples, 0 failures

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

module Raices_enteras where

import Test.Hspec (Spec, hspec, it, shouldBe)

import Test.QuickCheck (quickCheck)

-- 1ª solución

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

raizEnt1 :: Integer -> Integer -> Integer

raizEnt1 x n =

last (takeWhile (\y -> y^n <= x) [0..])

-- 2ª solución

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

raizEnt2 :: Integer -> Integer -> Integer

raizEnt2 x n =

floor ((fromIntegral x)**(1 / fromIntegral n))

-- Nota. La solución anterior falla para números grandes. Por ejemplo,

-- λ> raizEnt2 (10^50) 2 == 10^25

-- False

-- 3ª solución

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

raizEnt3 :: Integer -> Integer -> Integer

raizEnt3 x n = aux (1,x)

where aux (a,b) | d == x = c

| c == a = c

| d < x = aux (c,b)

| otherwise = aux (a,c)

where c = (a+b) `div` 2

d = c^n

-- Comparación de eficiencia

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

-- λ> raizEnt1 (10^14) 2

-- 10000000

-- (6.15 secs, 6,539,367,976 bytes)

-- λ> raizEnt2 (10^14) 2

-- 10000000

-- (0.00 secs, 0 bytes)

-- λ> raizEnt3 (10^14) 2

-- 10000000

-- (0.00 secs, 25,871,944 bytes)

-- λ> raizEnt2 (10^50) 2

-- 9999999999999998758486016

-- (0.00 secs, 0 bytes)

-- λ> raizEnt3 (10^50) 2

-- 10000000000000000000000000

-- (0.00 secs, 0 bytes)

-- Comprobación de la propiedad

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

-- La propiedad es

prop_raizEnt :: Integer -> Bool

prop_raizEnt n =

raizEnt3 (10^(2*m)) 2 == 10^m

where m = abs n

-- La comprobación es

-- λ> quickCheck prop_raizEnt

-- +++ OK, passed 100 tests.

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "e1" $

raizEnt1 8 3 `shouldBe` 2

it "e2" $

raizEnt1 9 3 `shouldBe` 2

it "e3" $

raizEnt1 26 3 `shouldBe` 2

it "e4" $

raizEnt1 27 3 `shouldBe` 3

it "e5" $

raizEnt2 8 3 `shouldBe` 2

it "e6" $

raizEnt2 9 3 `shouldBe` 2

it "e7" $

raizEnt2 26 3 `shouldBe` 2

it "e8" $

raizEnt2 27 3 `shouldBe` 3

it "e9" $

raizEnt3 8 3 `shouldBe` 2

it "e10" $

raizEnt3 9 3 `shouldBe` 2

it "e11" $

raizEnt3 26 3 `shouldBe` 2

it "e12" $

raizEnt3 27 3 `shouldBe` 3

-- La verificación es

-- λ> verifica

-- e1

-- e2

-- e3

-- e4

-- e5

-- e6

-- e7

-- e8

-- e9

-- e10

-- e11

-- e12

-- Finished in 0.0007 seconds

-- 12 examples, 0 failures

Soluciones en Python

from itertools import count, takewhile
from math import floor
from sys import setrecursionlimit
from timeit import Timer, default_timer

from hypothesis import given
from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª solución
# ===========

def raizEnt(x: int, n: int) -> int:
    return list(takewhile(lambda y : y ** n <= x, count(0)))[-1]

# 2ª solución
# ===========

def raizEnt2(x: int, n: int) -> int:
    return floor(x ** (1 / n))

# Nota. La solución anterior falla para números grandes. Por ejemplo,
#    >>> raizEnt2(10**50, 2) == 10 **25
#    False

# 3ª solución
# ===========

def raizEnt3(x: int, n: int) -> int:
    def aux(a: int, b: int) -> int:
        c = (a + b) // 2
        d = c ** n
        if d == x:
            return c
        if c == a:
            return c
        if d < x:
            return aux(c, b)
        return aux(a, c)
    return aux(1, x)

# Comparación de eficiencia
# =========================

def tiempo(e: str) -> None:
    """Tiempo (en segundos) de evaluar la expresión e."""
    t = Timer(e, "", default_timer, globals()).timeit(1)
    print(f"{t:0.2f} segundos")

# La comparación es
#    >>> tiempo('raizEnt(10**14, 2)')
#    2.71 segundos
#    >>> tiempo('raizEnt2(10**14, 2)')
#    0.00 segundos
#    >>> tiempo('raizEnt3(10**14, 2)')
#    0.00 segundos
#
#    >>> raizEnt2(10**50, 2)
#    10000000000000000905969664
#    >>> raizEnt3(10**50, 2)
#    10000000000000000000000000

# Comprobación de la propiedad
# ============================

# La propiedad es
@given(st.integers(min_value=0, max_value=1000))
def test_raizEntP(n: int) -> None:
    assert raizEnt3(10**(2*n), 2) == 10**n

# La comprobación es
#    >>> test_raizEnt)()
#    >>>

# Verificación
# ============

def test_raizEnt() -> None:
    assert raizEnt(8, 3) == 2
    assert raizEnt(9, 3) == 2
    assert raizEnt(26, 3) == 2
    assert raizEnt(27, 3) == 3
    assert raizEnt2(8, 3) == 2
    assert raizEnt2(9, 3) == 2
    assert raizEnt2(26, 3) == 2
    assert raizEnt2(27, 3) == 3
    assert raizEnt3(8, 3) == 2
    assert raizEnt3(9, 3) == 2
    assert raizEnt3(26, 3) == 2
    assert raizEnt3(27, 3) == 3
    print("Verificado")

# La comprobación es
#    >>> test_raizEnt()
#    Verificado

from itertools import count, takewhile

from math import floor

from sys import setrecursionlimit

from timeit import Timer, default_timer

from hypothesis import given

from hypothesis import strategies as st

setrecursionlimit(10**6)

# 1ª solución

# ===========

def raizEnt(x: int, n: int) -> int:

return list(takewhile(lambda y : y ** n <= x, count(0)))[-1]

# 2ª solución

# ===========

def raizEnt2(x: int, n: int) -> int:

return floor(x ** (1 / n))

# Nota. La solución anterior falla para números grandes. Por ejemplo,

# >>> raizEnt2(10**50, 2) == 10 **25

# False

# 3ª solución

# ===========

def raizEnt3(x: int, n: int) -> int:

def aux(a: int, b: int) -> int:

c = (a + b) // 2

d = c ** n

if d == x:

return c

if c == a:

return c

if d < x:

return aux(c, b)

return aux(a, c)

return aux(1, x)

# Comparación de eficiencia

# =========================

def tiempo(e: str) -> None:

"""Tiempo (en segundos) de evaluar la expresión e."""

t = Timer(e, "", default_timer, globals()).timeit(1)

print(f"{t:0.2f} segundos")

# La comparación es

# >>> tiempo('raizEnt(10**14, 2)')

# 2.71 segundos

# >>> tiempo('raizEnt2(10**14, 2)')

# 0.00 segundos

# >>> tiempo('raizEnt3(10**14, 2)')

# 0.00 segundos

# >>> raizEnt2(10**50, 2)

# 10000000000000000905969664

# >>> raizEnt3(10**50, 2)

# 10000000000000000000000000

# Comprobación de la propiedad

# ============================

# La propiedad es

@given(st.integers(min_value=0, max_value=1000))

def test_raizEntP(n: int) -> None:

assert raizEnt3(10**(2*n), 2) == 10**n

# La comprobación es

# >>> test_raizEnt)()

# >>>

# Verificación

# ============

def test_raizEnt() -> None:

assert raizEnt(8, 3) == 2

assert raizEnt(9, 3) == 2

assert raizEnt(26, 3) == 2

assert raizEnt(27, 3) == 3

assert raizEnt2(8, 3) == 2

assert raizEnt2(9, 3) == 2

assert raizEnt2(26, 3) == 2

assert raizEnt2(27, 3) == 3

assert raizEnt3(8, 3) == 2

assert raizEnt3(9, 3) == 2

assert raizEnt3(26, 3) == 2

assert raizEnt3(27, 3) == 3

print("Verificado")

# La comprobación es

# >>> test_raizEnt()

# Verificado

Método de bisección para calcular raíces de una función

PorJosé A. Alonso 19-noviembre-20238-noviembre-2023

El método de bisección para calcular una raíz 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$ 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 una raíz 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 = \frac{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

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

module Metodo_de_biseccion_para_calcular_ceros_de_una_funcion where

import Test.Hspec (Spec, hspec, it, shouldBe)

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

biseccion :: (Double -> Double) -> Double -> Double -> Double -> Double
biseccion f a b e
  | abs (f c) < e = c
  | f a * f c < 0 = biseccion f a c e
  | otherwise     = biseccion 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

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    biseccion (\x -> x^2 - 3) 0 5 0.01             `shouldBe`  1.7333984375
  it "e2" $
    biseccion (\x -> x^3 - x - 2) 0 4 0.01         `shouldBe`  1.521484375
  it "e3" $
    biseccion cos 0 2 0.01                         `shouldBe`  1.5625
  it "e4" $
    biseccion (\x -> log (50-x) - 4) (-10) 3 0.01  `shouldBe`  -5.125
  it "e5" $
    biseccion2 (\x -> x^2 - 3) 0 5 0.01             `shouldBe`  1.7333984375
  it "e6" $
    biseccion2 (\x -> x^3 - x - 2) 0 4 0.01         `shouldBe`  1.521484375
  it "e7" $
    biseccion2 cos 0 2 0.01                         `shouldBe`  1.5625
  it "e8" $
    biseccion2 (\x -> log (50-x) - 4) (-10) 3 0.01  `shouldBe`  -5.125

-- La verificación es
--    λ> verifica
--
--    e1
--    e2
--    e3
--    e4
--    e5
--    e6
--    e7
--    e8
--
--    Finished in 0.0008 seconds
--    8 examples, 0 failures

module Metodo_de_biseccion_para_calcular_ceros_de_una_funcion where

import Test.Hspec (Spec, hspec, it, shouldBe)

-- 1ª solución

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

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

biseccion f a b e

| abs (f c) < e = c

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

| otherwise = biseccion 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

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "e1" $

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

it "e2" $

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

it "e3" $

biseccion cos 0 2 0.01 `shouldBe` 1.5625

it "e4" $

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

it "e5" $

biseccion2 (\x -> x^2 - 3) 0 5 0.01 `shouldBe` 1.7333984375

it "e6" $

biseccion2 (\x -> x^3 - x - 2) 0 4 0.01 `shouldBe` 1.521484375

it "e7" $

biseccion2 cos 0 2 0.01 `shouldBe` 1.5625

it "e8" $

biseccion2 (\x -> log (50-x) - 4) (-10) 3 0.01 `shouldBe` -5.125

-- La verificación es

-- λ> verifica

-- e1

-- e2

-- e3

-- e4

-- e5

-- e6

-- e7

-- e8

-- Finished in 0.0008 seconds

-- 8 examples, 0 failures

Soluciones en Python

from math import cos, log
from typing import Callable

# 1ª solución
# ===========

def biseccion(f: Callable[[float], float],
              a: float,
              b: float,
              e: float) -> float:
    c = (a+b)/2
    if abs(f(c)) < e:
        return c
    if f(a) * f(c) < 0:
        return biseccion(f, a, c, e)
    return biseccion(f, c, b, e)

# 2ª solución
# ===========

def biseccion2(f: Callable[[float], float],
               a: float,
               b: float,
               e: float) -> float:
    def aux(a1: float,  b1: float) -> float:
        c = (a1+b1)/2
        if abs(f(c)) < e:
            return c
        if f(a1) * f(c) < 0:
            return aux(a1, c)
        return aux(c, b1)
    return aux(a, b)

# Verificación
# ============

def test_biseccion() -> None:
    assert biseccion(lambda x : x**2 - 3, 0, 5, 0.01) == 1.7333984375
    assert biseccion(lambda x : x**3 - x - 2, 0, 4, 0.01) == 1.521484375
    assert biseccion(cos, 0, 2, 0.01) == 1.5625
    assert biseccion(lambda x : log(50-x) - 4, -10, 3, 0.01) == -5.125
    assert biseccion2(lambda x : x**2 - 3, 0, 5, 0.01) == 1.7333984375
    assert biseccion2(lambda x : x**3 - x - 2, 0, 4, 0.01) == 1.521484375
    assert biseccion2(cos, 0, 2, 0.01) == 1.5625
    assert biseccion2(lambda x : log(50-x) - 4, -10, 3, 0.01) == -5.125
    print("Verificado")


# La comprobación es
#    >>> test_biseccion()
#    Verificado

from math import cos, log

from typing import Callable

# 1ª solución

# ===========

def biseccion(f: Callable[[float], float],

a: float,

b: float,

e: float) -> float:

c = (a+b)/2

if abs(f(c)) < e:

return c

if f(a) * f(c) < 0:

return biseccion(f, a, c, e)

return biseccion(f, c, b, e)

# 2ª solución

# ===========

def biseccion2(f: Callable[[float], float],

a: float,

b: float,

e: float) -> float:

def aux(a1: float, b1: float) -> float:

c = (a1+b1)/2

if abs(f(c)) < e:

return c

if f(a1) * f(c) < 0:

return aux(a1, c)

return aux(c, b1)

return aux(a, b)

# Verificación

# ============

def test_biseccion() -> None:

assert biseccion(lambda x : x**2 - 3, 0, 5, 0.01) == 1.7333984375

assert biseccion(lambda x : x**3 - x - 2, 0, 4, 0.01) == 1.521484375

assert biseccion(cos, 0, 2, 0.01) == 1.5625

assert biseccion(lambda x : log(50-x) - 4, -10, 3, 0.01) == -5.125

assert biseccion2(lambda x : x**2 - 3, 0, 5, 0.01) == 1.7333984375

assert biseccion2(lambda x : x**3 - x - 2, 0, 4, 0.01) == 1.521484375

assert biseccion2(cos, 0, 2, 0.01) == 1.5625

assert biseccion2(lambda x : log(50-x) - 4, -10, 3, 0.01) == -5.125

print("Verificado")

# La comprobación es

# >>> test_biseccion()

# Verificado

Límites de sucesiones

PorJosé A. Alonso 14-noviembre-20235-noviembre-2023

Definir la función

   limite :: (Double -> Double) -> Double -> Double

1	limite :: (Double -> Double) -> Double -> Double

tal que limite f a es el valor de f en el primer término x tal que, para todo y entre x+1 y x+100, el valor absoluto de la diferencia entre f(y) y f(x) es menor que a. Por ejemplo,

   limite (\n -> (2*n+1)/(n+5)) 0.001  ==  1.9900110987791344
   limite (\n -> (1+1/n)**n) 0.001     ==  2.714072874546881

1 2	limite (\n -> (2n+1)/(n+5)) 0.001 == 1.9900110987791344 limite (\n -> (1+1/n)*n) 0.001 == 2.714072874546881

Soluciones

A continuación se muestran las soluciones en Haskell y las soluciones en Python.

Soluciones en Haskell

module Limites_de_sucesiones where

import Test.Hspec (Spec, hspec, it, shouldBe)

limite :: (Double -> Double) -> Double -> Double
limite f a =
  head [f x | x <- [1..],
              maximum [abs (f y - f x) | y <- [x+1..x+100]] < a]

-- Verificación
-- ============

verifica :: IO ()
verifica = hspec spec

spec :: Spec
spec = do
  it "e1" $
    limite (\n -> (2*n+1)/(n+5)) 0.001  `shouldBe`  1.9900110987791344
  it "e2" $
    limite (\n -> (1+1/n)**n) 0.001     `shouldBe`  2.714072874546881

-- La verificación es
--    λ> verifica
--
--    e1
--    e2
--
--    Finished in 0.1927 seconds
--    2 examples, 0 failures

module Limites_de_sucesiones where

import Test.Hspec (Spec, hspec, it, shouldBe)

limite :: (Double -> Double) -> Double -> Double

limite f a =

head [f x | x <- [1..],

maximum [abs (f y - f x) | y <- [x+1..x+100]] < a]

-- Verificación

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

verifica :: IO ()

verifica = hspec spec

spec :: Spec

spec = do

it "e1" $

limite (\n -> (2*n+1)/(n+5)) 0.001 `shouldBe` 1.9900110987791344

it "e2" $

limite (\n -> (1+1/n)**n) 0.001 `shouldBe` 2.714072874546881

-- La verificación es

-- λ> verifica

-- e1

-- e2

-- Finished in 0.1927 seconds

-- 2 examples, 0 failures

Soluciones en Python

from itertools import count
from typing import Callable

# 1ª solución
# ===========

def limite(f: Callable[[float], float], a: float) -> float:
    x = 1
    while True:
        maximum_diff = max(abs(f(y) - f(x)) for y in range(x+1, x+101))
        if maximum_diff < a:
            return f(x)
        x += 1

# 2ª solución
# ===========

def limite2(f: Callable[[float], float], a: float) -> float:
    x = 1
    while True:
        y = f(x)
        if max(abs(y - f(x + i)) for i in range(1, 101)) < a:
            break
        x += 1
    return y

# 3ª solución
# ===========

def limite3(f: Callable[[float], float], a: float) -> float:
    for x in count(1):
        if max(abs(f(y) - f(x)) for y in range(x + 1, x + 101)) < a:
            r = f(x)
            break
    return r

# Verificación
# ============

def test_limite() -> None:
    assert limite(lambda n :  (2*n+1)/(n+5), 0.001) ==  1.9900110987791344
    assert limite(lambda n : (1+1/n)**n, 0.001)     ==  2.714072874546881
    assert limite2(lambda n :  (2*n+1)/(n+5), 0.001) ==  1.9900110987791344
    assert limite2(lambda n : (1+1/n)**n, 0.001)     ==  2.714072874546881
    assert limite3(lambda n :  (2*n+1)/(n+5), 0.001) ==  1.9900110987791344
    assert limite3(lambda n : (1+1/n)**n, 0.001)     ==  2.714072874546881
    print("Verificado")

# La comprobación es
#    >>> test_limite()
#    Verificado

from itertools import count

from typing import Callable

# 1ª solución

# ===========

def limite(f: Callable[[float], float], a: float) -> float:

x = 1

while True:

maximum_diff = max(abs(f(y) - f(x)) for y in range(x+1, x+101))

if maximum_diff < a:

return f(x)

x += 1

# 2ª solución

# ===========

def limite2(f: Callable[[float], float], a: float) -> float:

x = 1

while True:

y = f(x)

if max(abs(y - f(x + i)) for i in range(1, 101)) < a:

break

x += 1

return y

# 3ª solución

# ===========

def limite3(f: Callable[[float], float], a: float) -> float:

for x in count(1):

if max(abs(f(y) - f(x)) for y in range(x + 1, x + 101)) < a:

r = f(x)

break

return r

# Verificación

# ============

def test_limite() -> None:

assert limite(lambda n : (2*n+1)/(n+5), 0.001) == 1.9900110987791344

assert limite(lambda n : (1+1/n)**n, 0.001) == 2.714072874546881

assert limite2(lambda n : (2*n+1)/(n+5), 0.001) == 1.9900110987791344

assert limite2(lambda n : (1+1/n)**n, 0.001) == 2.714072874546881

assert limite3(lambda n : (2*n+1)/(n+5), 0.001) == 1.9900110987791344

assert limite3(lambda n : (1+1/n)**n, 0.001) == 2.714072874546881

print("Verificado")

# La comprobación es

# >>> test_limite()

# Verificado

Algoritmo de bajada para resolver un sistema triangular inferior

Integración por el método de los rectángulos

Raíces enteras

Método de bisección para calcular raíces de una función

Límites de sucesiones

Entradas recientes

Buscar

Correo electrónico