• Contenido Creative Commons Attribution license, CC-BY
  • Código con licencia MIT

creative commons attribution

  • (c) Kyle T. Mandli
  • (c) Modificado por Pablo Angulo para ETSIN@UPM

Método numérico

Un método numérico es un procedimiento para aproximar la solución de un problema matemático:

  • No todos los problemas admiten soluciones exactas.
  • El cálculo simbólico tampoco puede resolver todos los problemas.
  • Incluso aunque la solución exista, puede ser demasiado costoso encontrar esa solución.
  • Muchos métodos de cálculo exacto o simbólico no tienen garantías de terminar en tiempo razonable.

En general, un método numérico debe ser capaz de aproximar la solución al problema matemático con tanta precisión como queramos. En la práctica, tanta precisión es innecesaria, y muy costosa. Todos los métodos numéricos deben aproximar la solución buscando un equilibrio entre:

  • poco error
  • poco tiempo de cómputo

Tipos de error

  • Errores en los datos: error de medición, ruido.
  • Errores al simplificar la realidad mediante un modelo matemático: error de modelo.
  • Errores al reemplazar el objeto matemático de interés (número, vector, función, etc) por otro objeto matemático más fácil de calcular: error de truncamiento.
  • Errores al representar números reales con una cantidad finita de cifras decimales, y al hacer operaciones con números de coma flotante en vez de con números reales: error de redondeo.

En esta asignatura, hablaremos únicamente de los dos últimos tipos de error, pero es importante recordar que hay muchas fuentes de error:

  • No tiene sentido buscar una solución del problema matemático con precisión mayor que la permiten los errores de medición o de modelo.
  • Es interesante diseñar métodos robustos a errores.

Salvo en ejemplos de laboratorio, nunca podremos conocer exactamente el error, pero a veces podremos acotarlo.

Error de modelo

A menudo son simplificaciones conscientes, cuyo efecto en la solución se puede estimar

  • Ignoramos fuerzas que sabemos que actúan pero pensamos que tienen poca magnitud
  • Aproximamos el número de barcos por un número real
  • Simplificamos la geometría del buque
  • Asumimos que no hay mezcla entre capas de fluído
  • un largo etcétera

Pero también a veces el error de modelo es muy difícil de estimar

  • Efecto del metano en modelos climáticos
  • Burbujas y crisis son influenciados por percepciones humanas
  • otro largo etcétera

Error de truncamiento

Errores al sustituir un objeto matemático por otro más sencillo:

  • aproxima una función por otra: $sin(x) \approx x$ for $|x| \approx 0$.
  • aproxima un número por otro: $x$, solución de $10x + \sin(x)=1$ se puede aproximar por la solución de $10x=1$
  • veremos otros muchos ejemplos

Error absoluto y relativo

Si $x$ es el número que buscamos, $y$ es la aproximación:

  • Error absoluto: $$ e = |x-y|$$
  • Error relativo: $$ r = \frac{e}{|x|} = \frac{|x - y|}{|x|} $$ También se usa a menudo el error cuadrático, o error al cuadrado: $$ e = (x-y)^2$$ aunque a priori es una deformación del error real, y no tiene las unidades correctas, a menudo es más fácil de trabajar con este tipo de error.

Atención a las unidades de cada tipo de error.

Error de truncamiento en el Teorema de Taylor

Teorema de Taylor: Sea $f(x) \in C^{m+1}[a,b]$, $x_0 \in [a,b]$.

Para cualquier $x \in (a,b)$ existe un número $c = c(x)$ en el intervalo entre $x_0$ y $x$ tal que

$$ f(x) = T_N(x) + R_N(x)$$

donde $T_N(x)$ es el polinomio de Taylor de $f$ de grado $N$ en $x0$:

$$T_N(x) = \sum^N_{n=0} \frac{f^{(n)}(x_0)\cdot(x-x_0)^n}{n!}$$

y $R_N(x)$ es el error de Taylor

$$R_N(x) = \frac{f^{(n+1)}(c) \cdot (x - x_0)^{n+1}}{(n+1)!}$$

Es un error de truncamiento, independiente de cuántos dígitos de precisión usemos para representar los números reales.

El error de truncamiento se puede acotar si podemos acotar la derivada $f^{(n+1)}(c)$ por una constante $M$:

$$ R_N(x) = \frac{f^{(n+1)}(c) \cdot (x-x_0)^{n+1}}{(n+1)!} \leq M (x-x_0)^{n+1} $$

