Practicas de Analisis Matricial. Una Introduccion a MATLAB

Ion Zaballa

1. Introduccion

Aunque desarrollado con mi propio lenguaje y motivacion, el contenido e ideas de esta Intro-duccion a MATLAB esta sacadas de los Capıtulos 1 de los excelentes libros de Cleve Moler [1]y Charle van Loan [2]. El primero esta disponible por capıtulos separados totalmente gratis enla direccion:

http://www.mathworks.com/moler

Esta introduccion a MATLAB tiene como objetivo que el/la estudiante se familiarice rapidamen-te con algunas de las muchas posibilidades que ofrece MATLAB para trabajar las matematicas.La idea es presentar varios problemas que investigan algunos problemas elementales y al mismotiempo, espero que, interesantes de matematicas. Quienes tengan experiencia en otros lenguajesde programacion enseguida apreciaran la potencia y sencillez de manejo de MATLAB estudian-do estos ejemplos. Para todos, iniciados o no iniciados, se ha confeccionado una Guıa basica deMATLAB para el Curso de Analisis Matricial, disponible en la misma pagina web en la queesta esta introduccion. Pero para quienes esta sea la primera vez que tienen contacto con un len-guaje de programacion, hay multitud de manuales “on-line” y MATLAB, propiamente, disponede una guıa de ayuda para empezar a trabajar:

Ejercicio 1 Pinchando en Help en el menu principal y seleccionando Product Help apa-recera una nueva ventana con toda la ayuda de MATLAB. En la parte izquierda aparece unmenu que se puede desplegar pinchando en el sımbolo . Utilizando esta tecnica tomate unosminutos para familiarizarte con este sistema de ayuda. Por ejemplo, despliega los siguientesmenus: MATLAB → Getting Started → Matrices and Arrays → Matrices and Ma-gic Squares y lee los documentos que aparecen allı: About Matrices, Entering Matrices,sum,transpose and diag, Subscripts, The Colon Operator y The magic Function.Estan ordenados, de forma que se entiende mejor su contenido si se leen en el orden especifica-do.

Ejercicio 2 Tomate otros pocos minutos para ver una de las “demos” que vienen con MATLAB.Se accede a ellas a traves del menu Help que aparece en el menu principal tanto de la ventana decomandos como de la de ayuda. Una demostracion recomendable es Basic Matrix Operationsque aparece al seleccionar Mathematics en el menu principal.

1

Ejercicio 3 .- Utiliza helpwin para echar un vistazo a las funciones matematicas elementalesde MATLAB.

φ

φ − 1

1

1

Figura 1: El rectangulo aureo.

2. La razon aurea

Para ir aprendiendo algunos comandos de MATLAB vamos a tomar como disculpa la razonaurea. Se trata de una practica dirigida en la que debes ir haciendo lo que se indica parafamiliarizarte con algunos comandos basicos de MATLAB.

La razon aurea aparece en muchas partes de la matematica y aquı veremos unas pocas. Larazon aurea toma su nombre del reactangulo aureo que se muestra en la Figura 1. Este tienela propiedad caracterıstica siguiente: al suprimir un cuadrado unidad se obtiene un rectangulomenor con las mismas proporciones que el original. Es decir,

1φ

=φ− 1

1.

La razon aurea es el lado de tal rectangulo. Debe ser entonces la raız cuadrada positiva de lasiguiente ecuacion de segundo grado

φ2 − φ− 1 = 0.

Las dos soluciones de esta ecuacion son

φ =1±√

52

Ejercicio 4 Asigna el valor positivo a la variable phi y calcula el valor en formato short ylong.(Indicacion: Puedes utilizar el comando help format para consultar la ayuda sobre losposibles formatos en que se presentan los numeros).

2

Se puede usar MATLAB para hallar las raıces de un polinomio. MATLAB representa los poli-nomios como un vector: el de sus coeficientes. Ası, el vector

>> p=[1 -1 -1]

representa el polinomiop(x) = x2 − x− 1

Las raıces se calculan mediante la funcion roots:

>> r=roots(p)

produce

r =

-0.618033988749891.61803398874989

Ejercicio 5 Compruebalo.

Hay muchas formas de calcular la razon aurea con MATLAB. Una de ellas es definir la funcion

f(x) =1x− (x− 1)

y calcular sus ceros. Para definir una funcion a partir de la version 7 de MATLAB se puedeutilizar el siguiente codigo

>> f=@(x)1/x-(x-1)

Ası queda definida f como funcion de x. Se pueden definir funciones de varias variables. Porejemplo

>> f=@(x,y)y/x-(x-y)

definirıa f como funcion de las variables x, y. Y

>>f(2,1)

devuelve el valor de f en el punto (2,1).

Ejercicio 6 Escribe los comandos anteriores y comprueba los resultados.

3

Supongamos que tenemos definida la funcion f(x) = 1x − (x − 1). Queremos calcular las raıces

de esta funcion. Para ello necesitamos una estimacion de donde pueden estar aproximadamente.Una forma de obtener esta aproximacion es dibujar la grafica de esta funcion. El comando ezplotlo hace

>> ezplot(f,0,4)

0 0.5 1 1.5 2 2.5 3 3.5 4−3

−2

−1

0

1

2

3

4

5

6

7

x

1/x−(x−1)

Figura 2: Grafica de la fucnion f(x) = 1x − (x− 1).

produciendo la Figura 2 sin el cırculito que marca el cero de la funcion f . Los argumentos 0 y4 especifican el intervalo (0, 4) en el que se dibuja la grafica. Notese que f tiene una asıntotavertical en x = 0.

En la grafica se aprecia que f tiene un cero entre 1 y 2. Poniendo

>> phi=fzero(f,1.5)

se obtiene

phi =

1.61803398874989

que es una aproximacion a la razon aurea tan buena como se podrıa desear. Los comandos

>> hold on>> plot(phi,0,’o’)

sirven para colocar un circulito en en la posicion (φ, 0) en la misma figura donde se habıadibujado la curva (ver la Figura 2).

4

Ejercicio 7 Escribe los comandos anteriores y comprueba los resultados. Usa la ayuda (helphold, help plot) para saber lo que hacen dichas ordenes.

A continuacion escribiremos el programa de MATLAB que produce la Figura 1. Pero antes dehacerlo conviene tener en cuenta lo siguiente: un programa de MATLAB es un fichero con laextension .m que contiene varias lıneas de comandos validos de MATLAB y que se ejecutansucesivamente. Enseguida veremos un ejemplo. La distribucion de MATLAB instalada en elordenador tiene cientos o miles de tales programas. El Capıtulo 4 de la Guıa esta dedicado a laprogramacion en MATLAB y el objetivo de las practicas de este curso es llegar a confeccionaralgunos programas interesantes en el lenguaje de MATLAB. En el Capıtulo 3 se explican algunosconceptos previos necesarios para que MATLAB interactue correctamente con los programas quevayas creando. Esto es un pequeno resumen.

