"25provocandoerror_1.gif"
por José Luis Gómez Muñoz       
http://www.globalcomputing.com.mx/

Porqué debemos ser muy cuidadosos con el cálculo numérico en cualquier paquete o lenguaje computacional

Esta es la representación decimal de 1/3:

"25provocandoerror_2.gif"

"25provocandoerror_3.gif"

Tres veces un tercio es uno

"25provocandoerror_4.gif"

"25provocandoerror_5.gif"

Sin embargo, tres veces la representación decimal de 1/3 no es uno.
Casi da uno, pero no es exactamente uno:

"25provocandoerror_6.gif"

"25provocandoerror_7.gif"

Las computadoras pueden hacer cálculos efectuando millones cálculos numéricos en fracciones de segundo. En esos cálculos numéricos son por definición aproximaciones, pues números como 1/3 son representados como 0.333333, con cierto número de dígitos. Cada uno de estos cálculos tienes pequeñísimos errores, como obtener 0.999999 en lugar de uno. La acumulación de millones de esos errores pequeñísimos puede provocar al final del cálculo un error enorme, al grado de que el resultado sea totalmente erróneo.  

Este comportamiento no es un problema exclusivo de Mathematica. Es un problema de cualquier software que pueda hacer cálculos numéricos.

Aquí mostramos como un uso descuidado de Mathematica puede provocar malos resultados, sin que Mathematica tenga la capacidad de advertirnos que el resultado es erróneo. Esto no es una característica exclusiva de Mathematica, es una característica de cualquier programa computacional que haga cálculos numéricos.  

Cálculo correcto del determinante de una matriz de Hilbert

El cálculo del determinante de un tipo especial de matriz, la matriz de Hilbert, es un cálculo típico para poner a prueba los programas que hacen cálculos numéricos.
Esta es la matriz de Hilbert de 15 por 15, mostrada en un formato (MatrixForm) similar a la de los libros de texto:

"25provocandoerror_8.gif"

"25provocandoerror_9.gif"

Este es el resultado exacto del determinante de la matriz de Hilbert de 15 por 15, el resultado es un uno dividido por un entero enorme:

"25provocandoerror_10.gif"

"25provocandoerror_11.gif"

Esta es la representación decimal del determinante exacto de la matriz de Hilbert de 15 por 15.
Este número, que es el correcto, lo utilizaremos para compararlo con otros cálculos más abajo:

"25provocandoerror_12.gif"

"25provocandoerror_13.gif"

Provocando que Mathematica se equivoque al calcular el determinante de una matriz de Hilbert

Esta el la representación decimal de la matriz de Hilbert de 15 por 15:

"25provocandoerror_14.gif"

"25provocandoerror_15.gif"

A continuación se calcula el determinante de la representación decimal de la matriz.
Obsérvese que la diferencia entre este cálculo y el cálculo hecho más arriba en este documento es el uso de N[] para cambiar las fracciones como 1/4 en su equivalente decimal como 0.25
Uno podría creer que la respuesta debería ser la misma, ya que 1/4 "es lo mismo que" 0.25
El siguiente resultado puede ser distinto en tu propia computadora:

"25provocandoerror_16.gif"

"25provocandoerror_17.gif"

Podemos ver que los dos cálculos del determinante, el que se hizo inmediatamente arriba y se guardó en la variable numerico1, y el que se hizo más arriba y se guardó en la variable correcto, son números muy distintos. Observando los exponentes, puedes notar que la magnitud de uno de ellos es más de mil veces más grande que la magnitud del otro. Además, uno de ellos es negativo y el otro positivo:

"25provocandoerror_18.gif"

"25provocandoerror_19.gif"

¿Porqué falló de esa manera tan estrepitosa el segundo cálculo?
Una forma de explicarlo es notar que 1/3 no es exactamente lo mismo que 0.333333, porque se requerirían una cantidad infinita de 3 después del punto decimal, pero la memoria de la computadora es finita, y no puede contener un número infinito de digitos.
De forma similar, 1/7 (que aparece en la primera matriz) no es exactament lo mismo que 0.142857, para que fueran lo mismo, necesitaríamos un número infinito de dígitos. Las pequeñísimas diferencias se vuelven enormes en valor relativo (un cálculo resultó mil veces más grande que el otro, ver los exponentes) debido a las multiplicaciones involucradas en el cálculo del determinante (podemos decir que el error se multiplicó).

