Informática de Sistemas
Cálculo Numérico con Mathematica

Tema 1
Análisis de errores (3 ejercicios)

Nota: Se escribe al principio de la práctica la instrucción Clear["Global`*"], que permite borrar todos los símbolos almacenados (en concreto, valores de variables o funciones) desde que se inició la sesión con Mathematica. Así se facilita su uso siempre que sea necesario.

Clear["Global`*"]

La aritmética aproximada en Mathematica

Como se vio en la práctica introductoria, si se utilizan números con punto decimal, entonces Mathematica los redondeará si no tiene cifras decimales suficientes para mostrarlos:

{(4/7) + (1/2), ((Sqrt[8])/89672)}              

    {(4. /7) + (1/2),  ((Sqrt[8.]) /89672)}

Los otros dígitos que guarda Mathematica se pueden visualizar, como vimos, utilizando la instrucción InputForm (también vimos otras dos formas de visualizarlos)

{InputForm[(4. /7) + (1/2)], InputForm[((Sqrt[8.]) /89672)]}

También vimos que para obtener una aproximación numérica a partir de una exacta se puede utilizar la intrucción N. Viene a lo mismo que escribir la expresión con algún punto decimal:

{N[((Sqrt[8])/89672)], InputForm[N[((Sqrt[8])/89672)]], 
      InputForm[((Sqrt[8.]) /89672)]}

Cuando se desea que Mathematica proporcione una precisión mayor, esta debe especificarse como segundo argumento en la instrucción N.

N[((Sqrt[8])/89672), 20]

Parece que se nos ha dado el número (Sqrt[8])/89672 con 20 dígitos de precisión. Pero en realidad en este caso tenemos más precisión, como muestra la siguiente instrucción:

InputForm[N[((Sqrt[8])/89672), 20]]

Comprobemos cuántos de los dígitos anteriores son correctos:

N[((Sqrt[8])/89672), 55]

Es decir, todos menos el último (54 dígitos) eran correctos.
Otros ejemplos de uso de la instrucción N:

N[π, 500]

InputForm[N[π, 500]]

N[2^(-1025)]
N[1 - (2^(-53)), 60]
N[(4/7) + (1/2), 33]

No tiene sentido aplicar la orden N sobre un número ya redondeado. Si se hace no modifica nada el número de partida:

{InputForm[(4. /7) + (1/2)], 
InputForm[
        N[(4. /7) + (1/2)]], 
InputForm[
        N[(4. /7) + (1/2), 30]], 
InputForm[
        N[(4. /7) + (1/2), 8]]}

Al realizar operaciones, los redondeos pueden acumularse o cancelarse. Por ejemplo, la siguiente expresión se ha calculado con un máximo de 5 dígitos correctos en la mantisa decimal  (luego se verá que son 4):

N[12345678901234567890]*N[987654321987654321 (10^(-22))] - 
      N[1219326312460000]
InputForm[%]

En este caso se puede evitar el problema haciendo primero el cálculo exacto de la expresión y luego lo normal sería usar la instrucción N, aunque en este caso excepcional no hace falta:

12345678901234567890*987654321987654321 (10^(-22)) - 
      1219326312460000
N[%]
InputForm[%]

Como se ha visto, si se pueden hacer los cálculos de modo exacto (y solo al final pedir un valor aproximado), mucho mejor. Pero esto pocas veces será posible con problemas de Cálculo Numérico.

Números máquina

Un número máquina para un procesador determinado es cualquier número que admita una representación exacta en su aritmética de punto flotante.

Ejemplo

Muestre todos los números máquina de un procesador decimal de dos cifras de mantisa y exponente entre -3 y 3.

Solución

Una primera forma de resolver el problema es utilizando las instrucciones For y Print:

For[j = 0.1, j ≤ 0.99, j = j + 0.01, 
      For[i = -3, i ≤ 3, i++, 
        Print[j (10^i), "     ", -j (10^i)]]]

Pero quizá lo más sencillo es utilizar la instrucción Table.