MATLAB necesita conocer de antemano los directorios donde se encuentran los programas quequieres ejecutar. Si los programas vienen en la distribucion, no hay problema, ya sabe dondebuscarlos. Pero no tiene por que saber donde se encuentran los programas confeccionados porcada usuario; ası que este tiene tres opciones:

(a) ponerlos en la carpeta que utiliza MATLAB por defecto. Esto puede cambiar de unaversion a otra y de un sistema operativo a otro. En el momento actual Windows guarda,por defecto, los ficheros de los usuarios en la carpeta MATLAB de Mis Documentos.

(b) colocarlos en una carpeta personal (incluyendo un “lapiz de memoria” o disquete) y anadiresta carpeta personal al MATLAB search path, y

(c) colocarlos en una carpeta personal y cambiar el directorio de trabajo (ventana superiorizquierda) a dicha carpeta personal.

En tu propio ordenador personal cualquiera de las tres opciones es buena, pero en los ordenadoresde la Facultad, en los que muchas personas diferentes pueden usar MATLAB, no tienes garantıasde que un programa dejado en una carpeta del ordenador vaya a encontrarse allı la siguientevez que lo necesites. Es conveniente, por lo tanto, que cada cual venga provisto de un disqueteo “lapiz de memoria” para almacenar allı tus programas.

Por el momento comprueba que estas en la capeta MATLAB de Mis Documentos, y si no, cambiaa ella. Si no sabes como hacerlo preguntalo.

En la ventana de comandos de MATLAB escribe (es muy importante que no te olvides de laextension .m porque es la forma en que MATLAB sabe que se trata de un fichero ejecutable):

>> edit rectaureo.m

Se abrira, en otra ventana, el editor de MATLAB. Escribe en el las siguientes lıneas:

%RECTAUREO Rectangulo aureo% RECTAUREO dibuja el rectangulo aureo.

5

phi = (1+sqrt(5))/2;x = [0 phi phi 0 0];y = [0 0 1 1 0];u = [1 1];v = [0 1];plot(x,y,’b’,u,v,’b--’);text(phi/2,1.05,’\phi’);text((1+phi)/2,-.05,’\phi - 1’);text(-.05,.5,’1’);text(.5,-.05,’1’);axis equal;axis off;set(gcf,’color’,’white’);

A continuacion salvalo. En el directorio de trabajo debe aparecer ahora el fichero rectaureo.m.Este es un ejemplo tıpico de lo que en MATLAB se llama un script (guion, en traduccion literalal castellano). Los script son ficheros con la extension .m que contienen una serie de comandos deMATLAB. Cuando escribimos el nombre del fichero (sin extension) en la ventana de comandosse ejecutan sucesivamente las lıneas del fichero. Observa que despues de cada comando en elfichero rectaureo.m hemos puesto ; porque no queremos que los resultados aparezcan en laventana de comandos; solo queremos que aparezca la figura.

Ahora, en la ventana de comandos de MATLAB escribe:

>> rectaureo

MATLAB ejecutara todos los comandos escritos en el fichero y aparecera la Figura 1 en unanueva ventana.

Ejercicio 8 Trata de adivinar lo que hace cada una de las lıneas del programa y, una vez quelo hayas conseguido (usando la ayuda o preguntando, si es necesario), responde a las siguientescuestiones: En los paıses europeos la medida estandar del papel es DIN A4. Sus dimensiones son:210 mm de ancho y 297 mm de largo. La razon entre estas dos cantidades no es la razon aureapero es proxima a otro numero irracional. ¿Cual?. Si divides por la mitad una hoja DIN A4 ¿cuales la razon entre las dimensiones de cada una de las mitades? Modifica el fichero rectaureo.mpara ilustrar esta propiedad.

Ejercicio 9 Escribe un programa que se llame circunferencia.m que produzca la Figura 3.

La razon aurea φ se puede obtener como una fraccion continua; i.e. una expresion de la forma

a0 +1

a1 + 1a2+ 1

a3+···

.

6

−1 −0.5 0 0.5 1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

u

v

160o

Figura 3: Circunferencia de radio 1

En efecto, resulta que si todos los ai son iguales a 1 obtenemos φ:

φ = 1 +1

1 + 11+ 1

1+···

.

Esto permite confeccionar una funcion de MATLAB que genera y evalua la fraccion continuatruncada al numero de fracciones que se desee. Abre un editor y escribe lo que viene a conti-nuacion. Una vez escrito, guardalo con el nombre fracaurea.m.

function fracaurea(n)%FRACAUREA Fraccion continua para la razon aurea.% FRACAUREA(n) presenta n terminos de la fraccion continua% y varias formas de evaluar la razon aurea.

p = ’1’;for k = 1:n

p = [’1+1/(’ p ’)’];%Se debe dejar un espacio enter ’ y p y entre p y ’endp

p = 1;q = 1;for k = 1:n

s = p;p = p + q;q = s;

endp = sprintf(’%d/%d’,p,q)

format long

7

p = eval(p)

format shorterr = (1+sqrt(5))/2 - p

Uno de los objetivos de este curso es llegar a escribir algunas funciones y scripts en el lenguajede programacion de MATLAB. Es importante que leas atentamente el Capıtulo 4 de la Guıa;en particular, la seccion 4.2. En cualquier caso vamos a analizar un poco esta funcion.

Todas las funciones de MATLAB son ficheros con la extension .m y empiezan con la palabrafunction. Esta es una de las diferencias con los scripts. Otra diferencia es que, suelen recibirvalores (en este caso n) y, aunque no es el caso de esta funcion, tambien suelen devolver valores(enseguida veremos un ejemplo de funcion que cumple estos dos requisitos). Y una terceradiferencia muy importante es que las variables de las funciones son locales mientras que lasde los scripts son globales. Es decir, las variables en las funciones se almacenan en la memoriamientras se ejecuta la funcion; una vez terminada esta desaparecen de la memoria del ordenador.Para los scripts sin embargo, las variables se mantienen en la memoria con el valor que tomanen el script hasta que es sustituıdo por otro valor o se cierra MATLAB. Es decir, las variablesde los scripts actuan como si se ejecutaran en la ventana de comandos.

Analicemos esta funcion un poco. En primer lugar escribimos

>> fracaurea(6)

en la ventana de comandos de MATLAB. La respuesta es

p =

1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))

p =

21/13

p =

1.61538461538462

err =

0.0026

Las tres ps son diferentes formas de representar la misma aproximacion a φ.

8

La primera p es la fraccion continua truncada a 6 terminos. Se produce de la siguiente forma:se le asigna a la variable p el valor ’1’. Las tildes son importantes. Es la forma de decirle aMATLAB que se trata de una variable de tipo string, cadena de caracteres. De esta forma, alprincipio p es el string ’1’. Luego entramos en un bucle for ...end que repetidamente insertalos strings ’1+1/(’ y ’)’ por delante y por detras, respectivamente, del string almacenado en p.Independientemente de lo largo que este string llegue a ser, es una expresion valida en MATLAB.Finalmente, el comando p (sin ;) muestra el valor de p al salir del bucle.

El segundo p es la fraccion ordinaria que produce el truncamiento de la fraccion continua despuesde 6 terminos. Puede no parecerlo, pero si observamos que