Evitando el error numérico en Mathematica

A diferencia de otros programas, a Mathematica se le puede indicar que use una precisión en sus cálculos mayor que la precisión definida en la computadora donde se trabaja. Esto se hace con un segundo argumento en el comando N[].
Por ejemplo, sin segundo argumento en N[], el número Pi es:

"25provocandoerror_20.gif"

"25provocandoerror_21.gif"

Con mayor precisión

"25provocandoerror_22.gif"

"25provocandoerror_23.gif"

Con todavía mayor precisión

"25provocandoerror_24.gif"

"25provocandoerror_25.gif"

Podemos entonces calcular el determinante de la representación decimal de la matriz de Hilbert con más precisión:

"25provocandoerror_26.gif"

"25provocandoerror_27.gif"

Todavía con más precisión:

"25provocandoerror_28.gif"

"25provocandoerror_29.gif"

Podemos ver que con precisión de 80, se obtiene una respuesta que coincide con la respuesta correcta:

"25provocandoerror_30.gif"

"25provocandoerror_31.gif"

Entonces, ¿Porqué Mathematica no hace todos sus cálculos numéricos con precisión de 80? Primero porque la precisión de 80 no es suficiente para todos los cálculos, puede haber cálculos que requieran todavía más precisión (por ejemplo, el cálculo del determinante de una matriz de Hilbert de 30 por 30). Además, mientras más precisión se use, más lentos serán los cálculos y más memoria requieren. Para que los cálculos sean rápidos, debe manejarse la precisión natural de la computadora, que es lo que hace el comando N[] cuando se usa con un solo argumento, como en N[Pi]. Este conflicto entre lentitud y memoria por un lado y precisión por el otro, es una parte natural e inevitable del trabajo de quien desea hacer cálculos numéricos. Mathematica ofrece la posibilidad de hacer cálculos exactos (manejando las fracciones en lugar de su representación decimal). Otros programas, como las hoja de cálculo, no ofrecen esa posibilidad, y los resultados numéricos que se obtienen de una hoja de cálculo pueden simplemente ser erróneos, sin que el programa de hoja de cálculo pueda advertirnos de esa posibilidad. Por ejemplo, intente calcular el determinante de la matriz de Hilbert de 15 por 15 en una hoja de cálculo. El cálculo puede ser muy lento, y el resultado será totalmente erróneo, y el programa de hoja de cálculo no nos advertirá de ninguna manera que el número es erróneo. Lo mismo puede pasar si escribe su propio programa para calcular determinantes en algún lenguaje como C.

El comando Rationalize[]

Podemos transformar un número decimal en una razón (división de un entero entre otro entero) mediante el comando Rationalize[]:

"25provocandoerror_32.gif"

"25provocandoerror_33.gif"

Si ponemos demasiados decimales, entonces Rationalize con un sólo argumento no puede encontrar una razón adecuada y responde con el mismo número decimal, sin cambios:

"25provocandoerror_34.gif"

"25provocandoerror_35.gif"

Rationalize con dos argumentos busca una razón que se cuya diferencia con el primer argumento sea menor que el segundo argumento:

"25provocandoerror_36.gif"

"25provocandoerror_37.gif"

Si exigimos mayor exactitud la respuesta es distinta, por ejemplo la siguiente razón está más cercana de 0.333333 (con seis dígitos 3) que 1/3:

"25provocandoerror_38.gif"

"25provocandoerror_39.gif"

Si el segundo argumento es cero, forzamos a Mathematica a obtener una razón para cualquier número decimal:

"25provocandoerror_40.gif"

"25provocandoerror_41.gif"

En Mathematica Pi representa el número π, que da la razón entre el perímetro de cualquier círculo y su diámetro:

"25provocandoerror_42.gif"

"25provocandoerror_43.gif"

Esta es una representación aproximada de π, usando la precisión natural de la computadora (Machine Precision) donde el comando es evaluado:

"25provocandoerror_44.gif"

"25provocandoerror_45.gif"