lista = 
      Table[(((-1))^k) j (10^i), {j, 0.1, 0.99, 0.01}, {i, -3, 
          3}, {k, 2}]; TableForm[lista]

Ejemplo

Utilice Mathematica para encontrar el épsilon de la máquina donde trabaja, aplicando que el épsilon es un número de la forma (2^(-n)), con n ∈ ℕ. Compare el épsilon de la máquina con los menores números de la forma (2^(-n)) tales que en Mathematica resulte que (a) (10^(-6)) + (2^(-n)) sume más que (10^(-6)) y (b) (10^6) + (2^(-n)) sume más que (10^6).

Solución

Se puede empezar por la siguiente instrucción, y luego ajustar hasta llegar, por ejemplo, a la escrita inmediatamente después.

Table[InputForm[1.  + (2^(-n))], {n, 100}]

Table[InputForm[1.  + (2^(-n))], {n, 50, 55}]

Luego el épsilon de la máquina es (2^(-52)) (el doble del épsilon del ordenador típico considerado en la teoría de la asignatura). Su valor es aproximadamente igual a:

(2. ^(-52))

(a) Sustituyendo 1 por (10^(-6)):

Table[InputForm[(10. ^(-6)) + (2^(-n))], {n, 70, 75}]

Luego el menor número que suma algo a (10^(-6)) es (2^(-73)):

(2. ^(-73))

(b) Sustituyendo 1 por (10^6):

Table[InputForm[(10^6.)  + (2^(-n))], {n, 30, 35}]

Luego el menor número que suma algo a (10^6) es (2^(-33)):

(2. ^(-33))

Pérdida de dígitos significativos por cancelación

Uno de los ejemplos más comunes de la propagación de errores es el fenómeno conocido como pérdida o cancelación de cifras significativas. Esta situación se da, por ejemplo, cuando se opera con el valor de la diferencia entre dos cantidades muy similares.

Ejemplo

Considere el cálculo de valores de la función f(x)=(Sqrt[x + 1])-(Sqrt[x]) para x "grande". Si x=(10^i), para 3 ≤ i ≤ 20, las siguientes instrucciones elaboran una tabla los distintos valores de f((10^i)), con 3 ≤ i ≤ 20.

Solución

TableForm[Table[{i, (Sqrt[1.  + (10^i)]) - (Sqrt[(10^i)])}, {i, 3, 20}], 
      TableHeadings -> {None, {"Exponente base 10", "Por \
diferencia"}}]

O con la instrucción For:

lista = {}; 
For[i = 3, i ≤ 20, i++, 
      pordiferencia = (Sqrt[1.  + (10^i)]) - (Sqrt[(10^i)]); 
      AppendTo[lista, {i, pordiferencia}]]
TableForm[lista, 
      TableHeadings -> {None, {"Exponente base 10", "Por \
diferencia"}}]

La cancelación puede evitarse si en la expresión f(x)=(Sqrt[x + 1])-(Sqrt[x]) se multiplica y se divide por el conjugado, obteniendo: f(x)=(1/((Sqrt[x + 1])+(Sqrt[x]))).

Las siguientes instrucciones permiten comparar en una misma tabla los valores obtenidos al utilizar la dos expresiones de la función:

TableForm[
      Table[{i, (Sqrt[1.  + (10^i)]) - (Sqrt[(10^i)]), 
          (1/((Sqrt[1.  + (10^i)]) + (Sqrt[(10^i)])))}, {i, 3, 20}], 
      TableHeadings -> {None, {"Exponente base 10", "Por \
diferencia", "Por cociente"}}]

Ejemplo

Considere el cálculo de valores de la función f(x)=ln x-1 cuando x está próximo al número e. Intente diversos modos de resolver la cancelación que se produce.

Solución