1 +1pq

=p+ q

p

enseguida nos daremos cuenta de que empezando con los valores p = q = 1 (i.e., con la fraccion11) y reemplazando la fraccion

p

qpor

p+ q

p

obtenemos el valor de la fraccion continua truncada al termino deseado. Finalmente el comando

p = sprintf(’%d/%d’,p,q)

escribe la fraccion obtenida en la forma p/q con p y q en formato de numero decimal (aunque eneste caso sean enteros). Debe notarse que despues de p = sprintf(’ %d/ %d’,p,q), la variablep vuelve a ser un string.

El tercer p es el valor numerico del segundo p obtenido usando el comando eval. Esta funciondevuelve el valor numerico de un string formado por una operacion valida para MATLAB. Porejemplo, los comandos

>> A=rand(2);>> s=’det(A)’;>> eval(s)

producirıan el valor del det(A) (aunque, claro, esta es una forma muy tonta de calcularlo).

En el caso que estamos considerando, el segundo p es el string ’21/13’. Al hacer eval(p) secalcula el valor de la fraccion 21

13 y se escribe como numero decimal.

Todo lo anterior nos ha servido para aprender algunas funciones de MATLAB pero se podrıahaber obtenido el mismo resultado de forma mas rapida. El codigo

p = sprintf(’%d/%d’,p,q)format longp = eval(p)

9

se podrıa haber sustituıdo por

format ratp=p/qformat longp

De esta forma p serıa siempre una variable numerica, primero en formato racional, y despues,en formato decimal.

Finalmente la cantidad err es la diferencia entre p y φ; es decir, el error absoluto de la estimacion.Con 6 terminos, la aproximacion tiene menos de 3 dıgitos de precision.

Ejercicio 10 ¿Cuantos terminos se necesitarıan para conseguir una precision de 10 dıgitos? ¿Ypara que el error sea cero? ¿Que pasa, entonces, (teoricamente y en la practica) si aumentas elnumero de terminos?

3. Numeros de Fibonacci

La razon aurea aparece en muchos lugares de forma que podrıa considerarse misteriosa. Unejemplo son los numeros de Fibonacci. En un libro publicado en 1202, Liber Abaci, LeonardoPisano Fibonacci propuso el siguiente problema:

Un hombre coloca un par de conejos en un lugar rodeado por un muro de forma quequeden completamente aislados. Este par de conejos y los sucesivos que vayan na-ciendo engendran un nuevo par cada mes a partir del segundo mes de vida. ¿Cuantosconejos habra al cabo de un ano?

Si no se pusiera la condicion de que los conejos no son productivos hasta el segundo mes,la respuesta serıa sencilla y no habrıa habido numeros de Fibonacci. Si cada par de conejosengendrara otro par cada mes, el numero de estos se duplicarıa cada mes; y al cabo de n meseshabrıa 2n conejos. Pero la condicion de que solo se pueden producir nuevos pares de conejos apartir del segundo mes de vida, hace que el problema sea mas interesante.

Si fn representa el numero de pares de conejos despues de n meses, este numero debe ser igualal numero de conejos en el mes anterior fn−1 mas el de conejos que ha producido una nuevapareja (que son los que habıa un mes antes, fn−2). Ası

fn = fn−1 + fn−2

Los dos primeros meses habra f1 = 1 y f2 = 2 pares de conejos, respectivamente.

Con estos datos podemos producir una funcion de MATLAB que nos proporcione la secuenciade Fibonacci con tantos terminos como queramos. Usa el editor para escribirla y salvarla con elnombre de fibonacci.m

10

function f = fibonacci(n)%FIBONACCI secuencia de Fibonacci% f = FIBONACCI(n) genera los primeros n numeros de Fibonacci.f = zeros(n,1);f(1) = 1;f(2) = 2;for k = 3:n

f(k) = f(k-1) + f(k-2);end

Esta ya es una funcion tıpica de MATLAB. Es conveniente y altamente aconsejable que elnombre del fichero que usas para escribir esta funcion (en este caso fibonacci.m) coincida conel nombre de la funcion. Debes tener en cuenta que si llamas a la funcion escribiendo su nombreen la lınea de comandos, MATLAB buscara un fichero con ese nombre. Por lo tanto, si el nombrede la funcion y del fichero no coinciden, no obtendras el resultado que deseas.

fibonacci es una verdadera funcion de MATLAB, sencilla, pero con todos los ingredientes:empieza con la palabra function, admite datos de entrada: n, y produce una salida: f, quees un vector. La funcion define primero un vector (matriz) de tamano n × 1 con todas suscomponentes iguales a 0. A continuacion cambia la primera componente a 1 y la segunda a 2iniciando la secuencia de Fibonacci. Despues entramos en un bucle para sustituir las restantescomponentes del vector f .

Entre el nombre de la funcion y los comandos que la componen hay dos (o mas) lıneas quecomienzan por el sımbolo %. Son lıneas de comentarios que sirven para orientar a quien use lafuncion sobre lo que esta realiza. Si una vez escrita la funcion y salvado el fichero escribimos

>> help fibonacci

en la ventana de comandos, obtendremos

FIBONACCI secuencia de Fibonaccif = FIBONACCI(n) genera los primeros n numeros de Fibonacci.

Esta informacion es muy util para los posibles usuarios de la funcion. Es muy buena costumbrecolocar un par de lıneas de comentarios para ayudarnos a nosotros mismos y a los demas.

Veamos ahora la salida que produce la funcion. Si escribimos

>> g=fibonacci(12)

obtendremos

g =

11

123581321345589144233

La salida de la funcion, que es un vector, se ha almacenado en la variable g, que ahora contienelos 12 primero terminos de la sucesion de Fibonacci.

Por lo general, no necesitaremos todos los datos del vector, quiza solo unos pocos, o ni tansiquiera eso: quiza solo queramos tener el vector con los datos para usarlo posteriormente. Porejemplo, si hacemos

>> g=fibonacci(12);

No hay respuesta visible por parte de MATLAB, pero en g esta el vector con los 12 primerosnumeros de la sucesion de Fibonacci. Aunque no lo parezca estos estan relacionados con φ, larazon aurea. Veamos, escribamos

>> format rat>> g(7)/g(6)

La respuesta es

ans =

21/13

Este numero nos resulta familiar: es el numero que producıa fracaurea(6). ¿Es una casualidad?.

Ejercicio 11 Comprueba para varios valores de n (no mayores que 40, ¿por que?) que la fraccionque se obtiene como resultado de fracaurea(n) coincide con el valor de f(n+1)/f(n) donde fes el vector que devuelve fibonacci(n+1).

Esto hace plausible el siguiente resultado que es, en efecto, verdadero

lımn→∞

fn+1

fn= φ.

12

4. Cuadrados Magicos

Aunque en la actualidad MATLAB puede considerarse un software de proposito general, na-cio con el objetivo de proporcionar un entorno de trabajo para el calculo numerico con vectoresy matrices. De ahı su nombre que proviene de Matrix Laboratory. Los cuadrados magicos pro-porcionan un conjunto interesante de ejemplos de matrices. Si escribimos

