Los métodos numéricos son procedimientos matemáticos cuyo objetivo es la resolución numérica de problemas que carecen de expresión analítica para su resolución exacta. Estos procedimientos se expresan, en general, mediante algoritmos que especifican la secuencia de operaciones lógicas y aritméticas que conducen a la solución (normalmente aproximada) del problema planteado.
Los métodos numéricos se implementan usualmente mediante lenguajes de programación para su ejecución en un sistema computacional. Este sistema puede estar formado por un único ordenador, o por múltiples ordenadores conectados en red. Cualquiera que sea el caso, los sistemas computacionales tienen limitaciones inherentes que deben ser tenidas en cuenta en el diseño de los algoritmos.
Por ejemplo, los números irracionales como \(\sqrt{2}\), \(\pi\) o \(e\) tienen un desarrollo de infinitas cifras decimales, por lo que no pueden representarse exactamente en la memoria de un ordenador.
Normalmente cuando se aborda un problema científico, el primer paso es construir un modelo matemático que lo represente. Por ejemplo, si se quiere calcular el tiempo que tarda en llegar al suelo un objeto que se suelta desde una altura \(h\), utilizamos las ecuaciones del movimiento uniformemente acelerado:
\[h=\frac{1}{2}g t^2\] siendo \(g\) la aceleración de la gravedad.
Ahora bien, en la construcción de un modelo matemático se asumen siempre una serie de simplificaciones (distribución uniforme de temperaturas, material homogéneo, isótropo, regiones esféricas, ausencia de fricción con el aire …) que hacen que éste, ya desde su mismo planteamiento, se aparte (esperamos que ligeramente) de la situación real observada. Además, en el modelo pueden intervenir constantes (como \(g\) en el ejemplo anterior) cuyo valor no se determina de manera exacta. Depende del científico o del ingeniero afinar lo más posible en la elección del modelo matemático, de forma que la solución del mismo no se aparte significativamente de la solución del problema original.
Dado que los métodos numéricos se construyen a partir de los modelos matemáticos de los problemas que se pretende resolver, la primera fuente de error en la solución del problema depende de lo buena que sea la aproximación y calibración del modelo matemático empleado.
En esta asignatura no nos ocuparemos de esta fuente de error más que para citarla y para que el alumno-usuario-programador de los métodos numéricos sea consciente de su existencia.
Es el error debido a la propia forma en que está concebido el método numérico. En general, los métodos numéricos buscan dar soluciones aproximadas a los problemas, siendo el cálculo del error de aproximación una parte integrante del método.
Ejemplo :1 Aproximación del factorial de un número mediante la fórmula de Stirling. Se sabe que cuando \(n\) es grande, su valor puede aproximarse mediante:
\[n!\approx\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}\] estando la precisión de dicha aproximación acotada por: \[\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}e^{\frac{1}{12n+1}}<n!<\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}e^{\frac{1}{12n}}\]
Ejemplo 2: Aproximación de funciones mediante su desarrollo en serie de potencias. Por ejemplo, la función \(sen(x)\) puede aproximarse por su desarrollo de Taylor:
\[sen(x)=x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}-\frac{x^{7}}{7!}+\dots+\left(-1\right)^{n+1}\frac{x^{2n-1}}{\left(2n-1\right)!}\] La imposibilidad material de calcular los infinitos términos del desarrollo exige truncarlo con un número finito de términos, lo cual genera un tipo de error del método que, por su naturaleza, se denomina Error de Truncamiento. Por regla general, cuanto mayor sea el número de términos considerados, menor será el Error de Truncamiento.
El ordenador también introduce errores. La capacidad limitada de representación del ordenador impide representar de forma exacta constantes matemáticas (\(\pi\), \(e\)), ó números irracionales (\(\sqrt{2}\), \(\sqrt{3}\) ,\(\dots\)), que tienen un número infinito de decimales.
Aunque un número tiene infinitas cifras decimales, el ordenador sólo permite representarlo con un número finito de cifras significativas. En consecuencia, se produce un error que denominaremos Error de Representación o de redondeo
En el caso particular de octave/matlab, la representación interna de los números se hace con un total de 16 cifras decimales. Ello significa que dos números reales que sólo se diferencien a partir de la decimosexta cifra decimal serán considerados iguales a los efectos de su tratamiento con este programa:
>> a=0.123456789123456789;>> b=0.123456789123456788;>> a-bans = 0Así pues en el ordenador sólo puede ser representado un subconjunto \(\mathbb{F}\) del conjunto \(\mathbb{R}\) de los números reales. Los números de este subconjunto \(\mathbb{F}\) se denominan números de punto flotante o también números de coma flotante.
El error de representación o de redondeo también afecta a las operaciones aritméticas que se realicen con el ordenador; el resultado de cada operación se redondea de acuerdo con el número máximo de decimales que puede representar el ordenador; de esta forma en una sucesión de operaciones se pueden ir produciendo sucesivos errores, que se van arrastrando (propagando), lo que ocasiona una pérdida general de exactitud, puesto que se van perdiendo las cifras desechadas en cada paso.
Podemos ilustrar la importancia del error de propagación con un ejemplo muy conocido. Se trata del cálculo de la expresión (\(\left(1+10^{-20}-1\right)\cdot 10^{20}\)). Esta expresión, operada en el orden que se indica, da en la mayoría de los ordenadores como resultado 0, en contraste con su resultado correcto que es 1.
La explicación, como es fácil de adivinar, está en la operación \(1+10^{-20}\), donde, al tener los sumandos tanta diferencia en magnitud, para la capacidad de representación del ordenador su suma se redondea a 1, con las consiguientes consecuencias sobre el resto de las operaciones de la expresión, dando finalmente el resultado de 0.
Por tanto es importante tener en cuenta tanto el número de operaciones como el orden en que se realizan para evitar –o minimizar– el impacto del error de propagación.
Muchos métodos numéricos requieren la repetición indefinida de un proceso de cálculo. Es cuando la expresión matemática del método especifica que la solución del problema se obtiene como límite de una sucesión, esto es:
\[\textrm{Solucion =}\underset{n\rightarrow\infty}{\lim}\left\{x_i\right\}_{i=1}^n\]
La imposibilidad material de extender el proceso de cálculo indefinidamente, obliga a limitar el número de términos de la sucesión a un valor finito \(\left\{x_i\right\}_{i=1}^{n_0}\) considerando que para este valor, la solución aproximada alcanzada es un resultado suficientemente satisfactorio. Por tanto, se produce otro tipo de error que denominaremos, Error Residual, que representa la diferencia debida a los términos no calculados.
El Error numérico total será la suma de todos los errores cometidos. Limitándonos a analizar los errores, desde el modelo matemático hasta que se obtiene la solución numérica, las dos fuentes principales de error y las medidas para atenuarlas (que no eliminarlas) suelen ser:
Reducir el error del método, principalmente el error de truncamiento. Para reducir este error se procede a la mejora de la exactitud del método numérico empleado, lo que suele acarrear una formulación más compleja y, generalmente, un incremento en el número de operaciones aritméticas, lo que a su vez hace aumentar el error de redondeo, además de incrementar el tiempo de cálculo.
Reducir el error de redondeo. Como medidas para reducirlo, se presta atención al orden y a la forma en que se realizan las operaciones, y se incrementa el número de cifras significativas con las que opera el ordenador (por ejemplo, empleando números de doble precisión). En este último caso, esto supone aumentar las necesidades de memoria y almacenamiento del ordenador. En casos extremos, ello puede obligar a reformular completamente el método.
Como vemos, el comportamiento de los errores de truncamiento y redondeo responde en general a curvas con tendencias contrarias, por lo que en muchas ocasiones el remedio para un tipo de error incrementa el otro. Se hace necesario llegar a una situación de compromiso entre los dos tipos de error, lo que requiere cierto nivel de experiencia por parte del analista, y muchas veces obliga muchas a una serie de pruebas de ensayo‐error hasta llegar a la mejor solución.
En los problemas de cálculo numérico, si \(x\) es la solución (exacta) a un problema y \(x^{*}\) es la solución aproximada obtenida, el error absoluto cometido en la aproximación es simplemente \(\left|x-x^{*}\right|\). En la práctica suele ser mucho más interesante calcular el error relativo. Es evidente que un error de una milésima apenas tiene importancia si el valor de \(x\) se mide en la escala de los millones; pero el mismo error de una milésima es muy grave si el valor de \(x\) se mide en la escala de las millonésimas.
El error relativo se define como:
\[\varepsilon_{r}=\frac{\left|x-x^{*}\right|}{\left|x\right|}\]
En general los ordenadores almacenan los números de la forma siguiente:
\[x=\left(-1\right)^{s}\cdot\left(0.a_{1}a_{2}\dots a_{t}\right)\cdot\beta^{e}\] donde:
\(s\) vale 0 si el número \(x\) es positivo, ó 1 si \(x\) es negativo.
El número \(\left(0.a_{1}a_{2}\dots a_{t}\right)\) puede expresarse también como \((a_1a_2\dots a_t)\cdot \beta^{-t}\), donde \(0<a_1\le9\), y \(0\le a_i \le 9\) para \(i=2,\dots, t\). En esta representación, \((a_1a_2\dots a_t)\) un número entero llamado mantisa cuya longitud \(t\) es el número máximo de cifras de \(x\) que se almacenan (cifras significativas).
\(\beta\) es la base de numeración. En el sistema binario, la base de numeración es \(\beta=2\). En en sistema decimal es \(\beta=10\)
\(e\) es el exponente (valor entero). En matlab/octave, el valor de \(e\) puede estár entre -308 y +308.
Ejemplos:
\(0.00701=\left(-1\right)^{0}\cdot\left(0.701\right)\cdot10^{-2}\) tiene tres dígitos significativos.
\(1.00701=\left(-1\right)^{0}\cdot\left(0.100701\right)\cdot10^{1}\) tiene seis dígitos significativos.
\(-10.0701=\left(-1\right)^{1}\cdot\left(0.100701\right)\cdot10^{2}\) también tiene seis dígitos significativos.
En estos ejemplos vemos el significado del término valor de coma flotante: en realidad el mismo número se escribe de otra forma haciendo “flotar” la coma hasta colocarse delante de la primera cifra distinta de cero, ajustando a continuación convenientemente el valor del exponente.
En la práctica, aunque octave/matlab (u otros programas) guardan (y operan internamente con) números con mantisas de 16 cifras, normalmente los valores que se muestran por pantalla se muestran con menos dígitos. Por defecto, octave/matlab muestra 5 cifras significativas. Ello significa que si el número está guardado con más cifras, se redondea a cinco cifras: es decir, se muestran exactamente las cuatro primeras y la quinta se queda igual si la sexta cifra es menor que 5, y se le suma una unidad si la sexta cifra es mayor o igual que cinco.
Ejemplos:
>> x1=123.7639987;
x1x1 = 123.76>> x2=0.00038746534;
x2x2 = 3.8747e-04Si deseamos que octave/matlab muestren más cifras decimales podemos especificar alguno de los siguientes comandos:
format long: se muestran hasta 15 dígitos significativos.
format short: se muestran hasta 15 dígitos significativos
format long e: se muestran hasta 15 dígitos significativos, con notación exponencial.
format short e: se muestran hasta 5 dígitos significativos, con notación exponencial
El efecto de estos comandos es el siguiente:
>> x=0.12345123451234512345;
xx = 0.12345>> format long
xx = 0.123451234512345>> format short
xx = 0.12345>> format long e
xx = 1.23451234512345e-01>> format short e
xx = 1.2345e-01Si queremos especificar una precisión determinada, usamos el comando output_precision(n) donde n indica el número de cifras significativas que queremos que se muestren:
>> output_precision(12)
xx = 1.23451234512e-01Téngase en cuenta que la precisión de matlab/octave es de 16 cifras; si se especifica una precisión mayor con el comando anterior, el programa se “inventa” las cifras por encima de la décimosexta:
>> output_precision(20)
xx = 1.2345123451234511769e-01
Los valores más pequeño y más grande que se pueden representar en octave/matlab, pueden obtenerse mediante:
>> realmaxans = 1.7977e+308>> realminans = 2.2251e-308Si como resultado de alguna operación se obtiene algún valor inferior a realmin octave/matlab devuelve un 0; si se obtiene un valor mayor que realmax octave/matlab devuelve la cadena Inf.
Por lo expuesto en los apartados anteriores, las soluciones numéricas obtenidas con el ordenador son casi siempre soluciones aproximadas, incluso en muchos casos aquellas que corresponden a un procedimiento teóricamente exacto. Esta es la razón por la que a este tipo de método, teóricamente exacto, se le suele denominar método directo, evitando la denominación de “exacto”.
Por ejemplo, si se desea resolver la ecuación \(7x=5\), el método de resolverla es exacto; basta calcular \(x\) como \(x=\frac{5}{7}\). Sin embargo, la solución es aproximada:
>> x=5/7;
xx = 0.71429En general, lo que caracteriza a los métodos directos es que el número máximo de operaciones que se requiere para obtener la solución es conocido a priori. Tal vez, el ejemplo más conocido de método directo sea el Método de Eliminación de Gauss para resolver un sistema de ecuaciones lineales. Pocos más son los métodos de estas características.
Téngase en cuenta que, aunque obtengamos una solución siguiendo el procedimiento del método directo, la aritmética limitada del ordenador introduce distorsiones de forma que la solución final no es tan exacta como debería, o incluso aún peor, puede que no se parezca en nada a la solución de nuestro problema. Como la exactitud requerida en la solución varía de un problema a otro, resultan especialmente interesantes los métodos numéricos que permiten obtener los resultados con tanta exactitud como se desee.
Tal como se ha comentado, las soluciones obtenidas en los métodos numéricos son, en la gran mayoría de las ocasiones, soluciones aproximadas. Interesa por tanto medir el grado de aproximación alcanzado con nuestra solución, o lo que es lo mismo, el error de aproximación. Si \(a\) es la solución exacta del problema y \(x\) es la solución numérica, el error absoluto de aproximación es:
\[\varepsilon_{abs}=\left|a-x\right|\]
y el error relativo:
\[\varepsilon_{rel}=\frac{\left|a-x\right|}{\left|x\right|}\]
En general suele ser preferible controlar el error relativo, ya que el error absoluto será más o menos importante dependiendo de la magnitud de la solución \(a\).
Obviamente, en la práctica la solución exacta \(a\) no se conoce, por lo que no es posible conocer el valor del error. Por ello, en el caso de los métodos directos, para comprobar la validez de la solución obtenida, normalmente se recurre a sustituirla en las ecuaciones con las que se planteó el problema y comprobar que la solución obtenida las verifica.
La gran mayoría de los métodos numéricos que se estudian y aplican en la práctica, corresponden a procedimientos que proporcionan aproximaciones a la solución, pero cuya precisión se puede establecer arbitrariamente. A este tipo de métodos se les denomina métodos de aproximaciones sucesivas , ó también métodos iterativos. A cada intento de obtener una nueva aproximación se le denomina “iteración”, y a la acción, “iterar”. Con este tipo de métodos no es predecible el número total de operaciones a realizar porque, entre otras cosas, depende de la precisión con la que se quiere la solución.
Por esta razón es aconsejable poner algún límite al número total de operaciones que estamos dispuestos a permitir (por ejemplo, limitando el número máximo de iteraciones). En cuanto a la influencia de los errores de redondeo sobre las soluciones numéricas obtenidas con un método iterativo, recordemos que la precisión puede establecerse de forma arbitraria, por lo que, en caso de que la aproximación no sea satisfactoria, aumentaremos la exigencia de precisión, lo que supondrá continuar iterando hasta alcanzar una solución numérica más precisa (siempre que el crecimiento de los errores de redondeo se mantenga acotado).
Un método iterativo es un método que genera una sucesión de soluciones aproximadas:
\[x_0, x_1, x_2, \dots, x_n, \dots\]
tales que en el límite se alcance la solución exacta \(a\) del problema:
\[x_{i}\underset{_{i\rightarrow\infty}}{\longrightarrow}a\] En cada paso, el valor de \(x_i\) se obtiene a partir del valor de \(x_{i-1}\), esto es, el método requiere definir una función \(F\) tal que \(x_i = F(x_{i-1})\).
Obviamente, no es posible calcular infinitos términos de la sucesión \(x_i\). Por tanto debemos establecer una condición de parada del método. Dicha condición consiste en tomar como buena la solución \(x_n\) que cumpla que, dado un valor \(\varepsilon_a\) preespecificado:
\[\left|x_{n}-x_{n-1}\right|\le\varepsilon_{a}\] si nos interesa controlar el error absoluto, o bien que cumpla: \[\frac{\left|x_{n}-x_{n-1}\right|}{\left|x_n\right|}\le\varepsilon_{r}\] si lo que nos interesa es controlar el error relativo. El valor \(\varepsilon_r\) también se fija a priori. Los valores \(\varepsilon_a\) y \(\varepsilon_r\) reciben el nombre de tolerancia (absoluta y relativa, respectivamente).
El método iterativo se dice convergente si \(\underset{i\rightarrow\infty}{lim}x_{i}=a\). Llamando \(e_n=\left|x_{n}-x_{n-1}\right|\), para que el método sea convergente debe ocurrir que a partir de cierto valor \(n_0\) en adelante:
\[\left|e_{n}\right|<\left|e_{n-1}\right|,\;\forall n>n_0\]
Tal como hemos definido los métodos iterativos, la solución aproximada en cada paso, \(x_i\), depende de la aproximación obtenida en el paso anterior, \(x_{i-1}\). Ello significa que para poner en marcha el método es necesario especificar una aproximación inicial \(x_0\). La determinación de dicho valor requiere que el analista tenga un buen conocimiento del problema a resolver y que ello le permita conjeturar una aproximación inicial \(x_0\) razonable.
En muchos casos, el comportamiento del método iterativo depende fuertemente de la aproximación inicial \(x_0\), pudiendo llegar a darse en caso de que si dicha aproximación está mal elegida el método no converge a la solución del problema. Elegir la aproximación inicial no es siempre sencillo y depende del problema a resolver y de los conocimientos e intuición del analista-programador.
Un método iterativo no sólo tiene que funcionar, sino que además debe hacerlo rápidamente. Muchas veces es posible definir varios métodos alternativos para resolver un mismo problema y ello hace que sea necesario disponer de herramientas que permitan decidir qué método es más rápido.
Si \(a\) es la solución exacta del problema y \(x_i\) es la solución aproximada obtenida en la etapa \(i\), el error cometido en dicha etapa es \(e_i=\left|x_i-a\right|\). Para que el método sea convergente debe ocurrir que:
\[\underset{i\rightarrow\infty}{lim}\;{e_{i}}=0\] Como hemos señalado antes, ello requiere además que a partir de un \(n_0\) en adelante:
\[\left|e_{n}\right|<\left|e_{n-1}\right|,\;\forall n>n_0\]
El método se dice que tiene un orden de convergencia p o simplemente que el método es de orden p, siendo \(p\ge 1\) cuando: \[\left|e_n\right| \approx \lambda\left|e_{n-1}\right|^p\] siendo \(\lambda\) una constante independiente de \(n\). Obviamente cuánto mayor sea \(p\) más rápidamente disminuirá el error. Por tanto al elegir entre dos métodos iterativos alternativos para el mismo problema, será preferible el que tenga un orden de convergencia mayor.
Un método iterativo es estable si pequeñas variaciones en los datos iniciales producen pequeñas variaciones en los resultados. Si no sucede así el método es inestable.
La estabilidad numérica está relacionada con el problema en sí mismo (condicionamiento numérico del problema), y con la propagación y crecimiento del error de redondeo.
A modo de ejemplo, observemos los siguientes sistemas de ecuaciones lineales:
\[\left(\begin{array}{cc} 7 & 2\\ 1 & 0.28 \end{array}\right)\left(\begin{array}{c} x_{1}\\ x_{2} \end{array}\right)=\left(\begin{array}{c} 6\\ 2 \end{array}\right)\Rightarrow\textrm{Solucion: }\left(\begin{array}{c} x_{1}\\ x_{2} \end{array}\right)=\left(\begin{array}{c} 58\\ -200 \end{array}\right)\]
\[\left(\begin{array}{cc} 7 & 2\\ 1 & 0.28\boldsymbol{5} \end{array}\right)\left(\begin{array}{c} x_{1}\\ x_{2} \end{array}\right)=\left(\begin{array}{c} 6\\ 2 \end{array}\right)\Rightarrow\textrm{Solucion: }\left(\begin{array}{c} x_{1}\\ x_{2} \end{array}\right)=\left(\boldsymbol{\begin{array}{c} 458\\ -1600 \end{array}}\right)\]
Vemos como un pequeño cambio en la matriz de coeficientes del sistema induce grandes cambios en la solución. En este ejemplo las dos soluciones son exactas si sólo hay que resolver este problema. Ahora bien, si esta matriz de coeficientes apareciese en el curso de un método iterativo, y el vector de soluciones \(x\) ha de ir calculándose en cada iteración, pequeños errores de redondeo en alguno de los coeficientes se pueden ir amplificando y propagando, dando lugar a que la solución final proporcionada por el método sea incorrecta.
En casos como éste se dice que el problema está mal condicionado. La solución de los problemas mal condicionados es compleja y requieres muchas veces su replanteamiento para que no dependa de términos en los que los errores de redondeo puedan tener gran efecto. En general los errores de redondeo se consideran aceptables si crecen como mucho de manera lineal, esto es, existen un valor \(n_0\) y una constante \(C\) tales que:
\[\left|e_{n}\right|\approx Cn\left|e_{0}\right|,\,\,\,n\ge n_{0}\]