h1[x_] := Log[x] - 1; 
h2[x_] := Log[x/E];    
h3[x_] := (((((Log[x]))^2) - 1)/(Log[x] + 1));    
h4[x_] := ((x - E)/E)-((((x - E))^2)/(2  (E^2))) +((((x - E))^3)/(3  (E^3)))  (* 
          aproximación por desarrollo de Taylor de orden 3 (ver abajo) *)

Series[Log[x] - 1, {x, E, 3}]

Normal[%]

lista = 
      Table[{n, N[E + (10^(-n)), 20], 
          InputForm[N[h1[E + (10^(-n))]]], 
          InputForm[N[h2[E + (10^(-n))]]], 
          InputForm[N[h3[E + (10^(-n))]]], 
          InputForm[N[h4[E + (10^(-n))]]]}, {n, 1, 
          20}]; 
TableForm[lista, 
      TableHeadings -> {None, {"Exp.", "Variable", "Por \
diferencia", "Por prop log", "Por conjugado", "Por \
Taylor"}}]

La siguiente sustitución permite entender por qué no se produce la cancelación con h4 y por qué no es útil hacer un desarrollo de Taylor de orden mayor a 3. Es más, en este caso se obtendrían los mismos resultados con los desarrollos de Taylor de órdenes 1 y 2.

h4[E + (10^(-n))]

Ejemplo

Resuelva los problemas de cancelación que se producen al evaluar la función f(x)=1/((Sqrt[x + 1])-(Sqrt[x])) en valores "grandes" de la variable x (como por ejemplo los valores de la forma x=(10^i), con 10 ≤ i ≤ 20).

Solución

lista = 
      Table[{i, 1/((Sqrt[1.  + (10^i)]) - (Sqrt[(10^i)]))}, {i, 1, 
          20}]; 
TableForm[lista, 
      TableHeadings -> {None, {"Exponente base 10", "Por \
diferencia"}}]

La cancelación puede ser resuelta si en la expresión f(x)=1/((Sqrt[x + 1])-(Sqrt[x])) racionalizamos, obteniendo que f(x)=(Sqrt[x + 1])+(Sqrt[x]).

En la siguiente tabla se pueden comparar los valores obtenidos al utilizar las dos expresiones de la función:

lista2 = 
      Table[{i, 
          1/((Sqrt[1.  + (10^i)]) - (Sqrt[(10^i)])), (Sqrt[1.  + (10^i)]) + (Sqrt[(10^i)])}, {i, 1, 20}]; 
\(TableForm[lista2, 
      TableHeadings -> {None, {"Exp base 10", "Por diferencia", "Por cociente"}}]\)

Resolver los ejercicios siguientes:

Ejercicio 1

Considere el cálculo de la función f(x) = (Sqrt[(x^3) + 4]) - 2 para x próximo a 0.  Muestre cómo para valores suficientemente pequeños de x se produce la cancelación (anulación) del valor de f(x). Compruebe que esta situación puede evitarse si se multiplica y divide por el conjugado para obtener que f(x) = ((x^3)/((Sqrt[(x^3) + 4]) + 2)).

Solución

Ejercicio 2

Encuentre un modo de evaluar las siguientes expresiones en los valores que se indican, de modo que puedan evitarse errores de cancelación.
a) 1  -  cos x en valores de x próximos a cero.
b) (Sqrt[(x^4) + 4]) - 2 en valores de x próximos a cero.
c) x((Sqrt[x + 1]) - (Sqrt[x])) en valores grandes de x.
d) x - (1/2) - (Sqrt[(x^2) - x + 3]) en valores grandes de x.
e) 2 - (Sqrt[(x^2) - 6  x + 13]) en valores próximos a x = 3.
Construya en cada caso una tabla en la que se comparen los resultados obtenidos al evaluar directamente las expresiones de los apartados anteriores en una sucesión finita de puntos que "tienda" al punto establecido (o al infinito, en su caso) con los obtenidos por el método utilizado para evitar la cancelación.

Solución

Ejercicio 3

Considere el cálculo de valores de la función f(x)=(e^x)-e cuando x se aproxima a 1. Resuelva las cancelaciones que se producen.

Solución