>> help magic

obtenemos

MAGIC Magic square.MAGIC(N) is an N-by-N matrix constructed from the integers1 through N^2 with equal row, column, and diagonal sums.Produces valid magic squares for all N > 0 except N = 2.

Los cuadrados magicos eran conocidos en China 2000 anos antes de nuestra era, y el cuadradomagico de orden 3 se conoce con el nombre de Lo Shu. Cuenta la leyenda que fue descubierto enel caparazon de una tortuga en el rıo Lo en el siglo 23 antes de Cristo. Hay cantidad de mitologıay esoterismo asociado a los cuadrados magicos. Una rapida consulta en Google proporciona unsin fin de paginas relacionadas con los cuadrados magicos. Tal y como dice la ayuda de magic estecomando proporciona un cuadrado magico de cualquier orden, salvo de orden 2. Por ejemplo,Lo Shu se consigue de la siguiente forma:

>>A=magic(3)

que produce

A =8 1 63 5 74 9 2

El comando

>> sum(A)

devuelve un vector cuya i-esima componente es la sumas de los elementos de la i-esima columna:

ans =15 15 15

Para obtener la suma de los elementos en cada fila podemos usar el mismo comando sobre latranspuesta:

13

>> sum(A’)’ans =

151515

Para sumar los elementos de la diagonal principal (elementos en la posicion (i, i)):

>> sum(diag(A))ans =

15

El comando flipup pasa la primera fila a la ultima, la segunda a la penultima y ası sucesiva-mente. Si queremos comprobar que la contradiagonal (elementos en la posicion (i, n − i + 1))tambien suman lo mismo podemos utilizar el comando

>> sum(diag(flipud(A)))ans =

15

El mismo efecto habrıamos conseguido con el comando fliplr que hace lo mismo que flipudpero por columnas.

Todo lo anterior certifica que A es un cuadrado magico. Pero ¿por que la suma es 15?. Elcomando

>> sum(1:9)ans =

45

nos permite observar que la suma de los primeros 9 numeros enteros positivos es 45. Como estosestan colocados en A de tres en tres y sumando siempre lo mismo, la suma de cada fila, columna,etc. debe ser 15.

Hay 8 formas posibles de colocar una transparencia en un retroproyector (no olvidemos queuna transparencia se puede leer por delante y por detras). De la misma forma hay 8 manerasposibles de escribir el mismo cuadrado magico de orden 3 que corresponden a las posibles formasde rotar y reflejar el cuadrado. Se pueden obtener mediante los siguientes simples comandos quese explican por sı solos:

>> for k=0:3rot90(A,k)rot90(A’,k)end

14

Ahora un poco de algebra lineal:

>> det(A)ans =-360

nos da el determinante,

>> X=inv(A)X =

0.1472 -0.1444 0.0639-0.0611 0.0222 0.1056-0.0194 0.1889 -0.1028

nos da la inversa, aunque se aprecia mejor si usamos el formato racional

>> format rat>> XX =

53/360 -13/90 23/360-11/180 1/45 19/180-7/360 17/90 -37/360

Volvemos al formato decimal para calcular la norma espectral, los valores propios y los valoressingulares de A:

>> format short>> r=norm(A),e=eig(A),s=svd(A)r =

15e =

15.00004.8990-4.8990

s =15.00006.92823.4641

La suma “magica”, 15, aparece en los tres casos porque el vector (1, 1, 1) es vector propio yvector singular a izquierda y derecha de A.

En el grabado La Melancolıa de Albrecht Durero aparece, entre otros objetos matematicos, uncuadrado magico de orden 4. La distribucion de MATLAB contiene una copia electronica deeste cuadro. Se encuentra en el fichero durer.mat. Carguemos los datos del fichero:

15

>> load durer

Ahora veamos las variables que contiene. El comando whos devuelve las variables que han sidodefinidas en una sesion. Es posible que tengas otras variables, aparte de las definidas al cargarel fichero durer.mat, de modo que lo que aparece a continuacion puede ser solo parte de lo quese muestra en tu ventana de comandos.

>> whosName Size Bytes Class Attributes

X 648x509 2638656 doublecaption 2x28 112 charmap 128x3 3072 double

Hay libros especializados en como MATLAB maneja los colores y las imagenes, pero no estu-diaremos este asunto aquı. Como curiosidad veamos como conseguir el grabado La Melancolıa.Los tres siguientes comandos lo proporcionan:

>> image(X)>> colormap(map)>> axis image

Puedes teclearlo uno a uno para ver lo que vas obteniendo.

Hay una lupa en la barra de herramientas, utilızala para agrandar un cuadrado que se apreciadebajo de la campana. Observaras una matriz 4× 4. Quiza no la distingas muy bien. MATLABtiene otro fichero, detail.mat, en el que estan los datos de ese detalle del grabado con muchamejor resolucion. Los siguientes comandos te proporcionan la imagen:

>> load detail>> image(X);colormap(map);axis image;

Se trata de un cuadrado magico, pero no el que produce MATLAB con el comando magic(4):

>> A=magic(4)A =

16 2 3 135 11 10 89 7 6 124 14 15 1

El de Durero se obtiene de este permutando las columnas segunda y tercera:

>> A=A(:,[1 3 2 4])

16

A =16 3 2 135 10 11 89 6 7 124 15 14 1

Al permutar dos filas o dos columnas se puede destruir la “magia” del cuadrado, pero casualmen-te en este caso no. En efecto, las diagonales siguen sumando 34 (la sumas de filas y columnas, porsupuesto, no se modifican). Probablemente Durero escogio este cuadrado magico porque juntadolas cifras de los dos numeros que estan en el centro de la ultima fila se consigue el numero 1514,que fue el ano en que Durero hizo el grabado.

Hemos visto dos cuadrados magicos diferentes de orden 4; hay 880. Y de orden 5, 275305224.No se sabe cuantos hay de orden 6 o mas.

Curiosamente, o no, para nuestros cuadrados magicos de orden 4 tenemos que det(A) = 0. Y siintentamos calcular la inversa obtenemos

>> inv(A)Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 9.796086e-18.ans =

1.0e+15 *

0.1251 0.3753 -0.3753 -0.1251-0.3753 -1.1259 1.1259 0.37530.3753 1.1259 -1.1259 -0.3753-0.1251 -0.3753 0.3753 0.1251

Ya veremos el significado del mensaje obtenido. En cualquier caso, esto nos indica que haycuadrados magicos que son matrices singulares (no invertibles) ¿cuales?. El siguiente codigoproduce una tabla de dos columnas, en la primera se coloca el orden del cuadrado magicoproducido por MATLAB y en la segunda su rango:

>> for n=1:24; r(n)=rank(magic(n));end; [(1:24)’ r’]ans =

1 12 23 34 35 56 57 78 39 910 7

17

11 1112 313 1314 915 1516 317 1718 1119 1920 321 2122 1323 2324 3

Observa cuidadosamente esta tabla, ignorando el caso n = 2 que no es un cuadrado magi-co.¿Observas alguna regla?. Quiza un diagrama de barras puede ayudarte(Figura 4):