Ejemplo

$f(x) = e^x$, $x_0 = 0$

Podemos acotar el error que comete la aproximación $T_2$, en función de $x$.

Derivadas: $$\begin{aligned} f'(x) &= e^x \\ f''(x) &= e^x \\ f^{(n)}(x) &= e^x \end{aligned}$$

Polinomios de Taylor: $$\begin{aligned} T_N(x) &= \sum^N_{n=0} e^0 \frac{x^n}{n!} \Rightarrow \\ T_2(x) &= 1 + x + \frac{x^2}{2} \end{aligned}$$

Restos (usando $e<3$): $$\begin{aligned} R_N(x) &= e^c \frac{x^{n+1}}{(n+1)!}\Rightarrow \\ R_2(x) &= e^c \cdot \frac{x^3}{6} \leq \frac{e^x}{6} \end{aligned}$$

Es decir, aproximamos $e^1 = 2.718$ por $T_2(1) = 2.5$, y sabemos que comete un error menor que $\frac{e^1}{6}<0.5$. En realidad, el error es menor.

Números de coma flotante

Un sistema de números de coma flotante, representa algunos números reales en la forma siguiente:

$$F = \pm d_1 . d_2 d_3 d_4 \ldots d_p \times \beta^E$$

donde

  1. $\pm$ es un bit que representa el signo del número
  2. $d_1 . d_2 d_3 d_4 \ldots d_p$ se denomina la mantisa. El punto decimal siempre está detrás del primer dígito. Se asume que $d_1 \neq 0$ a menos que el número esa exactamente $0$.
  3. $\beta$ es la base. Puede ser base $\beta = 2$ para el sistema binario, o $\beta = 10$ para el decimal
  4. $E$ es el exponente, un número entero.

Características importantes

  1. La cantidad de números representables es finita
  2. Estos números no están distribuidos de forma uniforme
  3. La aritmética con números de coma flotante no es exacta, porque el resultado de una operación con dos números representables podría no ser representable (en es caso, se busca el número representable más cercano).

El desbordamiento ocurre cuando el resultado de una operación es:

  • overflow: mayor que el número más grande representable en el sistema
  • underflow: menor que el número positivo más pequeño representable en el sistema Un overflow da lugar a un número especial del sistema: $\inf$ ó $-\inf$. Las operaciones continúan, y por ejemplo $1/\inf=0$, pero si se pide una indeterminación, el resultado es otro "número" especial del sistema: undefined.

Real Systems - IEEE 754 Binary Floating Point Systems

Single Precision

  • Total storage alloted is 32 bits
  • Exponent is 8 bits $\Rightarrow E \in [-126, 127]$
  • Fraction 23 bits ($p = 24$)
s EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF
0 1      8 9                     31

Overflow $= 2^{127} \approx 3.4 \times 10^{38}$

Underflow $= 2^{-126} \approx 1.2 \times 10^{-38}$

$\epsilon_{\text{machine}} = 2^{-23} \approx 1.2 \times 10^{-7}$

Double Precision

  • Total storage alloted is 64 bits
  • Exponent is 11 bits $\Rightarrow E \in [-1022, 1024]$
  • Fraction 52 bits ($p = 53$)
s EEEEEEEEEE FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FF
0 1       11 12                                                      63

Overflow $= 2^{1024} \approx 1.8 \times 10^{308}$

Underflow $= 2^{-1022} \approx 2.2 \times 10^{-308}$

$\epsilon_{\text{machine}} = 2^{-52} \approx 2.2 \times 10^{-16}$

Cuando queramos distinguir las operaciones con números reales de las operaciones en el sistema de coma flotante, usaremos $\ominus$ para la resta, $ \oplus$ para la suma, etc.

Se pueden definir de esta forma:

  • $fl(x)$ es el número de coma flotante más próximo a $x$ (hay varias formas de redondear, pero no las vamos a discutir).
  • $x\oplus y = fl(x + y)$
  • $x\otimes y = fl(x \times y)$
  • etcétera...

Si los operandos $x$ e $y$ no son representables en el sistema, hay que definir:

  • $x\oplus y = fl(fl(x) + fl(y))$
  • $x\otimes y = fl(fl(x) \times fl(y))$
  • etcétera...

Ejemplo 1:

Elegimos un $\delta$ más pequeño que el $ \epsilon_{\text{machine}}$ (el número positivo más pequeño del sistema):

$$1\oplus \delta \ominus 1 = (1\oplus \delta) \ominus 1 = 0$$

