Gu´ıa de Pr´ acticas:
An´ alisis Num´ erico de
Ecuaciones Diferenciales usando MATLAB
2
´Indice general Notaciones
5
Introducci´ on on
6
I
9
Elem Elemen ento toss b´ asicos de MATLAB
1. Introducci´ Introducci´ on a MATLAB. 1.1. Vectores en MATLAB . . . . . . . . . . . . . . . . . . 1.2. Matrices en MATLAB . . . . . . . . . . . . . . . . . . 1.3. Funciones de vectores . . . . . . . . . . . . . . . . . . . 1.4. Bucles . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1. Relaciones . . . . . . . . . . . . . . . . . . . . . 1.5. La inst instruc rucci´ ci´ on if . . . . . . . . . . . . . . . . . . . . . . 1.6. Ficheros ejecutables . . . . . . . . . . . . . . . . . . . . 1.7. Subrutinas . . . . . . . . . . . . . . . . . . . . . . . . . 1.8. 1.8. Caden adenas as de texto, xto, mens mensaj ajees de error rror,, entrada radass . . . . . 1.9. 1.9. Co Comp mpar aran ando do la efici eficien enci ciaa de algo algori ritm tmos os:: cput cputim imee . . . . 1.10. Formatos de salida . . . . . . . . . . . . . . . . . . . . 1.11. Gr´aficos . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11.1. Representaciones en 2D . . . . . . . . . . . . . 1.11.2. Representaciones en 3D . . . . . . . . . . . . . 1.12. 1.12. Resu Resume men n de func funcio ione ness ele eleme ment ntal ales es y matr matric ices es espec especia iale less 3
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
11 12 14 20 25 28 29 29 31 33 34 34 35 35 37 39
4
II
Introducci´ on a la construcci´ on de algoritmos
2. Sistemas de numeraci´ on; errores y sus fuentes 2.1. Sistemas de n´ umeros y conversiones . . . . . . . . . . . . . . . . 2.1.1. Sistemas de numeraci´ on en base q . . . . . . . . . . . . . 2.1.2. Conversi´ on de base decimal a base q . . . . . . . . . . . 2.1.3. Est´ andar IEEE de representaci´o n de n´ umeros enteros y flotante . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Errores y sus causas . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1. Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2. Fuentes de errores . . . . . . . . . . . . . . . . . . . . . . 2.3. Propagaci´ on de errores: condici´on y estabilidad. . . . . . . . . . 2.4. Eficiencia de un algor´ıtmo num´erico . . . . . . . . . . . . . . . . 2.4.1. Ejemplo: Evaluaci´ o n de polinomios, m´ e todo de Horner . ´ 2.5. EJERCICIOS PRACTICOS . . . . . . . . . . . . . . . . . . . .
41 43 . . . . . . 43 . . . . . . 43 . . . . . . 44 en coma . . . . . . 47 . . . . . . 50 . . . . . . 50 . . . . . . 50 . . . . . . 54 . . . . . . 56 . . . . . . 57 . . . . . . 60
5
NOTACIONES Notaciones generales S´ımbolo x = (x1 , x2 ,...,xN ) r = x = (x21 + x22 + ... + x2N )1/2
| |
∇u =
∂u ∂u ∂u ∂x 1 , ∂x 2 ,..., ∂x N
Espacios de funciones
C C C C C C
S´ımbolo (Ω) 0 (Ω) k (Ω) k 0 (Ω) (Ω) (Ω) 0 (Ω) = ∞ ∞
D
Significado Elemento de IRN M´odulo de x Gradiente de u
Significado Funciones continuas en Ω Funciones continuas en Ω con soporte compacto Funciones de clase k en Ω Funciones de k (Ω) con soporte compacto Funciones indefinidamente diferenciables en Ω Funciones de (Ω) con soporte compacto
C C
∞
6
Introducci´ on
7
8
Parte I Elementos b´ asicos de MATLAB
9
Cap´ıtulo 1 Introducci´ on a MATLAB. Vamos a optar por MATLAB para programar los algoritmos que aparecen asociados a los distintos m´etodos num´ ericos que estudiaremos a lo largo del curso. La raz´ on de ello es que, para prop´ositos docentes, existen ventajas en la utilizaci´o n de este software de c´alculo matem´atico respecto a lenguajes de programaci´ on como C o Fortran: por un lado, la sintaxis de programaci´on de MATLAB es bastante similar a la de cualquier lenguaje de alto nivel , lo que nos permitir´a adaptarnos a las peculiaridades de su programaci´ on sin demasiados problemas; por otro lado, podremos utilizar como test de nuestros propios algoritmos las rutinas espec´ıficas que incorpora MATLAB; por u´ltimo, la integraci´on de una serie de paquetes gr´aficos en el entorno de programaci´on de MATLAB facilita extraordinariamente la tarea de presentaci´on de resultados. Los ejemplos que se incluyen en las secciones siguientes se han obtenido con la versi´on 6.0 (R12) de MATLAB. Al tratarse de ejercicios sencillos, no es esperable diferencias significativas cuando se reproduzcan en versiones posteriores de MATLAB. La aritm´etica de trabajo en MATLAB es doble precisi´on (concepto que aclararemos posteriormente) aunque la versi´on 7.0 ha introducido la posibilidad de trabajar con aritm´etica entera y real de precisi´on simple. Una vez hemos accedido a MATLAB y estando situados en la ventana de comandos, comenzamos nuestro recorrido introductorio. En el texto que sigue a continuaci´on, cualquier l´ınea que comienza con dos signos >> se utiliza para denotar una l´ınea de comando MATLAB. 11
12
1.1.
Vectores en MATLAB
Casi todos los comandos b´ asicos en MATLAB implican el uso de vectores. Podemos definir un vector tecleando cada uno de sus elementos o, alternativamente y cuando es posible, podemos hacerlo especificando el primer elemento, un incremento y el ´ultimo elemento. Por ejemplo, para crear un vector cuyos elementos son 0, 2, 4, 6 y 8, podemos teclear: >> 0:2:8 ans = 0
2
4
6
8
MATLAB tambi´en guarda el u ´ ltimo resultado. En el ejemplo previo, se ha creado una variable “ans”. Para obtener el vector traspuesto, tecleamos: >> ans’ ans = 0 2 4 6 8
Para guardar los vectores creados de modo que seamos capaces de trabajar con ellos posteriormente, debemos darles nombre. Por ejemplo, para crear el vector fila v, tecleamos: >> v = [0:2:8] v = 0
2
4
6
8
13 >> v v = 0
2
4
6
8
>> v; >> v’ ans = 0 2 4 6 8
Podemos darnos cuenta a partir del ejemplo anterior que si finalizamos una l´ınea con un punto y coma no se muestra el resultado. MATLAB permite tambi´en trabajar con elementos espec´ıficos del vector: >> v(1:3) ans = 0
2
>> v(1:2:4) ans = 0
4
>> v(1:2:4)’ ans =
4
14
0 4
Una vez especificada la notaci´on podemos realizar diversas operaciones. Por ejemplo: >> v(1:3)-v(2:4) ans = -2
1.2.
-2
-2
Matrices en MATLAB
Damos a continuaci´on una introducci´o n b´asica a la definici´on y manipulaci´ o n de matrices. La definici´on de una matriz es an´aloga a la definici´on de un vector. Podemos considerarla como una columna de vectores fila:
>> A = [ 1 2 3; 3 4 5; 6 7 8] A = 1 3 6
2 4 7
3 5 8
o como una fila de vectores columna:
>> B = [ [1 2 3]’ [2 4 7]’ [3 5 8]’] B = 1
2
3
15 2 3
4 7
5 8
A estas alturas de la introducci´on, si hemos estado reproduciendo los ejemplos anteriores tendremos, muy probablemente, una considerable cantidad de variables definidas en nuestro espacio de trabajo. Si queremos conocer esta informaci´on utilizaremos el comando whos: >> whos Name
Size
A B ans v
3 3 1 1
by by by by
3 3 3 5
Elements
Bytes
Density
9 9 3 5
72 72 24 40
Full Full Full Full
Complex No No No No
Volviendo a las matrices y a la hora de realizar operaciones b´asicas con las mismas, la notaci´ on que utiliza MATLAB es la usual en ´algebra lineal. Obviamente, las operaciones matriciales deben ser legales si no queremos tener problemas: >> v = [0:2:8] v = 0
2
4
6
8
>> A*v(1:3) ??? Error using ==> * Inner matrix dimensions must agree. >> A*v(1:3)’ ans = 16
16 28 46
Podemos trabajar con diferentes partes de una matriz, al igual que vimos que se pod´ıa hacer con vectores. De nuevo, debemos tener cuidado en hacer operaciones permitidas: >> A(1:2,3:4) ??? Index exceeds matrix dimensions. >> A(1:2,2:3) ans = 2 4
3 5
>> A(1:2,2:3)’ ans = 2 3
4 5
Aparte de las operaciones elementales (suma, producto), podemos hacer otras de inter´es como por ejemplo, obtener la inversa de una matriz (cuando esta operaci´on tenga sentido, si no obtendremos un error). Sin embargo, debemos tener cuidado puesto que las operaciones que se realizan pueden presentar errores de redondeo. En el ejemplo, la matriz A no es una matriz invertible, pero MATLAB devuelve una matriz como resultado. >> inv(A) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.565062e-18
ans =
17
1.0e+15 * -2.7022 5.4043 -2.7022
4.5036 -9.0072 4.5036
-1.8014 3.6029 -1.8014
Conviene hacer notar, en este punto, que MATLAB distingue entre may´usculas y min´ usculas. Esta es otra potencial fuente de problemas cuando trabajamos con algoritmos complicados: >> inv(a) ??? Undefined function or variable a.
Otra posible operaci´on es, por ejemplo, la obtenci´on de los valores propios aproximados de una matriz. Hay dos versiones de esta rutina: una encuentra los valores propios y la otra encuentra los valores y vectores propios. Si no recordamos cu´al es cu´al, podemos obtener m´as informaci´on tecleando eig en la l´ınea de comandos. >> eig(A) ans = 14.0664 -1.0664 0.0000 >> [v,e] = eig(A) v = -0.2656 -0.4912 -0.8295
0.7444 0.1907 -0.6399
-0.4082 0.8165 -0.4082
18
e = 14.0664 0 0
0 -1.0664 0
0 0 0.0000
>> diag(e) ans = 14.0664 -1.0664 0.0000
Existen tambi´en rutinas que permiten encontrar soluciones de ecuaciones. Por ejemplo, si Ax = b y queremos encontrar x, un modo “lento” de hacerlo es, simplemente, invertir A y realizar una multiplicaci´on por la izquierda sobre ambos lados de la ecuaci´on. Obviamente, hay m´etodos m´as eficientes y m´as estables para hacer esto (descomposiciones L/U con pivotes, por ejemplo). MATLAB tiene comandos especiales para hacer esto. MATLAB posee adem´as dos tipos diferentes de operadores / y . La acci´on del primer operador es la siguiente: x = A v A 1 v; la acci´on del segundo operador es: x = v/B vB 1 . Se dan ejemplos de su uso a continuaci´on:
\ ≡
−
−
\
>> v = [1 3 5]’ v = 1 3 5 >> x = A\v Warning: Matrix is close to singular or badly scaled.
≡
19 Results may be inaccurate. RCOND = 4.565062e-18
x = 1.0e+15 * 1.8014 -3.6029 1.8014 >> x = B\v x = 2 1 -1 >> B*x ans = 1 3 5 >> x1 = v’/B x1 = 4.0000
>> x1*B
-3.0000
1.0000
20 ans = 1.0000
3.0000
5.0000
Finalmente, si queremos borrar todos los datos del sistema y comenzar de nuevo utilizaremos el comando clear. Conviene utilizar este comando con cuidado porque MATLAB no pide una segunda opini´on ... >> clear >> whos
1.3.
Funciones de vectores
Es indudable que la gran ventaja de trabajar con MATLAB es la facilidad de manipulaci´ on de vectores y matrices. En este apartado comenzaremos con manipulaciones simples (suma, resta, multiplicaci´on). A continuaci´on mostramos c´omo se pueden definir operaciones relativamente complejas con un peque˜no esfuerzo. Comenzamos con la suma y resta de vectores. Definiremos dos vectores y a continuaci´on los sumaremos y restaremos: >> v = [1 2 3]’ v = 1 2 3 >> b = [2 4 6]’ b = 2
21 4 6 >> v+b ans = 3 6 9 >> v-b ans = -1 -2 -3
La multiplicaci´on de vectores y matrices sigue, l´ogicamente y como comentamos anteriormente, reglas estrictas, as´ı como la suma. En el ejemplo anterior los vectores son ambos vectores columna con tres componentes: >> v*b Error using ==> * Inner matrix dimensions must agree. >> v*b’ ans = 2 4 6 >> v’*b ans =
4 8 12
6 12 18
22
28
Hay ocasiones en las que queremos realizar una operaci´on sobre cada elemento de un vector o matriz. MATLAB permite hacer este tipo de operaciones. Por ejemplo, supongamos que queremos multiplicar cada componente de un vector v por la correspondiente componente del vector b. En otras palabras, Supongamos que queremos hallar v(1)*b(1), v(2)*b(2) y v(3)*b(3). Estar´ıa bien utilizar el s´ımbolo * puesto que estamos haciendo un tipo de multiplicaci´on. Sin embargo, como este s´ımbolo ha sido definido con otra funci´ on, debemos recurrir a otra cosa. Los programadores ocupados del desarrollo de MATLAB decidieron entonces utilizar los s´ımbolos .* para hacer esta operaci´on. De hecho, se puede emplear este s´ımbolo antes de cualquier s´ımbolo matem´atico para especificar a MATLAB que la operaci´on en cuesti´on debe tener lugar sobre cada entrada del vector. >> v.*b ans = 2 8 18 >> v./b ans = 0.5000 0.5000 0.5000
Puesto que hemos comenzado a hablar de operaciones no lineales continuemos con ellas: si pasamos un vector a una operaci´on matem´atica predefinida, obtendremos un vector del mismo tama˜no con entradas obtenidas realizando la operaci´on especificada sobre la correspondiente entrada del vector original: >> sin(v)
23 ans = 0.8415 0.9093 0.1411 >> log(v) ans = 0 0.6931 1.0986
La posibilidad de trabajar con estas funciones vectoriales es una de las ventajas de MATLAB. De este modo, podemos definir operaciones complejas r´apida y f´acilmente. En el siguiente ejemplo trabajamos con un vector con muchas componentes: >> x = [0:0.1:100] x = Columns 1 through 7 0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
99.7000
99.8000
99.9000
100.0000
(stuff deleted) Columns 995 through 1001 99.4000
99.5000
99.6000
>> y = sin(x).*x./(1+cos(x));
Adem´as de esta simple manipulaci´on de vectores, MATLAB nos permitir´a tambi´en representar gr´aficamente los resultados que obtengamos. Tecleando,
24 >> plot(x,y)
tendremos una representaci´on gr´afica de la funci´on antes considerada. Si tecleamos >> plot(x,y,’rx’)
obtenemos la misma gr´afica pero las l´ıneas son reempladas por puntos ro jos en forma de x. Para ver m´as opciones del commando plot, podemos teclear >> help plot
El comando help es, sin duda, el camino m´as corto para estar seguro de la sintaxis de un determinado comando de MATLAB. La notaci´on compacta permitir´a realizar un gran n´umero de operaciones utilizando pocos comandos. Veamos el siguiente ejemplo: >> >> >> >>
coef = zeros(1,1001); coef(1) = y(1); y = (y(2:1001)-y(1:1000))./(x(2:1001)-x(1:1000)); whos Name Size Elements Bytes ans b coef v x y
3 3 1 3 1 1
by by by by by by
1 1 1001 1 1001 1000
3 3 1001 3 1001 1000
Grand total is 3011 elements using 24088 bytes >> coef(2) = y(1); >> y(1) ans = 0.0500
24 24 8008 24 8008 8000
Density Full Full Full Full Full Full
Complex No No No No No No
25 >> y = (y(2:1000)-y(1:999))./(x(3:1001)-x(1:999)); >> coef(3) = y(1); >> >>
A partir de este algoritmo podemos encontrar el polinomio de Lagrange que interpola los puntos definidos anteriormente (vector x). Por supuesto, con tantos puntos el proceso puede resultar algo tedioso. Afortunadamente MATLAB dispone de un modo sencillo de hacer tareas mon´otonas, como veremos a continuaci´on.
1.4.
Bucles
En esta secci´on veremos c´omo utilizar los bucles for y while. En primer lugar, discutiremos el uso del bucle for con ejemplos para operaciones fila sobre matrices. A continuaci´ on, mostraremos el uso del bucle while. El bucle for permite repetir ciertos comandos. Todas las estructuras de bucles en MATLAB comienzan con una palabra clave (como for o while) y terminan con un end (parece sencillo, ¿no?). Veamos un ejemplo trivial: >> for j=1:4, j end
j = 1
j = 2
j =
26 3
j = 4 >>
Otro ejemplo: >> v = [1:3:10] v = 1
4
7
10
3
4
>> for j=1:4, v(j) = j; end >> v v = 1
2
´ Este es un ejemplo simple y una demostraci´on bonita de c´omo funcionan los bucles for. Sin embargo, no se debe utilizar en la pr´actica: la notaci´on utilizada en la primera l´ınea es mucho m´a s r´apida que el bucle. Un mejor ejemplo se presenta a continuaci´on, donde realizamos operaciones sobre las filas de una matriz. Si queremos comenzar en la segunda fila de una matriz y restar la fila previa y repetir esta operaci´on sobre las siguientes filas, un bucle for puede ocuparse de esto: >> A = [ [1 2 3]’ [3 2 1]’ [2 1 3]’] A =
27
1 2 3
3 2 1
2 1 3
>> B = A; >> for j=2:3, A(j,:) = A(j,:) - A(j-1,:) end A = 1 1 3
3 -1 1
2 -1 3
1 1 2
3 -1 2
2 -1 4
A =
Presentamos a continuaci´on un ejemplo m´as realista (implementaci´on de eliminaci´on Gaussiana): >> for j=2:3, for i=j:3, B(i,:) = B(i,:) - B(j-1,:)*B(i,j-1)/B(j-1,j-1) end end B = 1 0 3
3 -4 1
2 -3 3
28
B = 1 0 0
3 -4 -8
2 -3 -3
1 0 0
3 -4 0
2 -3 3
B =
La forma general de un bucle while es >> while (condiciones) (operaciones) end
Las operaciones se repetir´an mientras que las condiciones sean ciertas. Por ejemplo, dado un n´ umero a, el siguiente bloque permite calcular y mostrar el entero no negativo m´as peque˜ no tal que 2n a:
≥
>> n=0; >> while 2^n < a n=n+1; end >> n
1.4.1.
Relaciones
Los operadores de relaci´on en MATLAB son:
29
< > <= >= == =
∼
menor que mayor que menor o igual que mayor o igual que igual que distinto a
D´emonos cuenta que el s´ımbolo “==” se utiliza en lugar de “=” en una relaci´ on. Podemos conectar varias relaciones utilizando los operadores l´ogicos: & y o no
| ∼ 1.5.
La instrucci´ o n if
La forma general de una instrucci´on if simple es: >> if (condiciones) (operaciones) end
Las operaciones se realizar´an u ´nicamente si se cumplen las condiciones especificadas. Se admiten las ramificaciones m´ultiples como puede verse en el siguiente ejemplo: >> if n <0 a=1; elseif n<5 a=2; else a=3; end
1.6.
Ficheros ejecutables
En esta secci´on introducimos los conceptos b´asicos para crear ficheros ejecutables. Los ficheros ejecutables son ficheros de texto que incluyen una serie de comandos MATLAB.
30 Si una tarea MATLAB la vamos a ejecutar muchas veces, es buena idea escribir un fichero con estos comandos para poder ejecutarlos tantas veces como queramos. La edici´on del fichero ejecutable la realizamos con un editor cualquiera. Si nos resulta m´a s c´omodo, podemos utilizar el editor que incorpora MATLAB y al que invocaremos desde la l´ınea de comandos como: >>edit
Los ficheros ejecutables de MATLAB (llamados ficheros M) deben tener como extensi´on “.m”. En el ejemplo que damos a continuaci´on, crearemos un programa que calcula el factorial de 6: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Este es un programa no muy util, %que calcula el factorial de 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n=6; fac=1; for i=2:n fac=fac*i; end fac
Si guardamos esto en el fichero fac.m en el directorio de trabajo (o cualquier otro incluido en el “path”) y tecleamos el comando fac , obtenemos >> fac fac = 720
Las lineas tras el s´ımbolo “ %” son l´ıneas de comentario, que se conviene utilizar como explicaci´on del programa. El comando help sirve para mostrar esas l´ıneas: >> help fac %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Este es un programa no muy util, que calcula el factorial de 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31 En efecto, este no es un programa muy ´util, en primer lugar porque el propio MATLAB tiene su comando para calcular el factorial de n´umeros enteros:
factorial(6)
ans = 720
y en segundo lugar porque nuestro programa s´olo calcula el factorial de 6. Para poder calcular el factorial para distintos n´umeros deberemos crea una subrutina o funci´on MATLAB.
1.7.
Subrutinas
Si ahora escribimos en un fichero con nombre facf.m los siguientes comandos %Esta es una funcion para calcular %el factorial de n. %El valor de entrada es n function fac=facf(n) fac=1; for i=2:n fac=fac*i; end
habremos definido una funci´on que podemos utilizar tal como lo hacemos con los comandos intr´ınsecos de MATLAB. Por ejemplo, tecleando en la linea de comandos facf(6) tenemos: >> facf(6) ans = 720
Las funciones pueden tener varios par´ametros de entrada y/o salida. Por ejemplo, la siguiente es una funci´on que, dados dos vectores con la misma longitud, devuelve dos 2 valores (es decir, la subrutina implementa una funci´on f : Rn Rn R ).
× →
32 %Esta es una funcion que, dados %dos vectores de la misma longitud %calcula dos valores reales a y b %a=sin(x*y’), b=a*x*x´ function [a,b]=funci(x,y) a=sin(x*y’); b=a*x*x’;
Guardamos este fichero como funci.m y, como prueba, ejecutamos en la l´ınea de comandos: >> x=1:1:5 x = 1
2
3
4
5
>> y=0:0.1:0.4 y = 0
0.1000
0.2000
0.3000
0.4000
>> [f1,f2]=funci(x,y); >> f1 f1 = -0.7568 >> f2 f2 = -41.6241
Por supuesto, si las matrices de entrada x e y no son vectores de la misma longitud, el programa puede fallar. En el siguiente ejemplo, la primera llamada a funci es correcta, pero no as´ı la segunda:
33 >> [a,b]=funci(1:1:5,0:1:4); >> [a,b]=funci(1:1:5,0:1:5); ??? Error using ==> * Inner matrix dimensions must agree. Error in ==> c:\matlab\temporal\funci.m On line 6 a=sin(x*y’);
==>
Importante: – Fij´emonos que la sintaxis de las rutinas es: function [output1,output2,...]=nombre(input1,input2,...)
– Si queremos crear un fichero independiente que contenga a la rutina, el nombre que asignemos al fichero debe ser el mismo que el de la rutina con extensi´on “.m”. Por supuesto, las funciones as´ı definidas pueden ser utilizadas tantas veces como sea necesario dentro de otras funciones y programas. Una norma muy conveniente que conviene seguir es ser´a escribir todos los programas en ficheros de texto (utilizando un editor simple) y procurar que estos sean lo m´as estructurados posibles, utilizando subrutinas (funciones) para implementar tareas independientes. Se trata de programar de la forma m´as modular y estructurada posible.
1.8.
Cadenas de texto, mensajes de error, entradas
Se pueden introducir cadenas de texto en MATLAB si van entre comillas simples. Por ejemplo, la instrucci´on >> s=’Mi asignatura preferida es Calculo Numerico’
asigna a la variable s la anterior cadena de texto. Se pueden mostrar cadenas de texto utilizando la funci´on disp. Por ejemplo: >> disp(’En particular, me encantan las practicas de la asignatura’)
Los mensajes de error pueden (deben) mostrarse utilizando la funci´on error
34 >> error(’Fue un error no acabar las memorias a tiempo’)
puesto que cuando se utiliza en un fichero “.m”, interrumpe la ejecuci´on del mismo. Podemos, asimismo, en un fichero “.m” pedir que un usuario introduzca datos de forma interactiva utilizando la funci´on input. Cuando, por ejemplo, durante la ejecuci´ on de un programa aparece la l´ınea >>ifaltas =input(’Cuantas veces has faltado a practicas?’)
el citado mensaje de texto aparece en pantalla y la ejecuci´on del programa se interrumpe hasta que el usuario teclea el dato de entrada. Despu´es de presionar la tecla de “return”,el dato es asignado a la variable “ifaltas” y se reanuda la ejecuci´on del programa.
1.9.
Comparando la eficiencia de algoritmos: cputime
Una forma adecuada de determinar la eficiencia de un algoritmo es medir el tiempo de CPU transcurrido en la ejecuci´on del mismo. Este tiempo puede obtenerse utilizando el comando cputime. De hecho, la secuencia >> cpu1=cputime, (cualquier conjunto de operaciones), cpu2=cputime-cpu1
nos proporcionar´ a el tiempo de cpu invertido en la ejecuci´on del mencionado conjunto de operaciones.
1.10.
Formatos de salida
Mientras que todos los c´alculos en MATLAB se realizan en doble precisi´on, el formato de los datos de salida puede ser controlado con los siguientes comandos: Punto fijo y 4 decimales (es el que hay por defecto) f ormat short Punto fijo y 14 decimales format long format short e notaci´on cient´ıfica con 4 decimales format long e notaci´on cient´ıfica con 14 decimales
35
1.11.
Gr´ aficos
Aunque ya hemos mencionado anteriormente la utilizaci´on del comando plot, vamos a dar en esta secci´on alg´ un detalle adicional sobre las posibilidades gr´ aficas de MATLAB. MATLAB permite generar representaciones gr´aficas de curvas en 2D y 3D. Los comandos b´ asicos con los que nos manejaremos ser´an plot, plot3, mesh y surf .
1.11.1.
Representaciones en 2D
El comando plot crea gr´aficos de curvas en el plano x-y; si x e y son vectores de la misma longitud, el comando plot(x,y) abre una ventana gr´afica y dibuja en el plano x-y los elementos de x versus los elementos de y. Podemos, por ejemplo, dibujar el gr´afico de la funci´on seno en el intervalo -4 a 4 con los siguientes comandos: x=-4:.01:4; y=sin(x); plot(x,y)
El vector x es una partici´on del dominio en intervalos de longitud 0.01; el vector y es un vector que proporciona los valores del seno en los nodos de la partici´on. Como segundo ejemplo vamos a dibujar el gr´afico de y = e x en el intervalo -1.5 a 1.5: −
2
x=-1.5:.01:1.5; y=exp(-x.^2); plot(x,y)
D´emonos cuenta de que la operaci´on de elevar al cuadrado est´a precedida por un punto para que opere sobre cada una de las componentes del vector x. MATLAB posee adem´as la utilidad fplot para representar de forma eficiente y simple el gr´afico de una funci´on: para representar el ejemplo anterior podemos, de forma alternativa, definir la funci´on en un fichero M (que llamaremos, por ejemplo, expcu.m): function y=expcu(x) y=exp(-x.^2)
Entonces el comando fplot(’expcu’,[-1.5,1.5])
nos proporcionar´a el gr´afico en cuesti´on. Podemos generar tambi´en gr´aficos de curvas definidas param´etricamente. Por ejemplo,
36 t=0:.001:2*pi; x=cos(3*t); y=sin(2*t); plot(x,y)
De forma complementaria, podemos asignar a los gr´aficos: t´ıtulos, etiquetas en los ejes y texto (en la zona del gr´afico). Para ello utilizaremos los comandos
title xlabel ylabel gtext text
t´ıtulo del gr´afico etiqueta del eje x etiqueta del eje y sit´ ua texto sobre el gr´afico utilizando el rat´on sit´ ua texto en las coordenadas especificadas
Ejemplo, el comando: title(’Hola caracola’)
proporciona al gr´afico el t´ıtulo en cuesti´on. Los ejes del gr´afico son, por defecto, autoescalados. Para modificar esto podemos utilizar el comando axis: axis([xmin,xmax,ymin,ymax])
El comando axis debe utilizarse despu´es de plot. Es posible realizar distintas representaciones en un mismo gr´afico. Por ejemplo, x=0:.01:2*pi;y1=sin(x); y2=sin(2*x); y3=sin(4*x); plot(x,y1,x,y2,x,y3)
Otra posibilidad es utilizar hold. El comando hold on “congela” la terminal gr´afica en la que estamos trabajando, de modo que se pueden superponer diversos gr´aficos en ella. Los ejes pueden, sin embargo, ser reescalados. El comando hold off “descongela” la terminal gr´afica. Dentro de las opciones gr´aficas podemos elegir el tipo de l´ınea, el tipo de punto y el color. Por ejemplo, x=0:.01:2*pi;y1=sin(x); y2=sin(2*x); y3=sin(4*x); plot(x,y1,’--’,x,y2,’:’,x,y3,’+’)
37 dibuja una l´ınea discontinua y punteada, respectivamente, para los dos primeros gr´aficos mientras que el tercer gr´afico se muestra con el s´ımbolo “+”. Los tipos de l´ınea y marca son: Tipos de l´ınea: s´olida (-), discontinua (–), punteada (:), discontinua y punteada (-.) Tipos de marca: punto (.), mas (+), estrella (*), c´ırculo (o), x (x) Se pueden especificar colores para los distintos tipos de l´ınea y marca: Colores: amarillo (y), magenta (m), rojo (r), verde (g), azul (b), blanco (w), negro (k) El comando subplot puede utilizarse para hacer una partici´on de la terminal gr´afica, de modo que pueden situarse varios subgr´aficos en una misma figura.
1.11.2.
Representaciones en 3D
Gr´ aficos de l´ınea El comando plot3 en 3 dimensiones es el an´alogo al comando plot en 2 dimensiones: produce curvas en el espacio tridimensional. Si x, y y z son vectores del mismo tama˜ no, entonces el comando plot3(x,y,z) producir´a un gr´afico de perspectiva de la curva en el espacio tridimensional que pasa por los puntos especificados por x, y y z . Estos vectores suelen estar definidos de forma param´etrica. Por ejemplo, t=.01:.01:2*pi; x=cos(t); y=sin(t); z=t.^3; plot3(x,y,z)
proporciona una h´elice que est´a comprimida cerca del plano x-y.
Gr´ aficos de malla y de superficie Pueden obtenerse gr´aficos tridimensionales mallados de superficies utilizando el comando mesh. Si tecleamos mesh(z) obtenemos un gr´afico de perspectiva tridimensional de los elementos de la matriz z . La superficie representada est´a definida por las coordenadas z de los puntos sobre un ret´ıculo rectangular en el plano x-y. El siguiente ejemplo muestra un gr´afico de estas caracter´ısticas de los elementos de la matriz identidad 10 10 (comando eye(10)):
×
mesh(eye(10))
An´alogamente pueden obtener gr´aficos tridimensionales “compactos” de superficies utilizando el comando surf :
38 surf(eye(10))
Para dibujar el gr´afico de una funci´on z = f (x, y) sobre un rect´angulo debemos, en primer lugar, definir vectores xx e yy que proporcionan particiones sobre los lados del rect´angulo. Con la funci´on meshgrid creamos una matriz x, cada fila de la misma contiene las componentes del vector xx y cuyo n´umero de columnas es igual a la longitud del vector yy. De manera an´aloga creamos la matriz y, cuyas columnas contienen las componentes del vector yy. Esto lo conseguimos con la instrucci´on: [x,y]=meshgrid(xx,yy);
Una vez hecho esto, obtenemos la matriz z haciendo actuar f sobre las matrices x e y. La representaci´on de la matriz z se puede hacer acudiendo a los comandos mesh y surf . Veamos un ejemplo: Vamos a dibujar el gr´afico de z = e x y sobre el cuadrado [ 2, 2] [ 2, 2] del siguiente modo: −
2
−
2
−
×−
xx=-2:.2:2; yy=xx; [x,y]=meshgrid(xx,yy); z=exp(-x.^2-y.^2); mesh(z)
Las caracter´ısticas del comando axis introducido previamente son aplicables tambi´en a los gr´aficos tridimensionales, as´ı como lo son las de los comandos para t´ıtulos, etiquetado de ejes y el comando hold. El color de las superficies se ajusta utilizando el comando shading. Hay 3 opciones para este comando: faceted (el que est´a por defecto), interpolated y flat. Se accede a estas opciones tecleando: shading faceted, shading interp, shading flat
El comando shading debe introducirse despu´es del comando surf . La utilizaci´on de shading interp y shading flat causa la desparici´on del mallado en la superficie. El perfil de colores de una superficie se controla mediante el comando colormap. Mapas de colores predefinidos son: hsv (por defecto), hot, cool, jet, pink, copper, flag, gray, bone
Por ejemplo, el comando colormap(cool) proporciona un determinado perfil de colores en la figura en cuesti´on.
39
1.12 1.12..
Resum Resumen en de func funcio ione ness elem elemen enta tale less y matric matrices es especiales
En la tabla que se muestra a continuaci´on on se presentan algunas funciones de MATLAB MATLAB que nos pueden ser de utilidad:
abs sqrt real imag conj exp log log10 sin, asin, sinh, asinh cos, acos, cosh, acosh tan, atan, tanh, atanh cot, acot, coth, acoth sec, asec, sech, asech csc, acsc, csch, acsch
valor absoluto absoluto o m´ odulo odulo ra´ız ız cuadra cua drada da parte real parte imaginaria complejo conjugado exponencial logaritmo natural logaritmo logaritmo en base 10 seno, arcseno, arcseno, seno hiperb., hiperb., arcseno arcseno hiperb. idem para el coseno coseno idem para para la tangen tangente te idem para la cotangen cotangente te idem para la secante secante idem para la cosecan cosecante te
Finalmente, algunas matrices especiales de MATLAB:
zeros ones eye diag rand
matriz matriz de ceros matriz matriz de unos identidad diagonal n´umeros umeros aleatorios distribuidos unif.
40
Parte II Introducci´ on on a la construcci´ on on de algoritmos
41
Cap´ıtulo 2 Sistemas de numeraci´ on; errores y sus fuentes En este tema se introducen los distintos sistemas de numeraci´on y se describe c´omo realizar la conversi´on entre distintos sistemas, as´ı como el est´andar IEEE de representaci´ on de n´ umeros en el sistema binario. A continuaci´on, se introducen las nociones de error absoluto y relativo y se analizan las causas m´as frecuentes de error en un c´alculo computacional, ilustr´andolo con ejemplos. Finalmente, se introducen los conceptos de (in)estabilidad y condici´on.
2.1.
Sistemas de n´ umeros y conversiones
Es necesario conocer como se representan los n´umeros enteros y decimales en formato digital para entender las limitaciones intr´ınsecas que nos encontraremos al programar un algoritmo num´erico, tanto en cuanto a la mayor precisi´on alcanzable como en cuanto al rango de valores admisibles. Veremos adem´a s c´omo un m´ınimo conocimiento del sistema de representaci´on est´andar sirve para evitar la aparici´on de desastrosos errores de programaci´ on (como bucles infinitos).
2.1.1.
Sistemas de numeraci´ on en base q
El sistema de numeraci´on decimal es el sistema posicional de numeraci´o n m´as utilizado, que utiliza como base de numeraci´on el 10 (d´ıgitos 0...9). Como es bien sabido, se puede utilizar cualquier n´umero natural mayor que 1 como base de numeraci´on. 43
44 Un n´ umero en base q se denota como (an an 1 ...a1 a0 .b1 b2 ...bk ...)q donde ai y b j pertenecen al conjunto de los q d´ıgitos elementales 1 . Estos q d´ıgitos representar´an valores desde 0 hasta q 1. La conversi´on a decimales es, por definici´on: −
−
(an an 1 ...a1 a0 .b1 b2 ...bk ...)q = an q n + an 1 q n 1 + ... + a1 q + a0 q 0 +b1 q 1 + b2 q 2 + ... + bk q k + ... −
−
−
−
−
−
(2.1)
El sistema natural de numeraci´on digital es el binario (base 2), utilizando s´olo los d´ıgitos 0 y 1.
Ejemplo 2.1.1. Decimal: (123.25)10 = 1
2
× 10
+2
1
× 10
+3
0
× 10
+2
× 10
−1
+5
× 10
−2
.
Binario (base 2) con 2 d´ıgitos 0 y 1: 0 + 1 = 1, 1 + 1 = 10.
(1011.01)2 = 1
3
×2
+0
×2
2
+1
1
×2
0
+1
×2
+0
−1
×2
+1
−2
×2
= 11.25.
Hexadecimal (base 16) tiene 16 d´ıgitos: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E y F , que representan los valores desde 0 hasta 15 respectivamente. Por tanto: (7B.4)15 = 7
1
× 16
+ B
0
× 16
+4
−1
× 16
= 123.25.
Los sistemas de numeraci´on octal y hexadecimal, tambi´ en son muy utilizados, en parte debido a que es muy sencilla su conversi´on a binario con la ventaja de que un mismo n´ umero se representa con menos cifras cuanto mayor es la base.
2.1.2.
Conversi´ on de base decimal a base q
Ya sabemos como convertir un n´umero en base q a base decimal (2.1). La conversi´on de base decimal a base q se base en el hecho de que, acudiendo a la definici´on (2.1) observamos que (eliminamos el sub´ındice 10 para n´ umeros en base decimal): (an an 1 ...a1 a0 .b1 b2 ...bk ...)q (an an 1 ...a1 a0 .b1 b2 ...bk ...)q − −
1
× q = (a a ...a a b .b b ...b ...) × q = (a a ...a .a b b ...b ...) n n−1
−1
n n−1
1
2 3
k
q
0 1 2
k
q
0 1 1
(2.2)
Vamos a utilizar el punto decimal ( .) en lugar de la coma decimal ( ) para ser consistentes con los ejemplos MATLAB. ′
45 Utilizaremos esto para obtener la conversi´on de un n´ umero en base 10 a base 2. Para esto, conviene separar la conversi´on de la parte entera (antes del punto decimal) de la fraccionaria (despu´es del punto); observemos que un n´umero entero es necesariamente entero en cualquier base y que un n´umero fraccionario lo es en cualquier base.
Conversi´ on entera De (2.2) deducimos que: −1
(an ...a0 )q
= (an ...a1 )q + (.a0 )q
es decir, que (an ...a0 )q = (an ...a1 )q
× q + (.a ) × q = (a ...a ) × q + (a ) 0
q
n
1
q
0
q
(2.3)
que es una identidad entre n´umeros, que podemos evaluar en cualquier base, y en particular en base 10 (obs´ervese que de hecho mezclamos base n´ umeros en distintas bases). Vemos pues que al dividir un n´ umero entero entre q el resto es el d´ıgito menos significativo del n´ umero en base q. Dividiendo sucesivamente (hasta llegar al cociente 0) obtenemos los sucesivos d´ıgitos del n´umero en base q .
Conversi´ on fraccionaria A partir de (2.2) vemos que: (.b1 b2 ...bk )q
× q = (b ) + (.b ...b ) 1
q
2
k q
Vemos pues que la parte entera del resultado de multiplicar nuestro n´umero por q es el primer d´ıgito tras el punto decimal. Reteniendo la parte fraccionaria que resulta en cada paso y repitiendo sucesivamente el proceso, vamos obteniendo el resto de cifras tras el punto.
Ejemplo 2.1.2. Escribamos (26.1)10 en base 2. 1. Parte entera. Dividiendo sucesivamente, tenemos que: 26 = 2
× 13 + 0 ; 13 = 2 × 6 + 1 ; 6 = 2 × 3 + 0 ; 3 = 2 × 1 + 1 ; 1 = 2 × 0 + 1
Leyendo de izquierda a derecha los n´ umeros subrayados: (26)10 = (11010)2
46 2. Parte fraccionaria. Multiplicando sucesivamente por dos y separando la parte fraccionaria: 0.1 2 = 0.2 ; 0.2 2 = 0.4 ; 0.4 2 = 0.8 ; 0.8 0.2 = 1.6 ; 0.6 2 = 1.2 ; 0.2 2 = ...
×
×
×
×
×
×
Leyendo de izquierda a derecha tenemos los d´ıgitos de la parte fraccionaria (subrayados) luego: (0.1)10 = (0.00011)2 donde las cifras subrayadas son las cifras peri´ odicas. Obs´ervese que el n´ umero 0.1 tiene infinitas cifras distintas de cero cuando se escribe en base 2 (se repiten indefinidamente desde la segunda a la quinta cifra fraccionaria). Sumando los resultados de la parte entera y la fraccionaria tenemos que (26.1)10 = (11010.00011)2 En el anterior ejemplo hemos aprendido que un n´umero fraccionario con un n´umero finito de d´ıgitos en cierta base no tiene por qu´e tener un n´umero finito en otra base. En efecto, (0.1)10 es un n´umero peri´odico con infinitas cifras cuando se escribe en base 2. Esto es importante puesto que en un procesador digital los n´umeros se almacenar en binario utilizando un n´ umero finito de cifras. Esto quiere decir que el n´umero 0.1 no se representa exactamente en un ordenador. Esto hay que tenerlo en cuenta para evitar errores que conduzcan a bucles que se repiten indefinidamente. El siguiente algoritmo en Matlab es un ejemplo de bucle infinito, on finita binaria: Algoritmo 2.1. Bucle infinito en precisi´ x=0; while x =10 (Nota: en Matlab esto significa x = 0) x=x+0.1; end
∼
Esto nos deber´ıa prevenir de utilizar bucles que terminan cuando se alcanza cierto valor concreto de una variable que involucra n´umeros no enteros.
Ejercicio 2.2. Si en la tercera l´ınea del anterior algoritmo cambiamos 0.1 por 0.125, ¿seguir´ıamos teniendo un bucle infinito?.
47
2.1.3.
Est´ andar IEEE de representaci´ o n de n´ umeros enteros y en coma flotante
N´ umeros enteros Asumamos, como es com´ un en la mayor parte de los procesadores, que cada palabra tiene una longitud de 32 bits, es decir, que cada n´umero vendr´a representado por 32 d´ıgitos binarios. La forma m´ as sencilla de representar un n´umero entero es asignando uno de los 32 bits para designar el signo (0 para +, 1 para menos) y utilizando las 31 restantes posiciones para representar los d´ıgitos del m´odulo del n´ umero entero. En consecuencia, el 31 mayor n´ umero entero que se puede representar es 2 1. ´ Este es el est´andar IEEE salvo porque, por lo general, no se suele reservar un bit para el signo, sino que se utiliza una ordenaci´on distinta para guardar los enteros negativos. Sin embargo, en la pr´actica una y otra prescripci´on son pr´acticamente equivalentes. Tambi´ en pueden utilizarse enteros en longitud doble, utilizando dos palabras de 32 bits, con lo que se aumenta el rango de enteros admisibles.
−
N´ umeros no enteros Para representar n´umeros con parte fraccionaria se utiliza el formato de punto (o coma) flotante. Esta representaci´on es la correspondiente versi´on binaria de la conocida notaci´ on cient´ıfica o exponencial para los n´ umeros decimales: E
±M × 10
, 1
≤ M < 10
pero utilizando base 2: x =
E
±M × 2
, 1
≤ M < 2
donde M es la mantisa y E el exponente. Inevitablemente el n´umero de d´ıgitos que se pueden almacenar de la mantisa y exponente es limitado. Hay dos tipos de precisi´on, la simple y la doble, que difieren en el n´umero de bits de los que se dispone para almacenar las cifras; la distribuci´on est´andar es: 1. Precisi´on simple: 32 bits, de los cuales: a ) uno corresponde al signo, b) 8 al exponente, luego c ) 23 a la mantisa
8
−128 ≤ E ≤ 127 (2 × 128 = 2 )
48 2. Precisi´on doble a ) uno para el signo b) 11 al exponente ( 1022
−
c ) 52 para la mantisa
≤ E ≤ 1023)
Ejemplo 2.1.3. Veamos de forma simplificada como se guardar´ıa el n´ umero 5.5 en precisi´ on simple. Primero lo pasamos a binario: 5.5 = (101.1)2 y normalizamos para que la mantisa 1 < M 2, de forma que 5.5 = (1.011)2 22 , que tiene exponente E = 2 (que se guardar´ a en binario), signo + (bit 0) y mantisa 1.011. Normalmente, salvo para el caso del n´ umero cero, que se almacenan de distinta forma, el n´ umero delante del punto ser´ a siempre 1, por lo que no se suele almacenar.
≤
×
Ejemplo 2.1.4. El anterior es un ejemplo de un n´ umero decimal que se puede guardar de forma exacta en una palabra de 32 bits. Un ejemplo en el que esto no ser´ıa as´ı es el ya conocido caso de 0.1 = (0.00011)2 = (1.1001)2 2 4 , que necesariamente ha de redondearse antes de almacenarse; esquem´ aticamente (sin entrar en detalles de como se almacena el exponente), tendr´ıamos, en precisi´ on simple, la representaci´ on (recordemos que el 1 antes del punto decimal no se suele almacenar)
×
0 E =
|
−
−4|10011001100110011001101
donde la 23a cifra tras el punto se ha redondeado a 1 porque la siguiente era 1. De esta forma, lo que realmente se almacena es: (1.10011001100110011001101)2
4
×2
= 0.10000000149011611938
Esto muestra expl´ıcitamente que, efectivamente, el algoritmo 2.1 es un bucle infinito. Debemos resaltar que la especificaci´o n del n´ umero de bits que se utilizan para almacenar mantisa y exponente es lo ´unico que necesitamos para determinar cuales son los mayores y menores n´ umeros que se pueden almacenar en coma flotante as´ı como la m´axima precisi´on que se puede obtener. Veamos como obtener estos par´ametros en el caso de doble precisi´on (que es la precisi´on utilizada por Matlab). umeros de “overflow” y “underflow” as´ı como la precisi´ on Ejemplo 2.1.5. Obtener los n´ m´ aquina en doble precisi´ on, es decir:
49 1. Obtener el mayor n´ umero entero positivo representable en formato de coma flotante en el caso de doble precisi´ on (l´ımite de “overflow”) Puesto que el mayor exponente es 1023, el mayor n´ umero representable es (1.1....1)2 1023 308 2 1.8 10 (hay 52 unos detr´ as del punto “decimal”).
≃
×
×
2. Obtener el menor n´ umero entero positivo representable en formato de coma flotante en el caso de doble precisi´ on (l´ımite de “overflow”). El menor decimal representable con 52 d´ıgitos significativos binarios es (1.00...01)2 10 1022 2.23 10 308 . Sin embargo, el est´ andar IEEE admites n´ umeros m´ as peque˜ nos, aunque tengan menos cifras significativas (es lo que llaman n´ umeros sub-normales); estrictamente el n´ umero 1022 1074 m´ as peque˜ no representable en coma flotante es: (0.00..01)2 2 = 2 324 4.94 10 (52 d´ıgitos tras la coma decimal). Para no perder d´ıgitos significativos es mejor moverse en el rango (1/N ov , N ov ), donde N ov es el n´ umero de overflow 308 . N ov 2.23 10
×
× ≃
−
≃
×
−
×
×
−
−
−
≃
−
3. Obtener la diferencia entre 1 y el menor n´ umero positivo mayor que 1 en coma flotante (´ epsilon-m´ aquina). El n´ umero que se nos pide es (1.0...01)2 (1.0...0)2 = (0.0...01)2 = 2 52 2.2 16 10 . Este n´ umero representa el mejor error relativo que se puede manejar en doble precisi´ on (definimos este conocido concepto a continuaci´ on). Este valor de ´epsilon-maquina significa que todos los n´ umeros con 15 cifras significativas se pueden guardar de forma exacta en doble precisi´ on y que ocurre lo mismo para la mayor parte de los n´ umeros de 16 cifras. −
−
−
≃
×
on 2 Nota 2.3. Es importante tener en cuenta que MATLAB trabaja en doble precisi´ aunque por defecto s´ olo muestre 5 cifras significativas. Para que el sistema muestre m´ as cifras se puede utilizar el comando “format long”. Adem´as de los detalles se˜ nalados hay otras caracter´ısticas fundamentales del est´andar IEEE que no describiremos como son: 1. Redondeo correcto de la aritm´etica (volveremos m´as adelante sobre este punto)3 . 2
La versi´on 7.0 ya permite trabajar en precisi´on simple y con aritm´etica entera, aunque doble precisi´on sigue siendo la disponible por defecto. 3 Un error en el “hardware” de punto flotante de los Intel Pentium (1994) le dio una considerable mala publicidad a la compa˜n´ıa. El problema fue solventado reemplazando los procesadores defectuosos
50 2. Tratamiento de excepciones. Es decir, que un operaci´on del tipo log(0.0) no da error e interrumpe la ejecuci´on del programa, sino que se asigna un valor especial para representar esta excepci´on (para significar ).
−∞
3. Compatibilidad entre procesadores. Esta es por supuesto, una de las grandes ventajas de los est´andares.
2.2.
Errores y sus causas
El an´alisis num´erico trata de la construcci´on de m´etodos discretos para la resoluci´on de problemas continuos. Esta discretizaci´on, que se presenta tanto en la representaci´on num´erica en un ordenador como en los propios algoritmos de c´ alculo, necesariamente implica la aparici´on de errores cuyo origen y propagaci´on debemos estudiar 4 .
2.2.1.
Definiciones
Si x A es una aproximaci´on del verdadero valor x T , definimos entonces: Error absoluto: E abs = xT
| − x |
| − xx |,
Error relativo: E rel = 1
A
T
A
si xT = 0.
Tambi´en se pueden utilizar estas definiciones con signo (es decir, sin el valor absoluto). El error relativo en ocasiones se expresa tambi´ en en tanto por ciento.
2.2.2.
Fuentes de errores
Las fuentes de error en un c´alculo num´erico pueden ser de variada naturaleza, algunas de ellas independientes del m´etodo num´erico empleado. 1. Un m´etodos num´erico para resolver determinado problema cient´ıfico puede dar lugar a resultados err´oneos si el modelo no es una fiel descripci´on de la realidad. 4
No por ello se debe adoptar la visi´on reduccionista consistente en definir el an´alisis num´erico como el ´area de la matem´atica que analiza los errores de c´alculo
51 2. Un algoritmo num´erico puede depender de ciertos datos de entrada (por ejemplo, datos f´ısicos) afectados de cierto error. Se debe vigilar que el m´etodo num´erico no sea crucialmente dependiente de la exactitud de tales valores de entrada y que la precisi´ on de estos valores sea suficiente para nuestros prop´ositos. En el caso en que el propio modelo f´ısico/matem´atico sea muy sensible a estos par´ametros, tendremos un problema mal condicionado e imposible de resolver num´ ericamente de forma estable. 3. Los errores de programaci´on son pr´acticamente inevitables durante la construcci´on de un m´etodo num´erico. Una vez construido un algoritmo num´erico que en apariencia funciona correctamente es necesario certificar su funcionamiento mediante todos los tests que est´en a nuestro alcance. 4. Errores en la aritm´etica de punto flotante: errores de redondeo, p´erdida de cifras significativas por cancelaciones, problemas de “underflow/overflow”. Estos son los problemas derivados de la representaci´on en punto flotante, en la cual se dispone de un n´ umero finito de bits para representar la mantisa y el exponente. Como ya vimos, esto limita tanto la mejor precisi´on relativa alcanzable (15-16 d´ıgitos decimales en doble precisi´on) como el rango de valores disponible. Aunque la segunda de las limitaciones no suele presentar graves problemas, es necesario evitar la evaluaci´on de cantidades demasiado grandes o peque˜ nas. La limitaci´on de la precisi´on tiene por lo general mayores consecuencias. Por ejemplo, cuando se sustraen cantidades muy parecidas, el error relativo empeora dr´asticamente por p´erdida de cifras significativas. 5. Errores de truncamiento o de discretizaci´on, inherentes al hecho de aproximar un problema continuo mediante una aproximaci´on discreta. Por ejemplo, veremos que la regla trapezoidal sirve para aproximar integrales mediante: N −1
b
a
f (x)dx
≈
h (f (a)+f (b))+h f (a+kh) 2 k=1
≡ S (h) , donde h = b −N a , h = (b−1)/N
El significado de “ ” es “aproximadamente”, lo cual significa que
≈
b
a
f (x)dx = S (h) + ǫ(h)
52 donde ǫ(h) es el error de truncamiento, que es de esperar que sea menor cuanto menor sea h (N N lo mayor posible). Un c´alculo expl´ıcito de los errores de truncamiento es al menos igual de dif´ıcil que el c´alculo del problema original. Nos conformaremos con un conocimiento cualitativo de estos errores y con buscar acotaciones lo m´as finas posible.
∈
Trataremos de estimar los errores de truncamiento en cada algoritmo num´ erico que se presente. Discutamos en este punto acerca de los errores debidos a las limitaciones en la representaci´on de punto flotante, en particular por lo que respecta a la precisi´on limitada del sistema
Redondeo y p´ erdida de cifras significativas Puesto que la representaci´on en coma flotante es limitada en cuanto al n´umero de cifras significativas que se pueden almacenar, los n´umeros reales se redondear´an a un determinado n´ umero de cifras (siempre utilizando el sistema de numeraci´on binario) cuando el valor verdadero tenga m´as cifras de las que se pueden almacenar. Asimismo, tras cada operaci´on aritm´etica se vuelve a redondear el resultado. En el est´andar IEEE esto se hace de la mejor manera posible en el sentido de que el resultado de la operaci´on aritm´etica est´a redondeado de forma que se almacena el valor m´as pr´oximo al resultado exacto. As´ı, por ejemplo, supongamos que, para simplificar, la precisi´on fuese de tres d´ıgitos binarios tras el punto (y no consideremos limitaci´ on en el exponente), si se suman los n´umeros en punto flotante x = (1.010)2 , y = (1.001)2 2 4 tenemos que el resultado exacto es (1.0101001)2 que no ser´ıa un n´umero en punto flotante (s´olo disponemos de 3 posiciones tras el punto). El resultado IEEE ser´ıa redondear la u´ltima cifra, aument´andola si la siguiente fuese 1 y dej´andola como est´a en otro caso. As´ı, x + y se almacenar´ıa como (1.011)2 . Otra posibilidad (menos recomendable) es, sin m´as, eliminar las cifras que no quepan (truncamiento). Esta desafortunada prescripci´on se utilizaba en los supercomputadores Cray, pero est´a en total deshuso. A´un cuando en el est´andar IEEE se redondea al n´ umero en coma flotante m´as cercano, que es la opci´on m´as conveniente, es necesario estudiar como pueden afectar sucesivos redondeos al resultado de un algoritmo num´ erico. Por otra parte, la limitaci´o n en el n´ umero de bits para la mantisa tiene como consecuencia la posible p´erdida de cifras significativas en ciertos c´alculos donde dos cantidades
×
−
53 tienden a cancelarse entre s´ı 5 . Para evitar este tipo de errores, es recomendable:
a) O bien reescribir la f´ormula en cuesti´on de modo que se eviten las restas de cantidades de la misma magnitud. b) O bien utilizar (cuando sea posible) un desarrollo de Taylor para aproximar la f´ormula hasta la precisi´on requerida. Veamos algunos ejemplos en los que se da p´erdida de cifras significativas
Ejemplo 2.2.1. Obtener las ra´ıces de x2 106 x + 1. En una modesta calculadora con 10 decimales y aplicando la archiconocida f´ ormula 2 4ac , vamos obteniendo los resultados: b x = b 2a
− ±
√
−
−
x =
106
±
√
1012 2
6
6
− 4 ≃ 10 ± 10 2
que da x = 106 para la mayor ra´ız y x = 0 para la menor. Es evidente que hemos perdido todas las cifras significativas para la menor ra´ız. De hecho, la famosa f´ ormula deber´ıa ser desterrada de cualquier m´ etodo num´ erico. Es mucho mejor reescribir la soluci´ on de la ecuaci´ on de segundo grado como: x1 = L/2a , x2 = 2c/L, L =
√
−b − signo(b)
b2
− 4ac
De esta forma la mayor ra´ız es como antes y en la calculadora la menor saldr´ıa x = 10 6 , cuyo error relativo respecto a la soluci´ on verdadera es 10 12 . −
∼ √ √ Ejemplo 2.2.2. Consideremos por ejemplo f (x) = ( x + 1 − x) para x > 0. Cuando x √ es √ grande, se producen errores de redondeo importantes puesto que los valores de x + 1 y x se aproximan considerablemente. ¿C´ omo evitar´ıamos esta fuente de error?. Lo que −
podemos hacer es reescribir f (x) del siguiente modo:
√ x + 1 + √ x √ √ 1 = √ f (x) = ( x + 1 − x) √ √ √ . x + 1 + x x + 1 + x En la expresi´ on resultante podemos apreciar que no se producen cancelaciones de cantidades de la misma magnitud. 5
Un famoso y lamentable error de ´este tipo fue el que motiv´o le fallo de los misiles Patriot para derribar los misiles Skud iraqu´ıes durante la guerra del golfo: se med´ıan intervalos de tiempo restando tiempos desde el reinicio del sistema con la consiguiente p´ erdida progresiva de precisi´on
54 1 cos(x) . Cuando se eval´ ua g(x) para x << 1 Ejemplo 2.2.3. Consideremos g(x) = x2 se produce una p´erdida considerable de cifras significativas. En este caso, para remediar el problema consideramos el desarrollo de Taylor de cos(x) entorno a 0. De este modo:
−
g(x) =
1
2
| |
4
x 4 + . 4! 6!
− [1 − x /2! + x /4! + ...] ≈ 1 − x
2
2 x2 Esta expresi´ on es apropiada entonces para evaluar g(x) cuando x << 1.
Errores de overflow /underflow Aunque los l´ımites de overflow/underflow son considerablemente generosos, conviene escribir las expresiones que se utilicen de manera que se minimice esta posibilidad. Por ejemplo, es preferible calcular el m´o dulo de un n´ umero complejo, z = x + iy, z = x2 + y 2 , como z = x 1 + (y/x)2
||
|| ||
De esta forma, no hay que elevar al cuadrado n´umeros que pueden ser muy grandes. De la misma forma, hay que tomar precauciones para calcular la argumento de un n´umero complejo de manera que no se produzcan “overflows/underflows” en el c´alculo de y/x.
Ejercicio 2.4. Escribir un programa en Matlab que obtenga la representaci´ on polar de iθ un n´ umero complejo dado z = re , 0 θ < 2π, eliminando el riesgo de problemas de overflow-underflow.
≤
2.3.
Propagaci´ on de errores: condici´ on y estabilidad.
Un error en un c´alculo num´ erico contamina las sucesivas evaluaciones. Esta propagaci´on del error puede describirse en t´ erminos de dos conceptos relacionados, los de (in)estabilidad y condici´on. La condici´ o n de una funci´on f (x) mide la sensibilidad de los valores de f (x) a peque˜ nos cambios en x y se define como
E (f (x)) C = E (x) rel
rel
donde E rel (f (x)) es el error relativo de f (x) para un error relativo E rel (x) en x.
55 Entonces, como
6
′
f (xT )
′
− f (x ) ≃ f (x )(x − x ) → E A
t
luego
T
rel (f (x))
A
≃
f (xT ) (xT f (xT )
A)
− x
f (x ) C ≃ x f (x ) ′
T
T
T
Utilizaremos esta u ´ ltima expresi´on como definici´on de condici´on para funciones f (x) de una variable real. Definimos entonces los n´umeros de condici´ on como
f (x) C (x) = x f (x) ′
Cuando para un x dado 0 < C (x) < 1 para ese x se dir´a que el problema (c´alculo de f) est´a bien condicionado (y cuanto menor sea C mejor condicionado), mientras que si C (x) > 1 el problema estar´a mal condicionado. Si C (X ) = 1, el error relativo se mantiene.
√
Ejemplo 2.3.1. La funci´ on f (x) = x est´ a bien condicionada, pues C (x) = 1/2, luego el error relativo se reduce. 2 En cambio f (x) = x 2 1 est´ a mal condicionada para x 1 pues C (x) = 22x 1 x
−
≃
−
El concepto de condicionamiento se puede extender a situaciones m´as generales que la de una funci´on de una variable continua. Por ejemplo, un problema cl´asico que involucra funciones de una variable discreta, es el estudio del condicionamiento de relaciones de recurrencia: on de recurrencia: Ejemplo 2.3.2. Las funciones de Bessel J n (x) satisfacen la relaci´ 2n J n+1 (x) = J n 1 (x) + x J n (x), pero como las funciones de Bessel de segunda especia cumplen la misma relaci´ on y l´ımn alculo de las funciones J n a J n (x)/Y n (x) = 0, el c´ partir de J 0 y J 1 est´ a mal condicionado, pues una peque˜ na perturbaci´ on en los datos iniciales J 0 y J 1 contamina nuestra secuencia de funciones J n con la secuencia Y n , que crece m´ as r´ apido con n. As´ı, por ejemplo, empezando con los valores en precisi´ on simple J 0 (2) = 0.22389078 y J 1 (2) = 0.57672481, y aplicando la recurrencia obtenemos J 8 (2) = 4.00543213 10 5 que no tiene bien ni siquiera un cifra significativa: J 8 (2) = 2.2180 10 5 .
−
−
→∞
{ }
{ }
×
6
−
×
−
Recordemos el teorema del valor medio: Si g (x) continua en [ a, b] y derivable en ( a, b) entonces
∃c ∈ (a, b) : f (b) − f (a) = f (c)(b − a) ′
56 Un concepto relacionado, que no equivalente, es el de (in)estabilidad de un algoritmo, que describe la sensibilidad de un m´etodo num´erico espec´ıfico respecto a los inevitables errores de redondeo cometidos durante su ejecuci´on en aritm´ etica de precisi´on finita. Observemos que la condici´on no depende de errores de redondeo pero que la estabilidad de un algoritmo s´ı depende del condicionamiento de la funci´ on que queramos evaluar. Un ejemplo puede servir para distinguir estos dos conceptos relacionados:
Ejemplo 2.3.3. Dada la funci´ on f (x) =
√ x + 1 − √ x, su n´ umero de condici´ on es:
f (x) x = √ √ C (x) = x f (x) 2 x x + 1 ′
y vemos que C (x) < 1/2 para x > 0, luego la funci´ on est´ a bien condicionada (su error relativo es menor que el error relativo en x). Sin embargo, el algor´ıtmo para calcular x consistente en ir realizando las operaciones implicadas en f (x) = x + 1 x, a saber:
√
− √
1. Input: x 2. y = x + 1
√ x + 1 √ 4. f = x 3. f 1 = 2
5. f = f 1
− f
2
es inestable para x grande por culpa del paso 5 (hay cancelaciones entre n´umeros similares). Como sabemos, un algoritmo estable lo proporciona la siguiente reescritura de la funci´on: 1 f (x) = x + 1 + x
√
2.4.
√
Eficiencia de un algor´ıtmo num´ erico
Por supuesto, cualquier algoritmo num´erico sensato debe evitar ser inestable. Por otra parte, si existieran varios m´etodos para evaluar una misma funci´on, entonces convendr´a adoptar aquel m´etodo que sea m´ as eficiente, es decir, m´as r´apido. Con la mejora en proporci´on geom´etrica de la velocidad de los procesadores, podr´ıamos estar tentados en
57 despreocuparnos por la rapidez de c´alculo. Esta es, sin embargo, una p´esima filosof´ıa: se trata de aprovechar los recursos para poder resolver problemas m´as complejos y no para resolver peor problemas simples. La eficiencia es y siempre ser´a de importancia capital en el desarrollo de buenos m´etodos num´ericos. La diferencia en tiempos de ejecuci´on pueden llegar a ser muy considerables si ciertas operaciones elementales, que se pueden repetir miles y miles de veces, no se realizan con cuidado. Por ejemplo, para calcular x 4 para cierto valor de x, es muy mala idea calcular x4.0 (exponente en coma flotante); hay que tener cuidado en utilizar un exponente entero ya que los m´etododos de exponenciaci´on en coma flotante son distintos que los de enteros y mucho m´as lentos; a´ un es mejor idea considerar el c´alculo en dos pasos: x2 = x x, x4 = x 2 x2 , con lo que se economizar un producto frente a x 4 = x x x x.
∗
2.4.1.
∗ ∗ ∗
∗
Ejemplo: Evaluaci´ on de polinomios, m´ etodo de Horner
Otro ejemplo notable (y por desgracia no suficientemente) conocido es la evaluaci´on de polinomios. Por ejemplo, supongamos que nos planteamos evaluar el polinomio P (x) = 2 + 4x
2
− 5x
+ 2x3
4
− 6x
+ 8x5 + 10x6
Contando con que cada potencia de exponente k entero cuente como k 1 productos, tendr´ıamos que el n´umero total de productos para evaluar el polinomio de forma directa es:
−
1 + 2 + 3 + 4 + 5 + 6 = 21 y el n´umero de sumas es 6. Una forma mejor de proceder es ir calculando primero las potencias de forma sucesiva: x2 = x x , x3 = x x2 , x4 = x x3 , x5 = x x4 , x6 = x x5 con lo que s´olo se a˜ nade una nueva multiplicaci´on por potencia, par un total de 1 + 2 + 2 + 2 + 2 + 2 = 11 Pero a´ un se puede hacer mejor reescribiendo P (x) = 2 + x(4 + x( 5 + x(2 + x( 6 + x(8 + x10)))))
−
−
con lo que s´olo necesitamos 6 multiplicaciones (y el n´ umero de sumas no cambia). Vemos pues que para evaluar un polinomio de grado n en el que ninguno de los coeficientes
58 es cero, se necesitan n(n + 1)/2 multiplicaciones por el primer m´etodo, 2n 1 por el segundo y n para el tercero. De forma que se debe procurar utilizar el ´ultimo m´etodo, particularmente cuando n es grande. Este m´etodo es adem´as sencillo de programa, en efecto:
−
Algoritmo 2.5 (Algoritmo de Horner o de divisi´on sint´etica). Algoritmo de Horner Dado el polinomio P (x) = a 0 + a1 x + ... + an xn , an = 0
la evaluaci´ on de P (x) para cierto valor x = z se puede realizar en n pasos mediante
(1) bn = an (2) bn 1 = a n 1 + z bn −
−
(3) bn 2 = a n −
...
−2
∗ + z ∗ b
n−1
(n) b0 = a 0 + z b1
∗
donde P (z ) = b0 . Es f´acil programa este algoritmo en forma de bucle, sobretodo si no nos interesan los c´alculos intermedios. As´ı, se puede escribir:
(1) b = a n (2) Repetir mientras n > 0 (3) n = n
−1
(4) b = a n + z b
∗
(5) Volver a (2) (6) p(z ) = b.
59 Esta forma de evaluar polinomios es mucho mejor que el m´ etodo directo, especialmente para ´ordenes grandes. Dependiendo del ´orden de un polinomio y las veces que se repita el c´alculo, ser´a importante aplicar el m´etodo de Horner.
60
2.5.
´ EJERCICIOS PRACTICOS Introducci´ on a la construcci´ on de algoritmos
El objetivo de estos ejercicios es el de poner en pr´actica algunos conceptos introducidos en este tema desarrollando algunos ejemplos con MATLAB. Ejercicios: Se propone realizar los siguientes ejercicios y responder a las cuestiones que se plantean: −52
1. Comprobar que el epsilon-m´aquina es 2 comandos:
−16
= 2.220410
, tecleando en la l´ınea de
>> a=1+2^(-53);b=a-1
y comparando con >> a=1+2^(-52);b=a-1
2. Escribir la secuencia de comandos: x=0; while x~=10 x=x+0.1 end
en un fichero (con extensi´on .m) y ejecutarlo en MATLAB. Para interrumpir la ejecuci´on, pulsar CTRL+C. ¿Qu´e ocurre si en lugar de incrementarse la variable en 0.1 lo hace en 0.125?. ¿Por qu´e?. 3. Obtener el mayor y el menor n´umero positivo en punto flotante (n´umeros de over flow y underflow ). Para obtener el n´umero de overflow escribir un bucle que vaya calculando las sucesivas potencias de 2 y que finalice cuando se produce overflow. Se recomienda utilizar el comando isinf para detectar cuando se produce el overflow (teclear help isinf ) para obtener informaci´on sobre este comando. Otra instrucci´on que puede resultar u ´ til es break para interrumpir el bucle cuando se produce el over flow . El n´umero de underflow se puede obtener calculando las sucesivas potencias negativas de 2 hasta obtener un n´ umero indistinguible del cero en punto flotante.