>> bar(r); title(’Rango de los cuadrados magicos’)

0 5 10 15 20 250

5

10

15

20

25Rango de los cuadrados magicos

Figura 4: Rango de los cuadrados magicos producidos por MATLAB

Ejercicio 12 A=magic(4) es singular; i.e., sus columnas son linealmente dependientes. Los si-guientes comandos te informan sobre como es esa dependencia: null(A), null(A,’r’), rref(A).¿Que informacion te proporcionan?

Ejercicio 13 Sea A=magic(n) para n=3,4,5. ¿que efecto producen los comandos

>> p=randperm(n);q=randperm(n);A=A(p,q);

sobre sum(A),sum(A’)’,sum(diag(A)),sum(diag(flipud(A))),rank(A)?

18

5. Procesos Aleatorios

Muchas simulaciones que se realizan para predecir el comportamiento de objetos involucranprocesos aleatorios. El ejemplo que veremos mas abajo es el siguiente: Supongamos que queremosobtener una estimacion del numero π. Podrıamos idear un proceso aleatorio cuya probabilidadeste relacionada con el numero π. Una posibilidad es la siguiente: ¿Cual es la probabilidad deque al lanzar un dardo, sin apuntar, a una diana redonda inscrita en un cuadrado (Figura 5)acertemos en la diana y no en la superficie del cuadrado exterior a la misma? Sencillamente la

−1.5 −1 −0.5 0 0.5 1 1.5−1.5

−1

−0.5

0

0.5

1

1.5

Figura 5: Una diana.

relacion entre las areas de ambas figuras geometricas. Si el cuadrado lo pensamos centrado enel punto (0, 0) y de lado 2, la circunferencia tendra radio 1 y tal probabilidad es

π

4.

Con un ordenador podemos simular el lanzamiento de un dardo a la diana: cada lanzamientoes asociado a la eleccion aleatoria de un par de numeros (x, y) comprendidos entre −1 y 1.Este punto representarıa el punto del cuadrado donde se ha clavado el dardo. Realizando estasimulacion un numero muy alto de veces la frecuencia relativa de eleccion de puntos en el interiordel cırculo se debera aproximar a la probabilidad buscada:

π ≈ Numero de lanzamientos en el interior del cırculoNumero total de lanzamientos

· 4

El punto clave en esta simulacion es el de la eleccion aleatoria del par de puntos que deter-minan las coordenadas del cuadrado. Todos los puntos deberıan tener la misma probabilidadde ser elegidos; y esto nos indica que los numeros aleatorios que los determinan deben estaruniformemente distribuidos. MATLAB proporciona, mediante el uso repetido del comando randuna sucesion de puntos uniformemente distribuıdos entre 0 y 1. Como es generada tal sucesionno es un problema trivial, pero se puede modificar el comienzo de la misma con el comandorand(’seed’,n) con n un numero entero positivo.

19

Tambien se pueden obtener numeros aleatorios normalmente distribuıdos. El comando que lohace es randn. El siguiente script puede servir para comprender la diferencia entre ambos tiposde numeros aleatorios:

close allclc

subplot(2,1,1)x=rand(1000,1);hist(x,30)axis([-1 2 0 60])title(’Distribucion de los valores en rand(1000,1)’)xlabel(sprintf(’Media= %5.3f. Mediana=%5.3f.’, mean(x),median(x)))

subplot(2,1,2)x=randn(1000,1);hist(x,linspace(-2.9,2.9,100))title(’Distribucion de los valores en randn(1000,1)’)xlabel(sprintf(’Media= %5.3f. Desviacion standard=%5.3f.’, mean(x),std(x)))

Ejercicio 14 Copia este script en un fichero con el nombre randhisto.m y ejecutalo. Trata decomprender lo que hace cada comando.

Ejercicio 15 .- Para apreciar, de otra forma, la diferencia entre la distribucion uniforme y lanormal escribe un “script” que haga lo siguiente:1.- Escoja aleatoriamente 10000 pares de puntos, primero con distribucion uniforme y despuescon distribucion normal.2.- Haga dos graficas, cada una de ellas dibujando todos los puntos escogidos de cada una de lasdos formas.3.- Ponga un tıtulo significativo en cada grafica.4.- Los ejes de la primera grafica deben estar entre 0 y 1; y los de la segunda entre -3 y 3.

Una vez que entendemos mejor la diferencia entre los comandos rand y randn, podemos volvera nuestra simulacion para el calculo aproximado del numero π. Vamos a escribir una funcionque produzca la simulacion del lanzamiento de los dardos. Recordemos que la parte importantees la seleccion aleatoria de un par de puntos (x, y) en el interior del cuadrado (por lo tanto, condistribucion uniforme) y la determinacion de si tal punto esta o no en el interior del cırculo. Lafuncion admite un numero entero positivo

function pi=dardos(n)

%DARDOS es una funcion para calcular el numero pi% usando un metodo de Montecarlo

20

%DARDOS admite un numero entero n que es el numero de%cientos de lanzamientos y devuelve una%aproximacion a pi y la grafica de la convergencia a%pi del metodo

if (ceil(n)-n ~= 0 || n<0)error(’n debe ser entero’);return

elseclose all; clcrand(’seed’,0);ndentro=0;piestimada=zeros(n,1);for k=1:n

x=-1+2*rand(100,1);y=-1+2*rand(100,1);ndentro=ndentro+sum(x.^2+y.^2<=1);piestimada(k)=(ndentro)/(100*k)*4;

endplot(piestimada)pi=piestimada(n);title(sprintf(’Valor estimado de pi= %5.6f’,pi));xlabel(’Cientos de lanzamientos’);

end

La ejecucion de esta funcion con n = 500 (esto es, 50000 lanzamientos) produce el valor deπ = 3,1302 y la Figura 6. Debe observarse que el valor estimado de π mejora gradualmente a

0 50 100 150 200 250 300 350 400 450 5003.08

3.1

3.12

3.14

3.16

3.18

3.2

3.22

3.24

3.26Valor estimado de pi= 3.157

Cientos de lanzamientos

Figura 6: Una estimacion de π de tipo Montecarlo.

medida que n aumenta, pero que el “progreso” hacia 3’14159 no es ni estable ni, mucho menos,rapido. Simulaciones de este tipo se conocen con el nombre de Montecarlo (o Monte Carlo).

21