$$1 \ominus 1 \oplus \delta = (1 \ominus 1) \oplus \delta = \delta$$

In [2]:
eps

ans =    2.2204e-16
In [4]:
eps
1 + eps - 1 
1 - 1  + eps
'Ahora con un número menor que el epsilon de la máquina'
eps/2
1 + eps/2 - 1 
1 - 1  + eps/2

ans =    2.2204e-16
ans =    2.2204e-16
ans =    2.2204e-16
ans = Ahora con un número menor que el epsilon de la máquina
ans =    1.1102e-16
ans = 0
ans =    1.1102e-16

Cancelación

Si el ejemplo anterior era "inofensivo", observa qué ocurre cuando restamos números que son casi iguales.

Supóngase $x$ e $y$ dos números casi iguales con $x>y$, con sus representaciones en coma flotante de $k$ dígitos $$ \begin{align} fl(x)=0.d_1d_2\ldots d_p\alpha_{p+1}\alpha_{p+2}\ldots \alpha_{k}\times10^n \nonumber \\ fl(y)=0.d_1d_2\ldots d_p\beta_{p+1}\beta_{p+2}\ldots \beta_{k}\times10^n \nonumber \end{align} $$ Al restarlos $x\ominus y$: $$ fl(fl(x)-fl(y))=0.\sigma_{p+1}\sigma_{p+2}\ldots\sigma_{k}\times10^{n-p} $$ con $$ 0.\sigma_{p+1}\sigma_{p+2}\ldots\sigma_{k}=0.\alpha_{p+1}\alpha_{p+2}\ldots\alpha_{k}-0.\beta_{p+1}\beta_{p+2}\ldots\beta_{k} $$ La resta sólo tiene $k-p$ dígitos de precisión, aunque los operandos tenían $k$ dígitos de precisión.

Ejemplo: Función racional

Calculamos $f(x) = x + 1$ mediante $$F(x) = \frac{x^2 - 1}{x - 1}$$, y dibujamos el error $|f-F|$.

In [1]:
x = linspace(0.9, 1.1, 1001);
f_hat = (x.^2 - 1.0) ./ (x - 1.0);

plot(x, abs(f_hat - (x + 1.0)));
xlabel("$x$")
ylabel("Absolute Error")

Ejemplo 4: Evaluar un polinomio

Veremos que varias técnicas numéricas con mucha tradición se basan en polinomios de grado alto.

Pero evaluar un polinomio de grado alto introduce errores numéricos:

$$f(x) = x^7 - 7x^6 + 21 x^5 - 35 x^4 + 35x^3-21x^2 + 7x - 1$$

In [2]:
x = linspace(0.988, 1.012, 1000);
y = x.^7 - 7.0 * x.^6 + 21.0 * x.^5 - 35.0 * x.^4 + 35.0 * x.^3 - 21.0 * x.^2 + 7.0 * x - 1.0;

plot(x, y, 'r')
xlabel("x")
ylabel("y")

En este caso, el ruido numérico disminuye si evaluamos el polinomio de otra forma.

Esta evaluación anidada requiere algunas operaciones menos y acumula menos error que la evaluación naive de un polinomio.

El orden al hacer las operaciones importa.

In [3]:
x = linspace(0.988, 1.012, 1000);
y = ((((((x - 7.0).*x + 21.0).*x - 35.0).*x + 35.0).*x - 21.0).*x + 7.0).* x - 1.0;

plot(x, y, 'r')
xlabel("x")
ylabel("y")

Combinación de los dos tipos de error

En general tendremos que lidiar con los dos tipos de error a la vez. Con suerte, un error predomina, lo que simplifica el análisis, pero en general, los errores se propagan, se acumulan y a veces se amplifican.

Ejemplo

Aproximamos la derivada de $f(x) = e^x$ en $x=1$ por la diferencia:

$$f'(x) \approx \frac{f(x + \Delta x) - f(x)}{\Delta x}$$

Comparamos el error cuando $\Delta x$ decrece, entre la aproximación y la respuesta correcta $f'(1) = e$.

In [31]:
n = 60;
delta_x = zeros(1,n);
for i = 1:n
    delta_x(i) = 2^(-i);
end
x = 1.0;
f_hat = (exp(x + delta_x) - exp(x)) ./ (delta_x);

loglog(delta_x, abs(f_hat - exp(1)), 'o-')
xlabel("Delta x")
ylabel("Absolute Error")
In [ ]:
In [ ]: