Introducción al Cálculo numérico. Errores

Del algoritmo al programa

Una vez  escrito un algoritmo para resolver un problema es inmediato realizar el programa correspondiente en algún lenguaje. Mathematica es muy versátil porque tiene muchas órdenes que permiten resolver directamente un problema o parte del mismo y, además, permite estructurar el programa en pequeños módulos que podemos procesar por separado para verificar si funcionan correctamente y unirlos después, si se desea.

Vamos a calcular la suma de n números. Para ello los introducimos en un vector v y después sumamos las componentes de v mediante el algoritmo correspondiente.

\(v = {1, \(-3.6\), 4.2, 5}; 
n = Length[v]; 
s = 0; 

    For[i = 1, i ≤ n, \(i++\), \ s = s + v[\([i]\)]]; s\)

Podemos obtener el mismo resultado con la orden directa Sum.

Sum[v[[i]], {i, n}]

A partir de la versión 3.0 de Mathematica la orden Sum[v[[i]],{i,imin,imax}] puede expresarse de forma más tradicional (∑_(i = imin)^imax) v[[i]] o \(TraditionalForm\`∑\+\(i = imin\)\%imax v\_\(\([[\)\(i\)\(]]\)\)\), seleccionándose el símbolo de la suma en la paleta denominada BasicInput que se encuentra en el submenú Palettes del menú File. Utilizándolo, el cálculo precedente se puede efectuar de la siguiente manera:

(∑_(i = 1)^n) v[[i]]

Programamos a continuación el algoritmo permite calcular el producto de las componentes no nulas de un vector v.

\(v = {1, 0, 1.2, \(-1\), 7}; 
n = Length[v]; p = 1; 

    For[i = 1, i ≤ n, \(i++\), 
      If[v[\([i]\)] ≠ 0, p = p\ v[\([i]\)]]]; p\)

Comportamientos extraños debidos a los errores de redondeo

El ordenador puede trabajar con números que sean muy pequeños o muy grandes, pero, si trabaja con ambos simultáneamente, el número pequeño puede perder su valor debido al redondeo. Por ejemplo, para un valor grande de x es posible que 51.2 + x lo tome como x. Para simplificar, nos vamos a limitar a la sucesión de números 10, 100, 1000, etc., es decir, de la forma (10^i). Podemos buscar el menor de ellos que cumple 51.2 + x = x con el siguiente programa:

For[i = 1, i ≤ 100, i++,           
      
       
      If[51.2 + (10^i) == (10^i),  
        
            
        Print["El exponente de 10 con el que se consigue 
la igualdad es i = ", i]; 
            
        Print["Valor encontrado: x = ", 10^i]; 
        
            Break[], 
        
            
        Print["Para x = ", (10^i), " no se da la igualdad porque ", 
          51.2 + (10^i), " es diferente de ", (10^i)]; 
        
            Print[]]]; x = (10^i);

¿Con qué precisión trabaja nuestro ordenador?

No muestra más que 6 cifras, aunque ha necesitado llegar a (10^16) para que 51.2 deje de influir en la suma; después comentaremos lo de las cifras de salida. Podemos comprobar que para valores de x mayores o iguales al obtenido anteriormente se tiene (51.2 + x) + (-x) = 0, mientras que (51.2) + (x + (-x)) = 51.2, perdiéndose la asociatividad en la suma.

Print[(51.2 + x) + (-x), " diferente de ", 
      51.2 + (x + (-x))]

Print[(51.2 + 10  x) + (-10 x), " diferente de " , 
      51.2 + (10  x + (-10 x))]

Print[(51.2 + 100  x) + (-100 x), " diferente de ", 
      51.2 + (100  x + (-100 x))]

Si aproximamos la derivada de una función en el punto x = 1 por (1/h) (f(x + h) - f(x)), cuando h tome un valor similar al obtenido anteriormente, la respuesta va a ser cero, pero antes habrán entrado en juego los errores por truncatura y, después, por cancelación (éste producirá algún resultado extraño).

Error por truncatura, cancelación y por redondeo

Supongamos que el valor obtenido de i anteriormente fue menor que 20. Vamos a escribir aproximaciones de la derivada de (e^x) en el punto x = 1. Las almacenamos en una lista, aprox, que posteriormente mostramos.

La derivada de  (e^x) en el punto x = 1 es  e.

\(aprox = {}; 

    For[i = 1, i ≤ 20, \(i++\), 
      AppendTo[aprox, \(Exp[1. + 10\^\(-i\)] - Exp[1.]\)\/10\^\(-i\)]]\)aprox

Las primeras aproximaciones distan sustancialmente del valor exacto debido a que el error de la fórmula empleada es importante cuando h no es muy pequeño, ya que procede de cortar el desarrollo de Taylor por el tercer término y el error queda del mismo orden de h (error por truncatura). Cuando aumenta el valor de i la aproximación de la derivada cada vez se acerca más al valor 2.718281. Sin embargo, al final da cero, como era de esperar, puesto que resta dos números idénticos para el ordenador (redondea x + h y x al mismo valor); pero antes de dar cero hay varios resultados que empiezan a alejarse del valor exacto debido a que se está produciendo un importante error por cancelación.

Salvo cuando se opera en simbólico, el redondeo es inevitable en casi todas las operaciones, incluso en la entrada de datos. Para comprobar el efecto en la entrada de datos vamos a tratar de imprimir todos los números, a partir del 5, de 0.1 en 0.1, que son menores que 6 en un caso y menores o iguales a 6 en el otro. Si su ordenador no opera en decimal, sino en binario -que es lo usual-, dará la misma respuesta en ambos casos, por lo que uno de los resultados es erróneo.

\(lista = {}; 

    For[i = 0, i ≤ 1, i = i + 0.1, AppendTo[lista, 5 + i]]\)lista

\(lista = {}; 
For[i = 0, i < 1, i = i + 0.1, AppendTo[lista, 5 + i]]\)lista

En este caso, el error se podría haber evitado poniendo i = i + (1/10), ya que realizaría la operación en simbólico y por tanto no cometería error. Compruébelo haciendo el cambio correspondiente en la celdilla que ha fallado.

Nuevamente comprobamos el error de cancelación al procesar el valor de la diferencia entre dos cantidades muy similares, concretamente entre las raíces de y + 1 e y multiplicada por y, ((Sqrt[y + 1]) - (Sqrt[y])) y, para valores grandes de y, que equivale a (y/((Sqrt[y + 1]) + (Sqrt[y]))) y que vale aproximadamente ((Sqrt[y])/2). Tomamos y = 2.5 (10^i). Almacenamos los valores obtenidos en una lista y la mostramos en forma de tabla con cabeceras.

lista = {};
For[i = 1,i ≤ 20,i++,
pordiferencia=((Sqrt[1+2.5  (10^i)])-(Sqrt[2.5  (10^i)]))2.5  (10^i);
porcociente=((2.5  (10^i))/((Sqrt[1+2.5  (10^i)])+(Sqrt[2.5  (10^i)])));
AppendTo[lista, {pordiferencia,  porcociente}]]TableForm[lista, 
      TableHeadings 
        -> {None, {"Por diferencia", "Por cociente"}}]

Condicionamiento

Hay matrices que están bien condicionadas, por ejemplo las de Householder. Estas matrices son de la forma I - 2 u (u^T), donde u es un vector con norma euclídea igual a la unidad. Por ejemplo, sea u = ((4/7), (4/7), (4/7), (1/7)). Entonces la matriz de Householder correspondiente es

\(u = {4/7, 4/7, 4/7, 1/7}; 

    A = IdentityMatrix[4] - 
        2  Table[u[\([i]\)]*u[\([j]\)], {i, 4}, {j, 4}]\)

Veamos cómo cambia la solución de un sistema cuando se perturba su matriz de coeficientes, A. Por ejemplo, vamos a tomar como términos independientes (por comodidad) el vector \(TraditionalForm
      \`\((\(-\(55\/49\)\), \(-\(55\/49\)\), \(-\(55\/49\)\), 23\/49)\)\), con lo cual la solución del sistema correspondiente es (1, 1, 1, 1).

b = {-55/49, -55/49, -55/49, 23/49};

LinearSolve[A, b]

Perturbamos cada coeficiente de forma aleatoria, por ejemplo, con valores comprendidos entre -0.05 y 0.05, para lo que definimos, y mostramos, una matriz cuyos elementos son dichas perturbaciones.

B = Table[(1/10) Random[]  - 0.05, {i, 4}, {j, 4}]LinearSolve[A + B, b]

Podemos procesar varias veces y ver que las soluciones oscilan poco con respecto a la solución del sistema original. Pero si repetimos lo mismo con la matriz A que damos a continuación, los resultados cambian significativamente:

A = {{93, -80, 88, 25}, {-69, 88, 87, 35}, {7, 7, 6, 6}, {85, 
          -84, 83, 20}};

Para tener, de nuevo, por simplicidad, la solución (1, 1, 1, 1), tomamos como b el vector (126, 141, 26, 104).

\(b = {126, 141, 26, 104}; 
LinearSolve[A, b]\)

B = Table[(1/10) Random[]  - 0.05, {i, 4}, {j, 4}]LinearSolve[A + B, b]

Las matrices de Hilbert son otro ejemplo de mal condicionamiento. Definamos la de orden 4 y un vector de términos independientes.

\(A = Table[1\/\(i + j\), {i, 4}, {j, 4}]; 

    b = Table[Sum[A[\([i, j]\)], {j, 4}], {i, 4}]\)

La solución del sistema correspondiente es

LinearSolve[A, b]

Perturbamos la matriz A y resolvemos el nuevo sistema.

B = Table[(1/10) Random[]  - 0.05, {i, 4}, {j, 4}]LinearSolve[A + B, b]

Veamos ahora cómo se comporta la evaluación de diferentes funciones, por ejemplo suponemos dos funciones (Subscript[f, 1]) y (Subscript[f, 2]), que se van a evaluar en x = 0.05. Vamos a mostrar también los valores de  (Subscript[f, 1]) y (Subscript[f, 2]) en otros dos puntos próximos a 0.05, que distan de éste una milésima. y en dos puntos más: uno fijo.

(Subscript[f, 1])[x_] := Exp[200  x] + 20(Subscript[f, 2])[x_] := Exp[-200 x] + 20

{(Subscript[f, 1])[0.049], (Subscript[f, 1])[0.05], (Subscript[f, 1])[0.051]}

{(Subscript[f, 2])[0.049], (Subscript[f, 2])[0.05], (Subscript[f, 2])[0.051]}

En el primer caso las diferencias son enormes, mientras que en el segundo son insignificantes. Una variación de una milésima en el valor de x implica una variación igual a \(\(1\/1000\) \(f\_1\%'\) \((0.05)\) = 4405.29\) en la primera función y \(\(1\/1000\) \(f\_2\%'\) \((0.05)\) = 0.000009\).

Estabilidad

Vamos a generar 20 términos de la sucesión (1/(2^n)) por dos procedimientos, considrando una precisión de 6 dígitos. Uno estable y otro inestable cuando operamos con 6 cifras. Inicialmente no cometeremos error alguno.

\(x = 1\/2.; 
lista = {}; 

    For[i = 1, i ≤ 20, \(i++\), 
      
AppendTo[lista, x]; 

      x = N[Round[\(10\^6\) x\/2]\/10\^6]]\)lista

\(x = 1\/2; y = 1\/4.; 
lista = {}; 

    For[i = 1, i ≤ 30, \(i++\), 
\ \ \ \ \ \ \ 
      z = 23\/2\ y\  - \(11\/2\) x; \ 
\ \ \ \ \ \ 
      AppendTo[lista, z]; 
\ \ \ \ \ \ \ x = y; 
\ 
      y = N[Round[\(10\^6\) z]\/10\^6]]\)lista

Precisión con Mathematica

Cuando los números que se usan como datos son todos enteros, y no se indica que se utilice una aproximación numérica de los mismos, Mathematica opera en modo simbólico y, por tanto, con precisión infinita. Si alguno de ellos contiene cifras decimales lo considera aproximado y todas las operaciones en las que intervenga las hará de forma aproximada. Observe los siguientes ejemplos:

(2/5) + (32/10) - Log[2] - (Sqrt[3]) + π - E

(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E

Pero, ¿con qué precisión ha considerado el valor 3.2? Podemos pedir el resultado anterior con más cifras decimales. Veamos qué ocurre cuando le pedimos 10, 20 o 50 cifras.

N[(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E, 10]

N[(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E, 20]

N[(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E, 50]

Podemos conseguir más de 6 cifras con la orden SetPrecision.

SetPrecision[(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E, 8]

Incluso proporciona muchas más:

SetPrecision[(2/5) + 3.2 - Log[2] - (Sqrt[3]) + π - E, 52]

Pero, ¿son verdaderas?

SetPrecision[(2/5) + (32/10) - Log[2] - (Sqrt[3]) + π - E, 52]

Observe que, en este caso, da una cantidad algo menor que al usar el valor decimal en 3.2. Sólo coinciden las 16 primeras al usar 3.2 y (32/10).  Ello se debe a que cuando se le da un número con punto decimal, por ejemplo 0.7 o 3.2,  Mathematica lo considera con 16 cifras, en este caso 0.7000000000000000  o 3.2000000000000000, y lo redondea al número máquina más próximo a él (unas veces el redondeo será por defecto y otras por exceso). Veámoslo.

SetPrecision[0.7, 55]

SetPrecision[(7/10), 55]

SetPrecision[3.2, 55]

Otros ejemplos son los siguientes:

SetPrecision[123456.12, 20]

N[1234.56789012345678901234567890123456, 50]

SetPrecision[1234.56789012345678901234567890123456, 50]

\(a = N[7\/10, 100]; 
SetPrecision[a, 55]\)

SetPrecision[a, 115]

Así pues, cuando demos un número con cifras decimales, Mathematica considera que no tiene más de 16 cifras y lo aproxima por el número máquina más próximo a él. Además, los resultados de operar con ese número los muestra con 6 cifras. Sin embargo, con la orden SetPrecision se pueden rescatar más cifras, pero podemos caer en error al pedirle más cifras que las dadas en la entrada de los datos.

Hay que cuidar bien la entrada de datos. Observemos el resultado de la siguiente operación cuyo valor exacto es cero:

3.1  (10^50) - 31  (10^49)

Ejercicios

1.- Programe el paso de la numeración decimal a la binaria y aplíquelo para obtener la representación en binario de los números del sistema decimal 5500 y 0.1.

2.- Obtenga aproximaciones de la derivada de la función seno en el punto 0.25, a partir de la fórmula (1/h) (sen(0.25 + h) - sen(h)), con h de la forma (10^(-i)), con i = 1, 2, 3,  ... , deteniéndose cuando la diferencia entre dos aproximaciones consecutivas aumente.

3.- Programe la sucesión de término general (1/(3^n)), a partir de las siguientes leyes de recurrencia:
	a) (Subscript[a, 0]) = 1; (Subscript[a, n]) = (1/3) (Subscript[a, n - 1]),  
      n ≥ 1.
	b) (Subscript[a, 0]) = 1; (Subscript[a, 1]) = (1/3); (Subscript[a, n]) = (4/3) (Subscript[a, n - 1]) - (1/3) 
            (Subscript[a, n - 2]),  n ≥ 2.
	c) (Subscript[a, 0]) = 1, (Subscript[a, 1]) = (1/3), (Subscript[a, n]) = (10/3) (Subscript[a, n - 1]) - (Subscript[a, n - 2]).
Comprobará la mala propagación de los errores en algunos casos.