Hay una sentencia en el programa que merece atencion especial: sum(x.^2+y.^2<=1). x y y sonvectores de 100 componentes cada uno y x.^2 eleva cada componente de x al cuadrado (y lomismo y.^2). Ahora, x.^2+y.^2 es un vector cuyas componentes son la suma de los cuadrados delas componentes de x y de y. Finalmente, x.^2+y.^2<=1 tambien es un vector: sus componentesson 0 o 1. Es 0 si la correspondiente componente no cumple la condicion y 1 si la cumple. Porlo tanto, sum((x.^2+y.^2<=1) cuenta el numero de componentes del vector x.^2+y.^2 que sonmenores que 1; es decir, el numero de dardos que han impactado en el interior del cırculo.

Otros dos comandos utiles para trabajar con vectores de ceros y unos son any y all. any(x)devuelve un 1 (“true”) si alguna componente de x es distinta de cero; si no devuelve un 0(“false”). Y all(x) devuelve un 1 si todas las componentes de x son distintas de cero; si nodevuelve un cero. Si x es una matriz, any y all actua por columnas. Ası

>> x=[ 0 1 0; 1 -1 0; -1 1 0]x =

0 1 01 -1 0-1 1 0

>> any(x)ans =

1 1 0>> all(x)ans =

0 1 0

Ejercicio 16 ¿Cual es la probabilidad (aproximadamente) de que el polinomio cuadratico p(x) =ax2 + bx+ c tenga raıces complejas sabiendo que los coeficientes a, b y c son variables aleatoriascon distribucion uniforme (0, 1)? ¿Y si lo son con distribucion normal (0, 1)?. Se trata de fabricaruna funcion que lo simule.

6. Aritmetica en punto flotante

Dado que los ordenadores utilizan un numero finito de bits para representar los numeros reales,solo pueden representar un numero finito de ellos; es decir, un subconjunto finito de numerosreales (o de numeros complejos). Esta limitacion conlleva dos dificultades:

Los numeros representados no pueden ser arbitrariamente grandes (en valor absoluto).

Debe haber “agujeros” entre ellos.

Lo primero conduce a los fenomenos conocidos con los nombres de “overflow” y “underflow”(rebosamiento por arriba o por abajo) y lo segundo produce el “roundoff” (redondeo). Por logeneral, es posible usar MATLAB sin preocuparse de estas cuestiones, pero de vez en cuando

22

puede que nos encontremos con problemas que un somero conocimiento del funcionamiento dela aritmetica en punto flotante puede ayudar a solucionar.

Hace 20 anos la situacion era mucho mas complicada de lo que es ahora porque cada ordenadortenıa su propio sistema de numeros en punto flotante. Algunos eran binarios, otros decimales.Incluso hubo algunos ordenadores rusos que usaban aritmetica ternaria. Ademas, entre los bi-narios los habıa que usaban aritmetica de base 2, otros octal o hexadecimal. Y cada uno tenıasu propia precision. En 1985 la IEEE Standard Board y el American National Standard Insti-tute decidieron adoptar para la aritmetica binaria de punto flotante lo que se conocce como elANSI/IEEE Standard 754-1985. Esta decision fue la culminacion del trabajo de casi una deca-da de trabajo de un grupo de 92 personas integrado por matematicos, ingenieros e informaticosprocedentes de universidades y empresas dedicadas a la construccion de ordenadores y micropro-cesadores. Todos los ordenadores disenados a partir de 1985 usan la aritmetica IEEE de puntoflotante. Esto no significa que todos ellos producen los mismos resultados porque hay una ciertaflexibilidad en la norma ANSI/IEEE, sino que el modelo de funcionamiento de la aritmetica enpunto flotante es, hoy en dia, independiente del ordenador que se utilice.

MATLAB ha usado siempre el formato IEEE de doble precision. Hay un formato de precisionsimple que salva espacio pero que no supone una mejora en la velocidad con los ordenadoresactuales. Nosotros hablaremos solo de la precision doble. Para entender bien este conceptodebemos recordar la forma en la que los ordenadores trabajan con los numeros en punto flotante.La mayorıa estan normalizados. Si calculamos 1/

√120 en una calculadora de bolsillo con notacion

cientıfica obtendremos (esto es lo que marca la mıa)

9, 128709292e− 2 (1)

o algo similar. Es la notacion cientıfica para el numero

9, 12870929 · 10−2.

El numero en punto flotante de (1) no esta normalizado. Para considerarlo normalizado debemosescribirlo de la siguiente forma:

0, 912870929 · 10−1 o 0, 912870929e− 1.

Es decir, un numero en punto flotante es normalizado si la parte fraccionaria es de la forma0.x1x2 . . . con x1 6= 0. Para los numeros representados en base 10, 1 ≤ x1 ≤ 9, pero para losnumeros en base 2 x1 = 1 siempre, de modo que la parte fraccionaria tendra la forma 0,1x2x3 . . .con xi = 0 o 1. Notese que un numero de esta forma puede escribirse como (1 + f)2−1 siendof = 0.x2x3 . . .. Por lo tanto, en un ordenador que representa los numeros en base binaria (que esel que consideraremos nosotros), los numeros en punto flotante normalizados tienen la siguienteforma:

x = ±(1 + f) · 2e

con0 ≤ f < 1 y e un numero entero.

Al numero f se le llama fraccion o mantisa y a e exponente o caracterıstica.

Los numeros en punto flotante de precision simple se almacenan en “palabras” de 32 bits,mientras que los de doble precision lo hacen en “palabras” de 64 bits. De estos, 52 bits se usan

23

para almacenar f , 11 bits para e y 1 bit para el signo del numero, 0 para + y 1 para − (Figura7). Puesto que la fraccion debe estar entre 0 y 1 y se debe almacenar en 52 bits, debe resultar

s fracciónexponente

0 1 11 63

Figura 7: Palabra en punto flotante en doble precision.

que 252f es un numero entero que debe estar en el intervalo

0 ≤ 252f < 252

Y como 252 ≈ 4,5 · 1015, los numeros en punto flotante en doble precision se expresan, en basedecimal, con 15 decimales. Debe notarse, finalmente que la parte fraccional completa de losnumeros en punto flotante no es f sino 1+f , que necesita 53 bits. Sin embargo, como la primeracifra decimal siempre es 1 (si el numero esta normalizado), no es necesario almacenarlo. Dehecho, el formato IEEE empaqueta 65 bits de informacion en una “palabra” de 64 bits.

Por otra parte, para el exponente se reservan 11 bits, ası que e es un numero entero en el intervalo

0 ≤ e ≤ 211 = 2048.

Sin embargo, debemos permitir exponentes negativos almacenandolos en los 11 bits disponibles.Para conseguirlo lo que se hace es acomodar en los 11 bits el numero e+1023. Ası, si el numeroalmacenado en los 11 bits de e es 1250 entonces e = 1250 − 1023 = 227; mientras que si elnumero almacenado es 1000 entonces e = 1000− 1023 = −23. Los valores extremos posibles enestos 11 bits: 0 y 211 − 1 se reservan para definir numeros especiales que veremos mas adelante.Los restantes son valores admisibles para e. Por lo tanto el exponente debe encontrarse en elintervalo

−1022 ≤ e ≤ 1023.

Teniendo en cuenta que

(1 + 0,99999999999999) · 21023 ≈ 1,79 · 10308 y (1 + 0) · 2−1022 ≈ 2,23 · 10−308,

estos son el mayor y menor numeros representables en doble precision. En MATLAB estosnumeros reciben un nombre especial

realmax = (1 + 0,99999999999999) · 21023 = 1,7977 · 10308

realmin = 1 + 0) · 2−1022 = 2,2251 · 10−308.