Otra representación aproximada de π, usando mayor presición:

"25provocandoerror_46.gif"

"25provocandoerror_47.gif"

Mathematica se niega a racionalizar π, y hay buenas razones para eso, ya que los matemáticos han demostrado que π es un número irracional, es decir, que no hay razón (división entre dos enteros) que iguale exactamente a la razón (división) entre el perímetro y el diámetro de un círculo:

"25provocandoerror_48.gif"

"25provocandoerror_49.gif"

Sin embargo podemos obtener una aproximación racional a π forzando a Mathematica a racionalizar una aproximación decimal a π

"25provocandoerror_50.gif"

"25provocandoerror_51.gif"

Algunas culturas antiguas, como los griegos clásicos, que no habían inventado todavía la representación de los números decimales, sólo podian representar razones (divisiones de enteros). De esa manera, para el número π usaban aproximaciones con números racionales, incluso aunque sabían que sólo eran aproximaciones. Por ejemplo Arquímedes demostró que π es un número entre "25provocandoerror_52.gif" y "25provocandoerror_53.gif", es decir, "25provocandoerror_54.gif". Para su demostración, usó polígonos ligeramente mayores que el circulo (polígonos circumscritos) y polígonos ligeramente menores que el círculo (polígonos inscritos). Al ejecutar el comando Rationalize sobre la aproximación decimal a π y con una tolerancia de 0.01, Mathematica de hecho obtiene una de las cotas de Arquímedes:

"25provocandoerror_55.gif"

"25provocandoerror_56.gif"

Primer ejercicio

El siguiente mensaje real es una duda de un usuario de Mathematica. Reproduce en tu sesión de Mathematica los problemas que este usuario reporta:
> Dear All,
>
> I don't understand the following behaviour of Solve.
>
> Consider the following system :
>
> Solve[{x*y==(a+2*b)/(c+2*d),1/Sqrt[2]==Sqrt[(e*y)/(z*f*g*h)],
> 2*Pi*i==0.9/(z*f)},{x,y,z}]
>
> Everything fine, I obtain :
>
> {{z -> 0.1432394487827058/(f*i),
>   x -> (1.419827298426263*^-9*(9.83403688*^9*a +
>             1.966807376*^10*b)*e*i)/((c + 2.*d)*g*h),
>   y -> (0.0716197243913529*g*h)/(e*i)}}
>
> But if I ask for the answer of almost the same (only a 4 in the
> denominator of the second equation), Solve isn't abble anymore to
> manage without using inverse functions... why?
>
> Solve[{x*y == (a + 2*b)/(c + 2*d),
> 1/Sqrt[2] == Sqrt[(e*y)/(4*z*f*g*h)],
>   2*Pi*i == 0.9/(z*f)}, {x, y, z}]
>
> Worse:
>
> If I have numerical values for a, b, c, d, e, f, g, h, i:
> a = 65/10^6;
> b = 1/10^3;
> c = 1.9;
> d = 0.19;
> e = 1/(2.5/10^3);
> v = 18;
> w = 8;
> g = (2*v)/((c + 2*d)*w);
> i = 3000;
> h = 0.2;
>
> Then Solve isn't able anymore! Mathematica thinks there is no
> solution. But there is one. I have to use Reduce or give the numerical
> values as rules after the Solve to find them.
>
> In[247]:= Solve[{x*y == (a + 2*b)/(c + 2*d), 1/Sqrt[2] ==
> Sqrt[(e*y)/(z*f*g*h)],
>   2*Pi*i == 0.9/(z*f)}, {x, y, z}]
>
> Out[247]= {}
>
> Can somebody explain this?
>
> Thanks in advance,
>
> Regards
>
>
> F.Jaccard
>

Segundo ejercicio

El siguiente mensaje real es la respuesta de otro usuario de Mathematica a la duda del ejercicio anterior. Usa las recomendaciones de este usuario avanzado, y verifica si efectivamente resuelven los problemas planteados en el mensaje del ejercicio anterior:

Its basically the same old problem as has been discussed here many times
before. The methods used by Solve are not really suitable for dealing
with equations with MachinePrecision approximate numbers.
Rationalize all your approximate numbers and you will get solutions that
you expect. Note that when you use Reduce on your last system you obtain
the message:

Reduce::"ratnz" :  "Reduce was unable to solve the system with inexact \
coefficients. The answer was obtained by solving a corresponding exact
system \ and numericizing the result"

which gives you the same advice as above. Clearly Solve uses a different
approach from Reduce. I am sure that Solve also at some point replaces a
non-exact system by an exact one (for example it would have to do so when
it uses GrobnerBasis) but it does it at a different stage than Reduce,
and in this particular case this approach fails to produce an answer
(which suggests that your system is numerically unstable).

PREGUNTA Y RESPUESTA OBTENIDAS DEL SIGUIENTE FORO DE USUARIOS DE MATHEMATICA:
http://mathforum.org/kb/forum.jspa?forumID=79

Tercer ejercicio

El siguiente mensaje real es otra duda de un tercer usuario de Mathematica. Reproduce en tu sesión de Mathematica los problemas que este tercer usuario reporta:
Dear group members,

Consider

f[a_, b_, c_, k_, t_] :=With[{α = a k, β = b k}, (x - α Cos[t])^2/a^2 + (y - (β Sin[t] + c) - c)^2/b^2 - 1 == 0];

df[a_, b_, c_, k_, t_] := D[f[a, b, c, k, t], t];

and execute

{x, y} /. Simplify@PowerExpand@Simplify@Solve[{f[1, 2, 1/2, 4/5, t], df[1, 2, 1/2, 4/5, t]}, {x, y}] // Chop // N

and

numsol=({x, y} /. Simplify@PowerExpand@Simplify@NSolve[{f[1, 2, 1/2, 4/5, t], df[1, 2, 1/2, 4/5, t]}, {x, y}] // Chop)

Simplify@PowerExpand@Chop@Simplify[numsol /. Cos[2 t] -> (1 - 2 Sin[t]^2)]

respectively.

As you see, Solve and NSolve yield two different solutions, with the
solution from Solve being the correct one, as can be verified by
plugging in to the equations -- the solution from NSolve does not
satisfy the second equation, but only the first. Can anyone explain
this behaviour to me?

Best wishes,
Sigmund Vestergaard

Cuarto ejercicio

El siguiente mensaje real es la respuesta de un cuarto usuario de Mathematica a la duda del ejercicio anterior. Usa las recomendaciones de este usuario avanzado, y verifica si efectivamente resuelven los problemas planteados en el mensaje del ejercicio anterior:

The reason is that Solve solves the equations symbolically so it has no problems with numerical stability while NSolve with machine precision does. If you use NSolve with high extended precision you will get (essentially) the same answer.

first =
     N[Chop[{x, y} /.
      Simplify[PowerExpand[Simplify[Solve[{f[1, 2, 1/2, 4/5, t],
                     df[1, 2, 1/2, 4/5, t]}, {x, y}]]]]]];

second = Chop[{x, y} /. FullSimplify[PowerExpand[
              Simplify[NSolve[{f[1, 2, 1/2, 4/5, t], df[1, 2, 1/2, 4/5, t]},
                  {x, y}, WorkingPrecision -> 100]]]]];


Simplify[Rationalize[first - second, 0]]

{{0, 2*(Sqrt[Sin[t]^2] - Sin[t])}, {0, 2*Sin[t] - 2*Sqrt[Sin[t]^2]}}


The fact that we do not get 0 is due to the use of PowerExpand. So

PowerExpand[%]
( {{0, 0},{0, 0}} )


Andrzej Kozlowski


PREGUNTA Y RESPUESTA OBTENIDAS DEL SIGUIENTE FORO DE USUARIOS DE MATHEMATICA:
http://mathforum.org/kb/forum.jspa?forumID=79

"25provocandoerror_57.gif"
Autor: José Luis Gómez Muñoz

     Global Computing S. A. de C. V.
Florencia 57 Piso 10-01
Col. Juárez C.P. 06600
México D.F.
México
+52-(0)55-5525-2215
Fax: +52-(0)55-5514-4225

Adriana Vadillo avadillo@mx.inter.net

Hector Vadillo  hector.vadillo@prodigy.net.mx

http://www.globalcomputing.com.mx/

Spikey Created with Wolfram Mathematica 6