Si algun calculo intenta producir un valor mayor que realmax se dice que produce un rebosa-miento por arriba (“overflow”) dando como resultado un numero especial llamado “infinity”,Inf en MATLAB. Se representa tomando e = 1024 y f = 0 y cumple propiedades tales como1/Inf=0 o Inf+Inf=Inf. Por ejemplo

24

>> realmax^2ans =

Inf>> Inf*Infans =

Inf

Para operaciones tales como 0/0 o Inf-Inf el resultado es un valor excepcional denotado porel sımbolo NaN por “not a number”.

>> 0/0Warning: Divide by zero.ans =

NaN>> Inf-Infans =

NaN

Si algun calculo intenta producir un numero inferior a realmin, se dice que produce un re-bosamiento por abajo (“underflow”). En realidad, algunos sistemas como MATLAB permitennumeros mas pequenos llamados numeros en punto flotante denormales (o subnormales, aun-que esta denominacion en castellano tiene connotaciones negativas). El numero mas pequenoen MATLAB es 2 ^ (-1074) ≈ 0.494e-323. Veremos mas adelante de donde proviene estenumero.

>> 2^(-1074)ans =4.9407e-324

> 2^(-1075)ans =

0

El rebosamiento por arriba (que rara vez ocurre en la practica) puede, a veces, evitarse escalandoadecuadamente los datos del problema para convertirlo en un rebosamiento por abajo que notiene consecuencias importantes. Por ejemplo, supongamos que tenemos que calcular la cantidad

c =√a2 + b2

con a = 10200 y b = 1. Una respuesta aceptable serıa c = 10200 debido a que b es insignificanteen comparacion a a. Si hacemos este calculo directamente en MATLAB:

>> a=10^200;b=1; c=sqrt(a^2+b^2)c =

Inf

25

La respuesta es debida a que se ha producido un rebosamiento por arriba. Teoricamente, elmismo resultado se obtendrıa si hicieramos

c = s

√(as

)2+(b

s

)2

con s 6= 0. En particular, escogiendo s = max{|a|, |b|} = 10200 obtenemos lo siguiente:

> s=max(abs(a),abs(b));c=s*sqrt((a/s)^2+(b/s)^2)c =1.0000e+200

que es la respuesta que hemos dado como aceptable. Este resultado es debido al rebosamientopor abajo ocurrido al calcular (b/s)2:

>> (b/s)^2ans =

0

Los fenomenos de rebosamiento tienen su origen en las cotas superior e inferior de los exponentesadmisibles. Hay otro fenomeno importante que es consecuencia de la finitud de las posiblesmantisas o fracciones: el redondeo. Comencemos con la siguiente pregunta: ¿Cuantos numerosen coma flotante hay entre 1 y 2? Son los de la forma

+(1 + f) · 20.

Y como f esta representado por 52 bits (i.e., 52 posiciones cada una de ellas conteniendo un 0o un 1), el numero posible de fracciones es 252 = 4,503599627370496e+ 15 (las formas posiblesde colocar 0 y 1 en 52 posiciones distintas). Ası pues, los posibles valores de f son:

0, 2−52, 2 · 2−52, 3 · 2−52, 4 · 2−52, . . . , (252 − 1) · 2−52.

En consecuencia, los numeros en punto flotante en el intervalo [1, 2] son

1, 1 + 2−52, 1 + 2 · 2−52, 1 + 3 · 2−52, 1 + 4 · 2−52, . . . , 1 + (252 − 1) · 2−52, 2. (2)

De la misma forma, los numeros en punto flotante en el intervalo [2, 4) son de la forma

+(1 + f) · 21,

y son, concretamente:

2, 2 + 2−51, 2 + 2 · 2−51, 2 + 3 · 2−51, 2 + 4 · 2−51, . . . , 2 + (252 − 1) · 2−51, 4.

En general, los numeros en el intervalo [2j , 2j+1] son los de (2) multiplicados por 2j .

26

1/16 1/2 1 2 4 8−1/2

||||||||||||||||||||||||| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | || |

eps = 1/8

Figura 8: Floatgui.

Ejercicio 17 En el directorio matlab de tu cuenta personal encontraras un archivo tituladofloatgui. Ha sido escrito por Cleve Moler para poner de manifiesto la distribucion de losnumeros positivos en punto flotante en un sistema para el que es posible elegir los parametrost (numero de bits utilizado para almacenar f) y emax y emın (el intervalo en el que esta e).Inicialmente floatgui se presenta con t = 3, emın = −4 y emax = 3 y produce la distribucionde la Figura 8 Se puede observar que en cada intervalo binario [2j , 2j+1] los numeros estanigualmente espaciados con un incremento de 2j−t. Si j = 0 y t = 3, por ejemplo, el espacio entrelos numeros que estan entre 1 y 2 es 2−3; el espacio entre los que estan entre 2 y 3 es 2−2, etc.

Esta “equidistribucion” se hace mas evidente si utilizamos una escala logarıtmica. La Figura 9muestra el resultado de floatgui cuando se marca la casilla de log scale con t = 5, emın = −4y emax = 3.

1/16 1/4 1/2 1 2 4 16−1/4

||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||| |||||||||||||||||||||||||||||||| |

eps = 1/32

Figura 9: Floatgui.

Experimenta con floatgui para comprender la distribucion de los numeros en punto flotante yresponde a las siguientes preguntas: ¿Para un valores concretos de t, emın y emax cuantos numerospositivos son representables exactamente en punto flotante? ¿Cuantos numeros positivos contieneel conjunto de numeros en punto flotante del sistema IEEE?

Cualquier numero real que no coincida con uno de los numeros en punto flotante del sistemaIEEE debe ser aproximado por uno de ellos; concretamente por el mas proximo. A este fenomenose le llama redondeo. Si x es un numero real, se denota por fl(x), el numero en punto flotantemas proximo a x. Resulta que si realmin≤ x ≤realmax entonces

| fl(x)− x||x|

≤ 2−53.

En efecto, suponiendo que x ∈ [2j , 2j+1] y que fl(x) = (1 + a · 2−52)2j entonces

(1 + a · 2−52)2j ≤ x < (1 + (a+ 1) · 2−52)2j o (1 + (a− 1) · 2−52)2j < x ≤ (1 + a · 2−52)2j

27

Pero como fl(x) = (1 + a · 2−52)2j se tiene

x− (1 + a · 2−52)2j < (1 + (a+ 1) · 2−52)2j − x = (1 + a · 2−52)2j − x+ 2j−52

o(1 + a · 2−52)2j − x < x− (1 + (a− 1) · 2−52)2j = x− (1 + a · 2−52)2j − x+ 2j−52.

En cualquier caso

2|x− (1 + a · 2−52)2j | < 2j−52 ⇔ |x− (1 + a · 2−52)2j | < 2j−53

Como fl(x) = (1 + a · 2−52)2j y x > 2j concluımos que

| fl(x)− x||x|

≤ 2−53.

Un razonamiento similar nos permite ver que, en terminos relativos, la distancia entre dosnumeros en punto flotante consecutivos es 2−52. En particular, esta es la distancia de 1 alsiguiente numero en punto flotante tal y como hemos visto en (2). A este numero, se le conocecon el nombre de epsilon de la maquina y se denota en MATLAB por eps. Ası

eps = 2−52 ≈ 2,22 · 10−16

y se obtiene en MATLAB con el comando eps:

>> epsans =

2.220446049250313e-16.

Algunos autores prefieren llamar epsilon de la maquina al maximo error relativo posible alredondear un numero real por el numero en punto flotante mas proximo; i. e. eps/2=2−53. Encualquier caso, el nivel de error relativo es de aproximadamente 16 dıgitos decimales.

Notese que 2−52 · 2−1022 = 2−1074. Este numero ya lo hemos visto antes: era el menor valordenormal (o subnormal) admisible en MATLAB.

En el calculo binario del numero t=0.1 se produce un redondeo que suele afectar a muchasoperaciones aritmeticas. El valor almacenado en t no es exactamente 0,1 porque la fracciondecimal 1/10 expresada en base 2 requiere una serie infinita. Esto se entiende mejor si expresamos1/10 en base 16; es decir, si usamos la representacion hexadecimal :

110

= 2−4

(1 +

916

+9

162+

9163

+ · · ·). (3)

En efecto∞∑

j=1

16j

=1/6

1− 1/6=

115,

y

2−4

(1 +

915

)=

2416 · 15

=110.

28

Para pasar (3) a formato binario observamos que

2−4 916

=23 + 1

28=

125

+128,

2−4 9162

=23 + 1

212=

129

+1

212,

etc.

(4)

Por lo tanto110

=124

+125

+026

+027

+128

+129

+0

210+

0211

+1

212+ · · · .

La representacion en numero de punto flotante de 1/10 se obtiene tomando 52 terminos de estaserie (en formato binario), o 13 terminos de la serie en formato hexadecimal, y redondeando elultimo termino al entero mas proximo. En formato hexadecimal tendrıamos que si

t1 = 2−4

(1 +

916

+9

162+

9163

+ · · ·+ 91613

)

t2 = 2−4

(1 +

916

+9

162+

9163

+ · · ·+ 101613

)entonces

t1 < 1/10 < t2.

Pero resulta que 1/10 esta mas proximo a t2 que a t1. Es decir,

1/10 = (1 + f) · 2e

dondef =

916

+9

162+

9163

+ · · ·+ 101613

e = −4

(5)

Hemos escrito f en formato hexadecimal. Para pasarlo a binario y conseguir los 52 bits que lorepresentan se procederıa como en (4). En realidad, la conversion de numeros entre las bases 2 y16 es muy sencilla. Por ejemplo, en consonancia con (5), el comando de MATLAB format hexaplicado a t=0.1 produce:

>> t=0.1t =

0.1000>> format hex>> tt =

3fb999999999999a

29

Las letras a a f representan los numeros 10 a 15 en base hexadecimal. Para pasar t a binariobasta observar que cada cifra es una “palabra” de cuatro bits. Ası 3fb en binario serıa

0011 1111 1011,

que son los 12 primeros bits de t. El primer bit, 0, es el signo; por lo tanto t es positivo. Losprimeros 12 bits para −0,1 serıa

1011 1111 1011,

que en hexadecimal es bfb. Los restantes 11 bits representan el exponente (en realidad, recor-demos, 1023+e). Como 3fb en base diez es 3 · 162 + 15 · 16 + 11 = 1019, resulta que e = −4. Lasrestantes cifras 999999999999a representan la fraccion que en binario serıa

1001 1001 1001 · · · 1001 1010

ocupando los 52 bits correspondientes.

Para saber la representacion hexadecimal de un numero en MATLAB se puede usar el comandoformat hex. Una vez que se usa todas los resultados se presentaran en formato hexadecimal.Para volver al formato decimal se deben usar los comandos format short, o format long, ocualquiera de los otros disponibles. Ası, el siguiente codigo nos da la razon aurea en base 16:

>> format hex>> phi=(1+sqrt(5))/2phi =

3ff9e3779b97f4a8

Una vez mas, 3 es 0011 en binario, de modo que φ es positivo. Ademas 3ff en base 10 es 1023,por lo que e = 0. Las restantes 13 cifras hexadecimales contienen la fraccion f :

f =916

+14162

+ · · ·+ 101612

+8

1613.

Con estos valores de e y f(1 + f) · 2e ≈ φ.

Ejercicio 18 1. ¿Cual es la representacion hexadecimal de los numeros 3/10 y 3 · 0,1? ¿Porque no es la misma?. (Si eres capaz de deducir por que MATLAB produce esos resultados,mejor que mejor). ¿Cual es el error absoluto, tanto en formato hexadecimal como decimal,cometido al aproximar 3/10 por 3 · 0,1?.

2. Explica el resultado producido por el siguiente codigo de MATLAB

t = 0.1n = 1:10e = n/10 - n*t

30

Dos ultimos ejemplos en los que los efectos del redondeo se ponen de manifiesto. El primeroconsiste en la resolucion del siguiente sencillo sistema

17x1 + 5x2 = 221,7x1 + 0,5x2 = 2,2

(6)

Dado que la primera ecuacion es 10 veces la segunda, este sistema tiene infinitas soluciones y lamas simple es x1 = x2 = 1.

Ejercicio 19 1. Resuelve el sistema (6) con MATLAB y comprueba que la solucion obtenidano es x1 = x2 = 1.. (Indicacion: MATLAB utiliza el comando A\b para resolver el sistemaAx = b mediante eliminacion gaussiana.)

2. Utiliza MATLAB para hacer la eliminacion gaussiana en etapas con formato cientıfico; i.e.format short e:

>> A(2,:)=A(2,:)-1.7/17*A(1,:)>> b(2)=b(2)-1.7/17*b(1)

Observa que los elementos a22 y b2 no son cero. ¿Por que?

3. Como a22 y b2 no son cero, x2 = b2/a22. Calcula con MATLAB este valor y despues calculatambien x1.

El ultimo ejemplo de los efectos del redondeo lo vemos en el siguiente ejercicio

Ejercicio 20 El comando conv(p,q) multiplica los polinomios p y q. Por ejemplo, si p(x) =2x3 − x+ 3 y q(x) = −2x2 − x+ 2 entonces

>>p=[2 0 -1 3]; q=[-2 -1 2];>>conv(p,q)ans =

-4 -2 6 -5 -5 6

de modo que el polinomio p(x)q(x) = −4x5 − 2x4 + 6x3 + · · · .

Considera el siguiente polinomio p(x) = (x − 2)9 y utiliza el comando conv para calcular loscoeficientes de este polinomio, Llamalo q.

1. Dibuja la grafica de q(x) para x = 1,920, 1,921, 1,922, . . . , 2,080.

2. En la misma figura, dibuja en rojo p(x) en los mismos puntos. ¿Que conclusion sacas?.

Referencias

[1] C. B. Moler: Numerical computing with MATLAB. SIAM, 2004.

[2] Ch. F. van Loan: Introduction to scientific computing. Prentice Hall.

31