MatlabOctave Aprender MatlabOctave en 25 Horas

MATLAB + OCTAVE

APRENDER MATLAB/OCTAVE EN 25 HORAS

TERESA ARIAS-MARCO JOSÉ LUIS BRAVO IGNACIO OJEDA Mª ISABEL PARRA CONCHITA MARÍN © Editorial Tébar. Prohibida la reproducción sin la autorización expresa de la editorial

  

© Editorial Tébar. Prohibida la reproducción sin la autorización expresa de la editorial

             ISBN digital: 978-84-7360-472-7   

                                   ISBN digital: 978-84-7360-472-7      © Editorial Tébar. Prohibida la reproducción sin la autorización expresa de la editorial

Índice 1. Introducci´ on a MATLAB/Octave 1.1. Primeros pasos . . . . . . . . . . . . . . . . . 1.2. Atajos de teclado. El comando help . . . . . 1.3. Escalares, operaciones aritméticas y variables 1.4. Operadores y funciones . . . . . . . . . . . . 1.5. Gestión de ficheros . . . . . . . . . . . . . . . 1.6. Ejercicios . . . . . . . . . . . . . . . . . . . . 1.7. Soluciones . . . . . . . . . . . . . . . . . . . . 2. Matrices y vectores 2.1. Definición desde el teclado . . . . . . . . . . . 2.2. Otras formas de definir matrices y vectores . 2.3. Acceso a los elementos de una matriz . . . . . 2.4. Operaciones con vectores y matrices . . . . . 2.5. Cálculos con polinomios . . . . . . . . . . . . 2.6. Resoluci´ on de sistemas de ecuaciones lineales 2.7. Gestión de ficheros de hojas de c´ alculo . . . . 2.8. Ejercicios . . . . . . . . . . . . . . . . . . . . 2.9. Soluciones . . . . . . . . . . . . . . . . . . . . 3. Gr´ aficos 3.1. Creación y manipulaci´ on b´ asica de gráficos 3.2. Creaci´ on especializada de gr´ aficos . . . . . . 3.2.1. Gr´ aficos 2D . . . . . . . . . . . . . . 3.2.2. Gr´ aficos 3D . . . . . . . . . . . . . . 3.2.3. Gr´ aficos Estad´ısticos . . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

5


. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

. . . . . . . . .

. . . . .

. . . . . . .

9 9 12 13 17 20 24 26

. . . . . . . . .

29 29 31 34 35 40 42 45 50 53

. . . . .

59 59 72 73 82 85

6

Aprende MATLAB/OCTAVE en 25 horas 3.3. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4. Soluciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. Introducci´ on a la programaci´ on 4.1. “Scripts” . . . . . . . . . . . . 4.2. Funciones . . . . . . . . . . . . 4.3. Estructuras de control . . . . . 4.4. Entrada y salida . . . . . . . . 4.5. Manejo de funciones . . . . . . 4.6. Ejercicios . . . . . . . . . . . . 4.7. Soluciones . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

91 93 97 97 100 104 109 112 114 116

5. Ajuste de curvas a datos 119 5.1. La funci´ on polyfit . . . . . . . . . . . . . . . . . . . . . . . . 119 6. Autovalores y autovectores 125 6.1. La funci´ on eig . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.2. Un ejemplo “cl´ asico” . . . . . . . . . . . . . . . . . . . . . . . . 127 7. Animaciones con MATLAB 131 7.1. Creaci´ on de animaciones . . . . . . . . . . . . . . . . . . . . . . 131 7.2. Pel´ıculas o MOVIES . . . . . . . . . . . . . . . . . . . . . . . . 132 7.3. Ejemplos de pel´ıculas o MOVIES . . . . . . . . . . . . . . . . . 133 8. Ecuaciones diferenciales 137 8.1. Ecuaciones Diferenciales Ordinarias . . . . . . . . . . . . . . . . 137 8.2. Ecuaciones en derivadas parciales . . . . . . . . . . . . . . . . . 139 Bibliograf´ıa


143

Pr´ ologo Este libro surge a partir del material elaborado para impartir un curso intensivo de MATLAB/Octave para profesores de la Universidad de Extremadura. La primera edici´ on de este curso tuvo lugar en mayo de 2010 con una muy buena acogida por parte de los profesores participantes, lo que propició una segunda edición en 2011. Los contenidos expuestos se presentan como un curso de auto-aprendizaje de modo que el lector puede conseguir una visión global de las posibilidades de este software al tiempo que puede ir ejercitándose en el uso del mismo. As´ı, es fundamental que se realice una lectura atenta delante de un ordenador con MATLAB u Octave instalado. El libro consta de dos partes bien diferenciadas. La primera parte, cap´ıtulos del 1 al 4, corresponde a los aspectos generales e imprescindibles para el manejo de MATLAB/Octave. Cada uno de los cap´ıtulos de esta parte contiene una serie de ejercicios de auto-evaluación a fin de que el lector pueda poner en práctica los contenidos tratados. La segunda parte del libro muestra algunos ejemplos concretos del uso de este software en el ámbito académico universitario. Esperamos que esta obra contribuya a la difusión del uso del software matemático (representado aqu´ı por MATLAB/Octave) como instrumento eficaz de aprendizaje en ´ ambitos docentes y como herramienta fundamental para abordar sin temor las matem´ aticas.

Badajoz, 16 de julio de 2011 7


´ A MATLAB/OCTAVE CAPÍTULO 1. INTRODUCCION

9

Cap´ıtulo 1

Introducci´ on a MATLAB/Octave ´ Este y todos los dem´ as cap´ıtulos del manual están pensados para ser trabajados delante de un ordenador con MATLAB u Octave instalado1 , y no para ser le´ıdos como una novela. En vez de eso, cada vez que se presente un comando de MATLAB/Octave, se debe introducir el comando, pulsar la tecla para ejecutarlo y ver el resultado. Más a´ un, se desea que se verifique el resultado. Aseg´ urese de que se comprende perfectamente lo que se obtiene antes de continuar con la lectura.

1.1.

Primeros pasos

Los primeros p´ arrafos de esta secci´ on, as´ı como algunas frases aisladas en las siguientes, est´ an tomados, con ligeras variaciones, del libro de A. Quarteroni y F. Saleri “C´ alculo Cient´ıfico con MATLAB y Octave”, citado en la bibliograf´ıa; una fuente excelente para introducirse en el conocimiento de los métodos numéricos m´ as usados en el c´ alculo cient´ıfico. MATLAB y Octave, los programas, son entornos integrados para el cálculo cient´ıfico y la visualizaci´ on de datos. Ambos están escritos en los lenguajes C y C++. 1

Concretamente estas pr´ acticas est´ an realizadas para MATLAB 7.0 (R14) y Octave 3.2.3


10

Aprende MATLAB/OCTAVE en 25 horas

MATLAB est´ a distribuido por The MathWorks (http://www.mathworks. com). El nombre proviene de MATrix LABoratory, pues originariamente fue desarrollado para el c´ alculo matricial. Octave, también conocido como GNU Octave (http://www.gnu.org/software/octave/), es un software que se distribuye libremente sujeto a los términos de la licencia GPL (GNU Public License). A lo largo del manual haremos uso frecuente de la expresión “comando de MATLAB”; en ese caso, MATLAB deber´ a ser entendido como el lenguaje que es el subconjunto com´ un a ambos programas MATLAB y Octave. Una vez llevada a cabo la instalaci´ on, la ejecución de MATLAB y Octave nos da acceso a un entorno de trabajo caracterizado por los indicadores (prompt) >> y octave:1>, respectivamente. Después de presionar la tecla (enter, return), todo lo que esté escrito a continuaci´ on del indicador ser´ a interpretado2 . Si la instrucción es incorrecta, MATLAB devuelve un aviso de error. En caso contrario, el comando es ejecutado y si procede se mostrar´ a una salida. Ejemplo 1. Escribe 2 + 2 en el indicador de MATLAB, pulsa la tecla intro y observa el resultado: 2+2 Escribe ahora pepe y pulsa intro  pepe ¿Qué ha ocurrido? Finalmente, el sistema devuelve el indicador para poner de manifiesto que está preparado para un nuevo comando. En todo caso, conviene saber que cualquier comando puede abortarse pulsando simultáneamente las teclas control () y (). Para cerrar una sesión de MATLAB se debe escribir el comando quit (o exit) y pulsar la tecla intro . En adelante se entenderá que para ejecutar un programa o un comando se tiene que pulsar 2

Por tanto, un programa en MATLAB no necesita ser compilado como ocurre con otros lenguajes (Fortran o C, por ejemplo).



11

. Adem´ as, los términos programa, función o comando se utilizarán de forma equivalente. Más adelante veremos que la sintaxis de MATLAB no es muy distinta de la de cualquier otro lenguaje; las diferencias serán las habituales, los s´ımbolos de continuación, los comentarios, ... Veamos algunos de estos s´ımbolos: Las comillas simples o ap´ ostrofes ’-’ sirven para introducir un texto literal. Ejemplo 2. Escribe ’pepe’  ’pepe’ Obsérvese que MATLAB no nos ha devuelto ahora un mensaje de error, sino que ha devuelto como salida la variable por defecto ans (abreviatura de answer) al igual que ocurri´ o en el ejemplo 1 cuando escribimos 2 + 2. Si ahora escribimos una cadena de caracteres (n´ umero o comando) diferente, la variable ans tomará este nuevo valor. Ejemplo 3.  ’Bienvenido a la UEx’ El porcentaje % es el s´ımbolo usado para los comentarios. Todo lo que está por detrás de este s´ımbolo en una l´ınea es ignorado por el interprete. Ejemplo 4. 2 + 2

% Esto es un comentario

Si a˜ nadimos tres puntos ... al final de una l´ınea significa que se tomará la siguiente l´ınea como continuaci´ on. Ejemplo 5.  2 + 2 ... - 3


12


Se pueden escribir varios comandos en una misma l´ınea, separándolos por una coma , o por un punto y coma ; El punto y coma sirve además para inhibir que el resultado se muestre en la pantalla. Ejemplo 6.  2 + 2, 2 - 1, ans + 3, 3 + 5; Nótese que, aunque hemos realizado cuatro sumas, sólo se han mostrado los resultados de las tres primeras ya que hemos prohibido que se muestre la salida de la u ´ltima. No obstante, MATLAB ha interpretado la orden y ha guardado resultado en la variable por defecto. Ejemplo 7.  ans

1.2.

Atajos de teclado. El comando help

La consola de MATLAB es un entorno de trabajo más potente de lo que parece gracias a una serie de atajos de teclado de gran utilidad. Ejemplo 1. En la l´ınea de comandos de MATLAB, pulsa la flecha hacia arriba una vez ¿qué ha ocurrido? Pulsa ahora la flecha hacia arriba varias veces y observa lo que sucede, repite el experimento pulsando también la flecha hacia abajo . Puedes recuperar todas tus órdenes anteriores y volver a ejecutarlas si lo deseas. También se pueden usar las flechas izquierda y derecha del teclado para mover el cursor a derecha e izquierda y editar as´ı una l´ınea de comandos. De este modo se pueden corregir errores, recuperando un comando con la flecha hacia arriba y usando a continuación la edición para corregirlo. Como es natural se pueden cambiar tantos s´ımbolos como se quiera o utilizar la tecla de retroceso para borrar.



13

Ejemplo 2.  sqrt(5)

% Ra´ ız cuadrada de 5

Utilizando la flecha arriba podemos editar el comando sqrt(5), cambiar el 5 por un 4 y calcular as´ı el valor de la ra´ız cuadrada de 4.  sqrt(4) MATLAB posee una amplia colecci´ on de mensajes de ayuda. Si conoces el nombre de un comando, se puede encontrar la sintaxis correcta para el comando y ejemplos de su uso escribiendo help nombre comando. Ejemplo 3.  help sqrt  help help Una forma alternativa de obtener ayuda consiste en utilizar la tecla de tabulación : Ejemplo 4. Escribe sq en la l´ınea de comandos y pulsa para obtener la lista funciones y comandos que comienzan por “sq”. Obviamente sqrt es uno de ellos. Finalmente, el comando lookfor nos muestra una lista de comandos en cuyas explicaciones aparezca cadena de texto especificada. Ejemplo 5.  lookfor ’square root’  help lookfor

1.3.

Escalares, operaciones aritm´ eticas y variables

A decir verdad MATLAB no distingue entre escalares y matrices. Para mayor claridad, nosotros s´ı haremos tal distinción.


14

Aprende MATLAB/OCTAVE en 25 horas Esencialmente podemos utilizar n´ umeros enteros, reales y complejos: Ejemplo 1.    

15, -21 % Son n´ umeros enteros 0.32, -1.27, .065 % Son n´ umeros reales -5.32e+5, 4.78e-3 % Tambi´ en son n´ umeros reales 1-i, 2+3i, -3+j % Son n´ umeros complejos

A pesar de que en apariencia MATLAB usa cuatro cifras decimales para representar los n´ umeros reales y complejos, la representación interna del n´ umero se hace con 16 cifras decimales. De modo que lo que hemos visto es simplemente uno de los varios posibles formatos de salida de MATLAB (el utilizado por defecto). El mismo n´ umero puede tener diferentes expresiones decimales dependiendo de la declaraci´ on espec´ıfica de formato que se haga. Ejemplo 2. Veamos distintas representaciones de 1/7        

format 1/7 format 1/7 format 1/7 format 1/7

long short e long e rat

El formato por defecto es el formato corto (4 cifras decimales). El comando format a secas nos devuelve al formato por defecto.  format El n´ umero π est´ a predefinido en MATLAB,  pi Dos representaciones aritméticas predefinidas que pueden ser de utilidad conocer son Inf y NaN, la primera representa un n´ umero más grande que cual-



15

quier n´ umero que MATLAB pueda almacenar3 y juega el papel de infinito, mientras que la segunda representa lo que solemos denominar una indeterminación (Not a Number).

Ejemplo 3.  1/0  0/0  Inf + 1 Las convenciones para las operaciones aritméticas en MATLAB son similares a las de cualquier otro lenguaje de programación. El orden de asociación de las operaciones es también el mismo. Primero operan las funciones matemáticas elementales (senos, cosenos, logaritmos, ...), las multiplicaciones y divisiones y luego las sumas y restas. Como es natural, los paréntesis sirven para modificar el orden normal de las operaciones a realizar.

Ejemplo 4. Calculemos 2 0,11/2

1 −

0,4 21/3

− 3(1 − 2i)

 1/((2/0.1^(1/2))-(0.4/2^(1/3)))-3*(1-2i) El comando = permite la asignaci´ on de un valor (o de una cadena de caracteres) a una variable dada.

3

Los n´ umeros reales positivos m´ as grande xmax y m´ as peque˜ no xmin en MATLAB son realmax y realmin, respectivamente. El n´ umero m´ as pr´ oximo a uno dado x es x + eps(x). Se puede comprobar que los n´ umeros reales en MATLAB son m´ as densos cerca de xmin y menos densos cuando se aproximan a xmax .


16


Ejemplo 5. Para asignar la cadena ’Bienvenido a la UEx’ a la variable saludo podemos escribir:  saludo = ’Bienvenido a la UEx’ Para asignar el valor de resultado de la suma 2 + 3 a la variable a podemos escribir: a = 2 + 3 De este modo no hay necesidad de declarar el tipo de una variable, MATLAB lo hará automática y din´ amicamente. Por ejemplo, si escribimos  saludo = 3 la variable saludo contendr´ a ahora un n´ umero y no una cadena de caracteres. Esta flexibilidad no es gratis. Si ponemos una variable de nombre quit igual al n´ umero 5 estamos inhibiendo el comando de MATLAB quit. Por consiguiente, deber´ıamos tratar de evitar el uso de variables que tengan el nombre de comandos de MATLAB. Sin embargo, mediante el comando clear seguido del nombre de una variable, es posible cancelar cualquier asignación y restaurar, por ejemplo, el significado original del comando quit. Ejemplo 6.    

quit = 5 quit clear quit quit

La función clear a secas cancela todas las asignaciones de valores a variables que hayamos realizado. Los nombres de las variables pueden tener a lo sumo 63 caracteres alfanuméricos (no se admiten ni tildes ni la letra n ˜ y s´ı el guión bajo ). Es importante tener presente que el primer carácter siempre debe ser una letra y que para MATLAB las letras may´ usculas son distintas de las min´ usculas.



17

Ejemplo 7.       

a = 5 hola = 3 Hola hola a 1 = 12 2a = -3 a2 = -3

El comando clc limpia la consola de comandos manteniendo los valores asignados a variables. Ejemplo 8.  clc a Como se habr´ a observado para conocer el valor de una variable previamente definida basta con teclear su nombre. También se pueden conocer todas las variables definidas hasta el momento usando los comandos who y whos. Ejemplo 9.    

1.4.

who whos clear who

% Cancelamos los valores dados a las variables

Operadores de comparaci´ on y l´ ogicos. Algunas funciones elementales

Como cualquier lenguaje de programaci´ on, MATLAB dispone de operadores que nos permiten comparar valores.


18


Ejemplo 1. 2 < 3 2 > 3

Obsérvese que la respuesta ha sido 1 para la afirmación que es verdadera y 0 para la falsa. Otros operadores de comparación son <= (menor o igual que), >= (mayor o igual que), == (igual a) y ~= (distinto de). En todos los casos la salida es 1 si la afirmaci´ on es cierta y 0 si es falsa.

Ejemplo 2.    

2 >= 3 4 >= 3 Inf == Inf+1 3 ~= 4

Los operadores l´ ogicos de MATLAB son: & o and (conjunción), | o or (disyunción), ~ (negaci´ on l´ ogica), xor (disyunción exclusiva). Cuando se utilizan estos operadores, se eval´ uan siempre ambos operandos. Sin embargo, MATLAB ofrece la posibilidad de utilizar los operadores breves: &&, para la conjunción y || para la disyunci´ on, que simplifica la operación de comparación evitando operaciones innecesarias, y errores que se producir´ıan caso de evaluar incondicionalmente el segundo argumento (por ejemplo, para evitar la división por cero, podemos utilizar (b~=0)&&(a/b>0)). Sólo eval´ uan el segundo operando en caso de ser imprescindible, cuando el resultado del primer operando sea 1 para la conjunci´ on o 0 para la disyunción.



19

Ejemplo 3.       

2 <= 3 && 2 >= 3 && 2 >= 3 || ~(2 <= 3) ~(2 >= 3) xor(2 <3, xor(2 <3,

4 <5 4 <5 4 <5

2 >3) 2 >1)

MATLAB posee una importante biblioteca de funciones elementales cuyos argumentos pueden ser n´ umeros enteros, reales o complejos, siempre que tenga sentido. El resultado se devuelve en el mismo tipo del argumento. Algunas funciones elementales soportadas por MATLAB son4 : sqrt abs exp log log10 sign

ra´ız cuadrada valor absoluto o m´ odulo exponencial logaritmo natural logaritmo decimal signo

sin cos tan asin acos atan

seno cos tangente arcoseno arcocoseno arcotangente

Algunas funciones relacionadas con los n´ umeros complejos son conj, real, imag, angle y compass. Ejemplo 4.       4

z = 2+3i conj(z) real(z) imag(z) angle(z) compass(z)

Todas las funciones trigonométricas de la tabla operan en radianes.


20


Ahora podemos utilizar las operaciones y las funciones elementales introducidas arriba para construir nuevas funciones. Para ello definimos lo que se conoce como una funci´ on an´ onima. Las funciones anónimas responden a la siguiente sintaxis: nombre funcion = @(argumentos) expresi´ on Veamos algunos ejemplos: Ejemplo 5.     

fun1 = @(x) x^2-1+exp(x) fun1(0) fun2 = @(x,y) cos(x)*sin(y) fun2(pi,pi/3) fun2(pi/3,pi/2)

Un ejemplo de una funci´ on a trozos podr´ıa ser Ejemplo 6.  fun3 = @(x) exp(x)*(x<=0) + x*(x>0)

1.5.

Gesti´ on de ficheros

Todo trabajo que se realiza durante una sesión MATLAB, se puede escribir en un archivo a modo de diario. Este archivo contendrá una transcripción completa y literal de todos los comandos que ha interpretado MATLAB y sus correspondientes respuestas. El contenido del archivo se puede editar con cualquier procesador de textos. Para iniciar un archivo de diario llamado, por ejemplo, prueba.txt , debemos escribir



21

 diary prueba.txt en la l´ınea de comandos de MATLAB y, como siempre, pulsar . Probablemente no se percibir´ a ning´ un cambio, pero el ordenador ha comenzado a registrar todo lo que ocurre en la consola de MATLAB en el fichero prueba.txt  %Esto se incluir´ a en nuestro diario En cualquier momento se puede suspender esta orden escribiendo  diary off Escribiendo ahora diary on (o diary prueba.txt) se volverá a guardar toda la sesión en el archivo de diario que previamente se hab´ıa creado, a˜ nadiendo la informaci´ on nueva a continuaci´ on de la ya registrada. Si hay un diario activo al salir de MATLAB, el archivo de diario será cerrado automáticamente. Al salir de una sesi´ on de MATLAB, los valores de todas las variables se pierden. El archivo de diario s´ olo mantiene un registro de lo que ha aparecido en la pantalla, no se almacenan los valores de las variables. Si se desea, antes de salir de MATLAB, se pueden guardar los valores actuales de las variables en un archivo para su uso en una sesi´ on futura. El siguiente comando guarda los valores de todas las variables existentes en un archivo que llamaremos mis datos.mat; por defecto, MATLAB (pero no Octave) a˜ nade la extensión .mat si no ha sido especificada.  save mis datos.mat Los valores de todas las variables se guardarán en forma binaria (ilegible, por tanto, para los seres humanos) en un archivo llamado mis datos.mat dentro del directorio de trabajo. Para cargar estos valores de las variables en una futura sesión de MATLAB se escribe


22


 load mis datos Es posible salvar el valor de una serie de variables determinadas sin necesidad de guardar el resto.

Ejemplo 1.         

a = 1; b = sqrt(2); c = ’no se guarda’; save var ab.mat a b clear a b c load var ab a, b c

Finalmente, podemos pedir a MATLAB que nos indique cuál es el directorio de trabajo actual con la orden pwd y su contenido con la orden dir. El contenido de un fichero concreto se muestra con el comando type seguido del nombre del fichero. También podemos borrar ficheros del directorio de trabajo usando el comando delete y cambiar de directorio con el comando cd de forma similar a como se hac´ıa en MS-DOS.

Ejemplo 2.    

dir trab=pwd dir type var ab.mat delete var ab.mat



 dir cd(’C:\’)  cd(’home’)  dir  pwd  cd(dir trab)  pwd  dir

23

% Si est´ as usando Windows y hay unidad C. % Si est´ as usando Linux.

En MATLAB existen funciones especificas para importar datos de hojas de c´ alculo generadas con otros programas. El abanico de formatos que reconoce MATLAB es mucho m´ as amplio que en Octave. Esencialmente, Octave reconoce bases de datos separados por delimitadores, por ejemplo comas (ficheros .csv). Mientras que MATLAB puede leer los datos de hojas de cálculo más sofisticadas, por ejemplo de Excel (ficheros .xls). En la siguiente sesi´ on mostraremos brevemente el uso de las funciones csvread y dlmread (v´ alidas en MATLAB y Octave) y de xlsread (exclusiva de MATLAB). De forma an´ aloga, los datos numéricos generados en MATLAB se pueden exportar a hojas de c´ alculo para su lectura y manipulación con otros programas. En este caso, la funciones que veremos serán csvwrite y dlmwrite (para MATLAB y Octave) y xlswrite (s´ olo para MATLAB).


24


1.6.

Ejercicios

Antes de comenzar escribe lo siguiente5  clear all  diary capitulo1.txt Ejercicio 1. Asigna a una variable a el n´ umero uno. Sobre esa variable calcula a = 2a. Repite esta operaci´ on hasta alcanzar el valor de 128. Ejercicio 2. 1. Escribe un n´ umero entero menor que 777. 2. Div´ıdelo entre 7 (si la divisi´ on es exacta elige otro n´ umero y vuelve a dividirlo). 3. Suma los valores de los seis primeros decimales. Si no nos hemos equivocado te ha salido 27. Ejercicio 3. Calcula el valor de ((1 + x) − 1)/x para x = 10−15 (obviamente el resultado es 1 para todo x = 0). Ejercicio 4. Define tres variables a, b y c con valores 10308 , 1.1·10308 y −1.001· 10308 , respectivamente. Compara, usando los operadores de comparación, el valor de a + (b + c) con el de (a + b) + c (evidentemente, por la propiedad asociativa de la suma, son iguales para todo a, b y c). Ejercicio 5. Define una variable z cuyo valor sea un n´ umero complejo no real. √ Calcula su módulo usando la conocida f´ ormula |z| = z · z¯ (la barra superior es la conjugación). El m´ odulo de un n´ umero complejo se puede calcular de forma directa usando la funci´ on abs. Compara el valor obtenido antes con el valor que te da la funci´ on abs. 5

Si no deseas hacer un diario de la sesi´ on, no hace falta que lo escribas. En otro caso, el fichero que se genere se guaradar´ a en el directorio en el que te encuentres; recuerda que puedes saber cu´ al es escribiendo pwd .



25

Ejercicio 6. Usando la ayuda de MATLAB (comando help), explica la diferencia entre los comandos floor, ceil, round y fix. Escribe tus explicaciones con comentarios, %, en la ventada de comandos e il´ ustrarla con ejemplos en MATLAB. Ejercicio 7. Calcula el coseno de 43 grados [Pista: usa el comando help de MATLAB para encontrar la funci´ on apropiada]  diary off


26

1.7.


Soluciones

Ejercicio 1. a = 1 a = 1  a = 2*a a = 2  a = 2*a a = 4  a = 2*a a = 8  a = 2*a a = 16  a = 2*a a = 32  a = 2*a a = 64  a = 2*a a = 128

Ejercicio 2.  format long  776/7 ans=110.857142857143  8+5+7+1+4+2 ans=27  format short Ejercicio 3.  x = 10^(-15);  ((1+x)-1)/x ans = 1.1102 Observamos que el resultado se desv´ıa en un 11 por ciento del valor esperado. El problema es que estamos dividiendo por un n´ umero muy próximo a cero.



27

Ejercicio 4.  a = 10^308; b = 1.1*10^308; c = -1.001*10^308;  a + (b + c) == (a + b) + c ¿Falso? Deberia ser cierto. Veamos que ha ocurrido:  a + (b + c) ans = 1.0990e+308  (a + b) + c ans = Inf Al sumar a y b hemos sobrepasado el infinito de MATLAB. Ejercicio 5.  z = 3+i;  mod z = sqrt(z*conj(z)) mod z = 3.1623  abs(z) ans = 3.1623  abs(z) == mod z ans = 1 Ejercicio 6.  help floor floor(x) el mayor n´ umero entero menor o igual que x  help ceil ceil(x) el menor n´ umero entero mayor o igual que x  help round round(x) es el n´ umero entero m´ as pr´ oximo a x  help fix fix(x) n´ umero (entero) que est´ a a la izquierda de la coma en la expresión decimal de x Ejercicio 7.  help cos cos(x) calcula el coseno de x en radianes  help cosd cosd(x) calcula el coseno de x en grados


28


 cosd(43) ans = 0.73135


CAPÍTULO 2. MATRICES Y VECTORES

29

Cap´ıtulo 2

Matrices y vectores Como ya mencionamos en el cap´ıtulo anterior, MATLAB u Octave son fundamentalmente programas para el c´ alculo matricial. En realidad, todos los datos son matrices, siendo los vectores y escalares casos particulares. Comenzaremos esta segunda sesi´ on aprendiendo a introducir matrices y vectores desde el teclado, para continuar, m´ as adelante, con otras formas más potentes de definir matrices. As´ımismo, experimentaremos con algunas de la funciones incorporadas en MATLAB, que lo convierten en una potente calculadora matricial.

2.1.

Definici´ on desde el teclado

Para definir una matriz no hace falta declararla o establecer de antemano su tama˜ no. De hecho, se puede definir un tama˜ no y cambiarlo posteriormente. Las matrices se definen o introducen por filas1 ; los elementos de una fila se separan por espacios en blanco o por comas. Para cambiar de fila se puede escribir punto y coma o . 1

Sin embargo, se almacenan por columnas. Este dato ser´ a de utilidad m´ as adelante, para comprender c´ omo act´ uan ciertas funciones.


30


Ejemplo 1. Introduce el siguiente comando, que define una matriz 4x3:  A=[1,2,3;4,5,6;7,8,9;10,11,12] O bien  A = [ 1 2 3 4 5 6 7 8 9 ] Una sencilla operaci´ on con matrices es calcular su traspuesta. Para ello, basta con escribir su nombre seguido del apóstrofo (’). Ejemplo 2.  A’ Como el resultado de la operaci´ on no ha sido asignado a ninguna variable, MATLAB utiliza, por defecto, la variable ans, que puede ser usada posteriormente. Sin embargo, una opci´ on m´ as adecuada generalmente hubiera sido asignarle un nombre: Ejemplo 3.  B=A’ MATLAB determina el n´ umero de filas y de columnas en función del n´ umero de elementos que se proporcionan (o utilizan). Ejemplo 4.  C(4,2)=-1 En el ejemplo, hemos definido una matriz de orden 4 × 2 con todos los elementos nulos excepto el segundo elemento de la cuarta fila, cuyo valor es -1.



31

De forma an´ aloga a las matrices, es posible definir un vector fila x (con los n´ umeros separados por comas o espacios en blanco) o un vector columna y (con los n´ umeros separados por puntos y coma). Para pasar de un vector fila a uno columna, basta con usar la funci´ on traspuesta. Ejemplo 5.  x=[1 2 3]  x’  y=[1;4;7]  y’

2.2.

Otras formas de definir matrices y vectores

Introducir matrices a través del teclado sólo es práctico en algunas situaciones, por ejemplo cuando su tama˜ no es peque˜ no o no es necesario repetirlo muchas veces. MATLAB dispone de varias funciones orientadas a definir matrices y vectores de tipos particulares2 , algunas de las cuales son las siguientes: eye zeros ones rand diag linspace :

matriz unidad matriz de ceros matriz de unos matriz de n´ umeros aleatorios entre 0 y 1 matriz diagonal o diagonal de una matriz vector de partici´ on de un intervalo vector de sucesi´ on aritmética

2

Para obtener una lista de todas las matrices elementales escribir help elmat. Este comando no funciona en Octave. Para obtener informaci´ on detallada sobre una de dichas funciones en concreto, tanto en MATLAB como en Octave, escribir help seguido del nombre de la funci´ on.


32


Ejemplo 1. Construir algunas matrices, utilizando funciones:  A=eye(5)  B=ones(3)  D=rand(3,6)  E=zeros(3,4)  F=diag([-1,0,1]) Para construir vectores, podemos tener en cuenta que se trata de un caso particular de matrices, o bien utilizar alguna de las funciones espec´ıficas para ello. El operador : es muy importante en MATLAB, y puede usarse de varias formas, que iremos practicando a lo largo de la sesión. Para empezar, nos permite generar vectores cuyos elementos son los valores contenidos entre dos valores, inferior y superior, a incrementos regulares, que por defecto serán la unidad. Ejemplo 2.  x=1:100

% incremento por defecto

 x1=1:2:100

% incremento de dos unidades

 x2=2:2:100 Puede utilizarse con valores enteros y reales, positivos o negativos. Ejemplo 3.  y=100:-10:1  z=0:pi/4:2*pi Una opción alternativa para generar vectores con valores igualmente espaciados es utilizar la funci´ on linspace.



33

Ejemplo 4. Veamos c´ omo generar los vectores de los ejemplos anteriores, utilizando linspace:  x=linspace(1,100)

% por defecto 99 particiones

 x1=linspace(1,99,50)

% 50 particiones

 x2=linspace(2,100,50)  y=linspace(100,10,10)  z=linspace(10,2*pi,9) Ejemplo 5. Compara los vectores resultantes de las siguientes instrucciones:  x=1:3:10  y=linspace(1,10,4) y de estas otras  x=1:4:10  y=linspace(1,10,3) Finalmente, MATLAB también ofrece la posibilidad de definir una matriz a partir de otra u otras, definidas previamente, recibiendo alguna de sus propiedades (por ejemplo el tama˜ no, su diagonal, ...), por composición de varias submatrices o modificando de alguna forma otra matriz. Ejemplo 6.  A=rand(4)  B=diag(diag(A))  C=[A B]  D=[A;B]  E=[A,eye(4);zeros(4),B]


34


Una caracter´ıstica peculiar de algunas funciones en MATLAB es que dependiendo del tipo de argumentos que reciban, as´ı serán los valores de retorno, determinándose en tiempo de ejecuci´ on. Un ejemplo claro es la función diag, que hemos utilizado en el ejemplo anterior para construir la matriz B. Al aplicarla sobre una matriz, devuelve un vector con su diagonal principal, mientras que al aplicarla sobre un vector, nos devuelve la matriz diagonal. Por defecto, dicha diagonal ser´ a la principal, aunque admite un segundo argumento para referirse a las diagonales secundarias.

2.3.

Acceso a los elementos de una matriz

En MATLAB se accede a los elementos de un vector o matriz escribiendo el ´ındice o los ´ındices correspondientes entre paréntesis. Como ya mencionamos, las matrices se almacenan por columnas, aunque se introduzcan por filas, por lo que se puede acceder a cualquiera de sus elementos con un u ńico ´ındice. Ejemplo 1.  A=rand(5)  A(3,2)  A(7)  A(8)  x=[-2,-1,0,1,2]  x(3) Observar que se ha obtenido el mismo valor escribiendo A(3,2) que A(8). También es posible acceder a un bloque de componentes de una matriz o vector, indicando los ´ındices m´ınimo y m´ aximo o indicando un subconjunto de ´ındices, En este contexto, el s´ımbolo : indica todo el rango.



35

Ejemplo 2.  x(2:4)  A(1:3,[1,3,5])  A(:,1)

% Primera columna de la matriz A

 A(3,:)

% Tercera fila de la matriz A

Acceder a elementos de un vector o matriz de esta forma tiene abundantes aplicaciones, por ejemplo la de modificar selectivamente esos elementos. Ejemplo 3.  A=rand(6)  A(:,2)=1

% Transforma en 1 los elementos de la segunda columna

 A(1:2:6,:)=0

% Transforma en 0 los elementos de las filas impares

Como ya hemos mencionado, podemos redimensionar las matrices en cualquier momento, a˜ nadiendo o eliminando componentes. Para esa u ´ltima acción utilizaremos la matriz vac´ıa (representada en MATLAB por [ ]). Ejemplo 4.  B=[]

% Define una matriz vac´ ıa llamada B

 A(:,7)=7

% A~ nade una columna de sietes

 A(2:2:6,:)=[]

2.4.

% Elimina las filas pares

Operaciones con vectores y matrices

MATLAB puede operar con matrices por medio de operadores y de funciones.


36


Los operadores aritm´ eticos son: + (adición o suma), - (sustracción o resta), * (multiplicaci´ on), ^ (potenciaci´ on). Todos estos operadores son coherentes con las correspondientes operaciones matriciales. Ejemplo 1.  A=[1 2;3 4]  A^2  A*A  A/2 Los cuatro u ´ltimos operadores, precedidos de un punto (.), representan la misma operación pero llevada a cabo elemento a elemento. Ejemplo 2. Comparar con los resultados anteriores:  A.^2  A.*A Como ya vimos en la pr´ actica anterior, también se aplican a variables o valores escalares, y, por supuesto, de modo mixto operando un escalar con una matriz o vector. Ejemplo 3.  A=[0 -2 3 5;7 2 -3 4;1 5 -6 -7]  A+2  A*2 Si se usan de modo incorrecto (por ejemplo, intentando operar matrices con dimensiones incompatibles) se obtiene un mensaje de error.



37

Ejemplo 4.  B=[-1,0,1,2;3,2,1,0;2,4,6,8]  x=[0;1;-1;0]  B*x  x*B

Al igual que con los n´ umeros, los operadores relacionales son: < (menor), <= (menor o igual), > (mayor), >= (mayor o igual), == (igual), ~= (distinto) y los operadores l´ ogicos son: & o and (conjunción), | o or (disyunción), ~ (negación lógica), xor (disyunci´ on exclusiva). Cuando se aplican a dos vectores o matrices del mismo tama˜ no, la operación se realiza elemento a elemento y el resultado es otra matriz, del mismo tama˜ no, de ceros y unos, que recoge el resultado de cada operaci´ on.

Ejemplo 5.  A=[1,2;3,4]  B=[-1,2;3,-4]  C=(A==B)  D=(A~=B)  C+D

% Debe ser una matriz de unos

Los operadores l´ ogicos se combinan con los relacionales para poder comprobar el cumplimiento de condiciones m´ ultiples.


38


Ejemplo 6. Analiza el resultado obtenido al ejecutar las siguientes instrucciones. (Téngase en cuenta que los operadores lógicos breves && y || no funcionan con matrices)  A=rand(4);A(:,5)=1;A(:,6)=1/2;A(:,7)=0  (A>1/4)&&(A<3/4)

% Dar´ a error

 (A>1/4)&(A<3/4)  (1/4
% Observar que el resultado no es el esperado

Finalmente, trabajaremos con algunas de las m´ ultiples funciones intr´ınsecas de MATLAB, especialmente aquellas que son espec´ıficas de vectores y matrices. Las funciones matem´ aticas elementales, algunas de las cuales mostramos en el cap´ıtulo anterior, al aplicarlas sobre una matriz, act´ uan sobre cada elemento de la matriz como si se tratase de un escalar. Las siguientes funciones s´ olo act´ uan sobre vectores (no sobre matrices3 , ni sobre escalares). sum prod dot cross mean max min length

suma los elementos de un vector multiplica los elementos de un vector producto escalar de vectores producto vectorial de vectores calcula la media aritmética valor m´ aximo de las componentes de un vector valor m´ınimo de las componentes de un vector n´ umero de componentes de un vector

3

Cuando una de estas funciones (excepto dot y cross) se aplica sobre una matriz, en realidad, la respuesta es un vector fila cuyas componentes son el resultado de aplicarla por separado a cada columna de la matriz.



Ejemplo 7.  A=magic(5)

39

% Matriz definida en MATLAB

 sum(A)  x=[-5,3,0,2,1,-2]  sum(x)  mean(x)  mean(x’)  exp(x) Cuando una funci´ on tiene en su definici´ on más de un argumento de salida, puede ser utilizada de varias formas. Las funciones max y min nos proporcionan un ejemplo. Si utilizamos la funci´ on en la forma max(x) MATLAB devuelve el máximo valor de las componentes del vector, pero si utilizamos la forma [m,i]=max(x), el valor m´ aximo quedar´ a almacenado en la variable m y la posición que ocupa en la variable i.

Ejemplo 8.  max(x)  [m,i]=max(x) Algunas funciones espec´ıficas para el manejo de matrices son:

size inv rank diag det

da el n´ umero de filas y columnas calcula la inversa de una matriz rango obtiene una diagonal calcula el determinante de la matriz


40


Ejemplo 9.  D=eye(5)  B=diag([1,2,3,4],1)  C=diag([5,6,7,8],-1)  A=D+B+C  IA=inv(A)  A*IA  det(A)  rank(A)  size(A)

2.5.

C´ alculos con polinomios

Para MATLAB, un polinomio se puede definir mediante un vector con sus coeficientes, permitiendo adem´ as realizar diversas operaciones sobre él, como por ejemplo: roots poly

polyval conv deconv polyder

Calcula las ra´ıces del polinomio que recibe como argumento. En general, calcula aproximaciones. Si recibe como argumento un vector columna, devuelve un polinomio que tiene dichas ra´ıces. Si el argumento es una matriz cuadrada, la salida será el polinomio caracter´ıstico. Eval´ ua el polinomio para un determinado valor. Calcula el producto de dos polinomios. Calcula la divisi´ on de dos polinomios. Devuelve tanto el cociente como el resto. Devuelve la derivada del polinomio.



41

Ejemplo 1. Por ejemplo, el polinomio 3x4 − 8x2 + 6x − 1 se puede representar por el vector p  p=[3,0,-8,6,-1]  r=roots(p)  polyval(p,[0,r(1),r(2),5])  poly(r)  3*poly(r)

% Volvemos al polinomio original

Tomemos ahora como ejemplo los polinomios: p(x) = 4x3 − x y q(x) = x2 + 6 y operemos con ellos. Para sumar y restar polinomios, hay que buscar una estrategia para que los dos vectores que los representan tengan el mismo tama˜ no, puesto que en otro caso obtenemos un mensaje de error. Ejemplo 2.  p=[4,0,-1,0]  q=[1,0,6]  q1=[0,1,0,6]  p+q  p+q1  conv(p,q)

% Multiplicaci´ on de polinomios

 [c,r]=deconv(p,q)

% Cociente de polinomios

 polyder (p)


42


2.6.

Resoluci´ on de sistemas de ecuaciones lineales

MATLAB utiliza el operador divisi´ on a izquierda (\) (división-derecha) para calcular una soluci´ on particular de cualquier sistema de ecuaciones lineales, ¡incluso si el sistema es incompatible! En este caso, calculará una solución aproximada; concretamente, el vector de norma m´ınima que esté a menor distancia del espacio de soluciones del sistema. Un sistema de ecuaciones lineales, con m ecuaciones y n incógnitas,

a11 x1 a21 x1 .. .

+ +

a12 x2 a22 x2 .. .

+ ··· + + ··· + .. .

a1n xn a2n xn .. .

= =

b1 b2 .. .

am1 x1 + am2 x2 + · · · + amn xn = bm

          

se puede escribir en forma matricial

Ax = b

donde x y b son vectores columna y A una matriz. Una formulación tan sencilla, puede esconder complicaciones. Aparentemente, la solución es: x = A−1 b

que en MATLAB ser´ıa x=inv(A)*b. Sin embargo no siempre es posible calcular la inversa de la matriz A. Una opción alternativa es utilizar la división izquierda (x=A\b) que calcula una solución del sistema, aunque la matriz no tenga inversa, o incluso no sea cuadrada, tras realizar un an´ alisis previo de la matriz A para decidir cuál es el método más adecuado.



43

Ejemplo 1. Resolver el siguiente sistema de ecuaciones:  4x + 3y + 2z = 5   2x − 7y + 2z = −7  3x + 2y + 8z = −3 

 A=[4 3 2; 2 -7 2; 3 2 8]  b=[5; -7; -3]  sol1=inv(A)*b  sol2=A\b

Para comprobar que efectivamente se trata de la solución, podemos calcular:  A*sol1-b  A*sol2-b

Por otra parte, en el caso de un sistema de ecuaciones incompatible, la solución particular es el punto m´ as cercano, de norma m´ınima, aunque no cumpla exactamente ninguna de las ecuaciones.

Ejemplo 2. Analizar qué ocurre al resolver el sistema de ecuaciones:  x + 2y = 2   x = 0  y = 0   A=[1 2; 1 0; 0 1]  b=[2; 0; 0]  sol=A\b


44


Como el sistema es incompatible la solución obtenida es una solución aproximada. Para ver el error cometido:  norm(A*sol-b) En otras ocasiones, MATLAB nos avisa que la solución puede ser incorrecta.

Ejemplo 3.  A=[1 2 3;4 5 6;7 8 9], b=[1;0;0]  A\b

Una vez sabemos c´ omo calcular una solución particular de un sistema de ecuaciones lineales podemos hallar todas las soluciones del sistema usando la función null. La funci´ on null devuelve una base del espacio de soluciones del sistema homogéneo. As´ı, teniendo en cuenta que todas las soluciones del sistema son la suma de una soluci´ on particular más una solución del sistema homogéneo, podremos resolver cualquier sistema de ecuaciones lineales (de forma aproximada en el caso de los sistema incompatibles).

Ejemplo 4. Calculemos todas las soluciones del sistema:  x + 2y + 3z = 6 x+y+z = 1  A = [1, 2, 3; 1, 1, 1], b=[6;1]  x=A\b  H=null(A)



45

Las soluciones de nuestro sistema son aquellas de la forma x + H, por ejemplo:  y = x + rand()*H Comprobémoslo, por ejemplo, escribiendo:  A*y - b Hay que tener en cuenta que todo este trabajo previo del operador \ para devolver una soluci´ on en distintas situaciones conlleva un coste de tiempo de computación.

2.7.

Gesti´ on de ficheros de hojas de c´ alculo

Como se coment´ o en el primer cap´ıtulo, vamos a mostrar el uso de una serie de funciones de MATLAB para importar/exportar datos de/a ficheros de hojas de cálculo. Las primeras csvread, csvwrite, dlmread y dlmwrite, funcionan tanto en MATLAB como Octave, mientras que las dos u ´ltimas, xlsread y xlswrite, s´ olo funcionan en MATLAB.

Ejemplo 1. Vamos a leer los datos del fichero fichero.csv que puedes descargar en http://www.editorialtebar.com.  csvread(’fichero.csv’) Si quieres puedes ver su contenido abriendo el fichero con un editor de texto cualquiera.

Como es natural, podemos asignar el contenido de una hoja de cálculo a una variable.


46


Ejemplo 2.  A = csvread(’fichero.csv’) La variable A es una matriz, y como tal admite las funciones usuales de las matrices. Por ejemplo podemos preguntarnos su tama˜ no:  size(A) Es posible que no queramos importar todos los datos de un fichero .csv. Por ejemplo, ocurre a menudo que las primeras filas o columnas corresponden a los encabezados de una tabla. En casos como este, podemos fijar a partir de qué fila y columna exportaremos los datos. Ejemplo 3. Supongamos que s´ olo queremos importar las entradas de nuestro fichero fichero.csv que se encuentran por debajo de la primera fila y a la derecha de la segunda columna, en este caso escribiremos:  csvread(’fichero.csv’,1,2) O si, por ejemplo, queremos excluir los datos de las dos primeras columnas, escribiremos:  csvread(’fichero.csv’,0,2) También podemos delimitar hasta qué fila y columna se exportarán los datos de un fichero .csv. Ejemplo 4. Para importar las entradas de nuestro fichero fichero.csv que se encuentran entre la segunda y la cuarta fila (ambas inclusive) y la primera y la tercera columnas (ambas inclusive), escribiremos:  csvread(’fichero.csv’,1,0,[1,0,3,2]) En Octave:  csvread(’fichero.csv’,[1,0,3,2])



47

La función csvwrite no entra˜ na gran complicación, con ella podemos exportar los datos de una matriz de MATLAB a un fichero .csv.

Ejemplo 5. Comencemos generando una matriz aleatoria con tres filas y cuatro columnas que llamaremos A.  A = rand(3,4) Guardemos ahora los datos de A en un fichero .csv  csvwrite(’datosA.csv’,A)

Quizá la u ńica observaci´ on a tener en cuenta es que podemos exportar nuestros datos dejando filas o columnas vac´ıas arriba o a la izquierda, respectivamente; de tal modo que luego podamos editar el fichero .csv creado para a˜ nadir nombres a las filas o columnas. Ejemplo 6. Si, por ejemplo, quisiéramos guardar el contenido de la variable A reservando espacio para a˜ nadir una fila y una columna posteriormente, bastar´ıa escribir  csvwrite(’datosA.csv’,A,1,1)

Recuerda que si lo deseas, puedes ver el contenido de los ficheros creados abriéndolos con cualquier editor de texto o importándolos a un gestor de hojas de cálculo. La funciones csvread y csvwrite son casos especiales de la funciones ´ dlmread y dlmwrite, respectivamente. Estas u ´ltimas permiten especificar el tipo de delimitador del fichero de datos (puntos y comas, espacios, ... y por supuesto comas) a diferencia de las primeras donde el delimitador es siempre una coma. Veamos un ejemplo:


48


Ejemplo 7. Definamos una matriz aleatoria con cinco filas y tres columnas  B = rand(5,3) Guardemos ahora los datos de B en un fichero .txt separando las columnas por puntos y comas (también puedes utilizar tu s´ımbolo favorito).  dlmwrite(’datosB.txt’,B,’;’) Puedes ver el contenido del fichero creado abriéndolo con cualquier editor de textos y recuperar su contenido en MATLAB con la función dlmread especificando el delimitador usado; en nuestro caso ’;’.  A=dlmread(’datosB.txt’,’;’) Si quieres puedes repetir el ejemplo anterior cambiando de delimitador. Como ya se apunt´ o anteriormente, las funciones xlsread y xlswrite solamente están implementadas en MATLAB (no en Octave). Su uso no es muy diferente del de las funciones csvread y csvwrite, aunque xlsread permite además importar el texto contenido en una hoja de cálculo, por ejemplo los nombres de las filas y columnas. Mostremos algunos ejemplos: Ejemplo 8. Descarga el fichero excel.xls de http://www. editorialtebar.com y col´ ocalo donde MATLAB pueda leerlo. Una vez descargado puedes abrir el fichero .xls para ver su contenido usando Excel u otro programa de gesti´ on de hojas de cálculo. A continuación escribe:  xlsread(’excel’) El contenido numérico del fichero .xls se mostrará en la ventana de MATLAB. Generalmente no nos interesar´ a ver el contenido del fichero. Recuerda que puedes inhibir la salida en pantalla escribiendo punto y coma ; al final de la frase.



49

Al igual que ocurr´ıa con csvread y dlmread, podemos asignar a una variable el contenido numérico de una hoja de cálculo de Excel. Ejemplo 9. En nuestro ejemplo, para guardar los datos numéricos de nuestra hoja de c´ alculo en una variable de MATLAB que llamaremos datos, basta escribir:  datos = xlsread(’excel’); También podemos importar el texto de una hoja de cálculo de Excel. Ejemplo 10. Si queremos importar el contenido numérico y el texto de nuestra hoja de c´ alculo excel.xls y guardarlo en sendas variables llamadas datos y texto, respectivamente, escribiremos:  [datos, texto] = xlsread(’excel’); La primera variable contiene los datos numéricos y la segunda el texto contenido en nuestra hoja de c´ alculo.  texto Al igual que ocurr´ıa con la funci´ on csvread podemos importar determinadas filas o columnas de una hoja de c´ alculo de Excel utilizando la función xlsread. Dado que el uso es totalmente an´ alogo al de csvread no pondremos ning´ un ejemplo ilustrativo. De igual modo, como el comando xlswrite funciona exactamente igual que csvwrite, tampoco propondremos ning´ un ejemplo para mostrar su uso. La u ńica diferencia es la salida, la función csvwrite genera un fichero de datos separados por comas y la función xlswrite una hoja de cálculo de Excel. En todo caso, animamos al lector interesado a inventarse y jugar con sus propios ejemplos.


50


2.8.

Ejercicios

Antes de comenzar escribe lo siguiente4  clear all  diary capitulo2.txt Ejercicio 1. Obtén las siguientes matrices sin escribirlas expl´ıcitamente y utilizando en menor n´ umero de comandos posible: 



1 0 0 0 0   (i)  0 1 0 0 0  0 0 1 0 0



  (ii)  

0 1 0 0

0 0 2 0

0 0 0 3

0 0 0 0

    



  (iii)  

0 0 0 0

1 0 0 1

0 1 0 2

0 0 1 3

    

Ejercicio 2. Escribe un comando MATLAB que genere una matriz 5 × 4 con valores aleatorios enteros entre −5 y 5. Ejercicio 3. Calcula el rango, el determinante y la inversa de cada una de las siguientes matrices: 

 1 2 1   A =  −1 −2 −1  2 4 2



 1 2 1   B =  −1 4 −1  2 4 2



 1 2 1   C =  −1 4 −1  2 4 5

Ejercicio 4. Construir una matriz 9 × 3 cuya primera columna sean valores entre 0 y 2π, igualmente espaciados, y las columnas segunda y tercera, el seno y el coseno de dichos valores.

Ejercicio 5. Crear una matriz cuadrada 100 × 100 de valores aleatorios y determinar cuántos de ellos son mayores que 1/2. Puesto que rand genera n´ umeros aleatorios distribuidos uniformemente en el intervalo (0,1), este n´ umero ser´ a aproximadamente 100 × 100/2 4

Si no deseas hacer un diario de la sesi´ on, no hace falta que lo escribas. En otro caso, el fichero que se genere se guaradar´ a en el directorio en el que te encuentres; recuerda que puedes saber cu´ al es escribiendo pwd .



51

Ejercicio 6. Realiza las operaciones propuestas sobre los siguientes polinomios: p(x) = 3x3 + 2x2 − x + 1

q(x) = −7x3 + 3x2 + 2x − 15 t(x) = 3x2 + 2x − 15

1. p(x) + q(x) + t(x) 2. q(x) ∗ t(x) 3.

p(x) q(x)

4. Calcular las ra´ıces de p(x) y evaluar cuánto vale el polinomio en ellas. Ejercicio 7. Obtén la soluci´ on del sistema de ecuaciones y comprueba que es correcta.  3x + 2y − z = 1   5x + y + 3z = −2  3y − 4z = 3 

Ejercicio 8. Calcula el m´ aximo autovalor y el autovector asociado (para ello, utilizar el comando eig) a la matriz de coeficientes del sistema anterior. Ejercicio 9. Descarga el fichero comprimido hipotecas.zip de http://www. editorialtebar.com y descompr´ımelo donde MATLAB pueda encontrarlo. Verás que hay seis ficheros, tres est´ an en formato .csv y otros tres en formato .xls. Los ficheros contienen los datos del n´ umero e importe total (en miles de euros) de hipotecas suscritas mensualmente entre enero de 2003 y diciembre de 2008 en Badajoz (BA), C´ aceres (CA) y Espa˜ na (ES). Se pide: 1. Guarda en sendas variables de MATLAB los datos las hipotecas en Badajoz, Cáceres y Espa˜ na.


52

Aprende MATLAB/OCTAVE en 25 horas 2. Calcula el precio medio mensual de una hipoteca en Badajoz, Cáceres y Espa˜ na dividiendo los datos de la segunda columna (importe total en miles de euros) por los de la primera (n´ umero total de hipotecas suscritas). 3. En cada caso, guarda el resultado en un vector columna y a˜ nádelo como u ´ltima columna a la variable correspondiente. 4. Exporta el contenido de las variables modificas en sendas hojas de cálculo con los nombres hipoBA, hipoCA y hipoES.

Recuerda que si usas Octave debes elegir los ficheros .csv y a˜ nadir la extensión que corresponda, en otro caso, elige el formato que prefieras para hacer el ejercicio.  diary off



2.9.

Soluciones

Ejercicio 1.  eye(3,5) ans = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0  diag([1,2,3],-1) ans = 0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0  [zeros(4,1),[eye(3);1:3]] ans = 0 1 0 0 0 0 1 0 0 0 0 1 0 1 2 3 Ejercicio 2.  round(10*rand(5,4)-5) Ejercicio 3.  A = [1 2 1; -1 -2 -1; 2 4 2] A = 1 2 1 −1 −2 −1 2 4 2  rank(A), det(A), inv(A) ans = 1 ans = 0 Warning: Matrix is singular to working precision.


53

54


ans = Inf Inf Inf Inf Inf Inf Inf Inf Inf  B = [1, 2, 1; -1, 4, -1; 2, 4, 2] B = 1 2 1 −1 4 −1 2 4 2  rank(B), det(B), inv(B) ans = 2 ans = 0 Warning: Matrix is singular to working precision. ans = Inf Inf Inf Inf Inf Inf Inf Inf Inf  C = [1, 2, 1; -1, 4, -1; 2, 4, 5] C = 1 2 1 −1 4 −1 2 4 5  rank(C), det(C), inv(C) ans = 3 ans = 18 ans = 1,3333 −0,3333 −0,3333 0,1667 0,1667 0 −0,6667 0 0,3333 Ejercicio 4.  x = linspace(0,2*pi,9)’;  D = [x, sin(x), cos(x)]


CAPÍTULO 2. MATRICES Y VECTORES D = 0 0,7854 1,5708 2,3562 3,1416 3,9270 4,7124 5,4978 6,2832

0 0,7071 1,0000 0,7071 0,0000 −0,7071 −1,0000 −0,7071 −0,0000

1,0000 0,7071 0,0000 −0,7071 −1,0000 −0,7071 −0,0000 0,7071 1,0000

Ejercicio 5.  M = rand(100);  N = M >= 1/2;  sum(sum(N)) ans = 5105 Ejercicio 6.  p = [3, 2, -1, 1]; q = [-7,3, 2, -15]; t = [3, 2, -15];  tt = [0, 3, 2, -15] tt = 0 3 2 −15  p+q+tt ans = −4 8 3 −29  conv(q,t) ans = −21 −5 117 −86 −60 225  [c,r] = deconv(p,q) c = -0.4286 r = 0 3,2857 −0,1429 −5,4286  R = roots(p) R = -1.1852


55

56


0.2593+0.4626i 0.2593-0.4626i  polyval(p,R) ans = 1.0e-014* 0.8882 0.0222+0.0888i 0.0222-0.0888i Ejercicio 7.  A = [3, 2, -1; 5, 1, 3; 0, 3, -4];  b = [1,-2,3]’;  u = A\b u = 0,1429 −0,1429 −0,8571  A*u ans = 1,0000 −2,0000 3,0000  A*u - b ans = 1.0e-015 * 0,222000 Ejercicio 8.  help eig  [V,D] = eig(A) V = −0,5515 0,3940 0,2101 −0,7949 −0,7616 −0,5160 −0,2530 −0,5146 0,8304 D =


CAPÍTULO 2. MATRICES Y VECTORES 5,4239 0 0 0 0,4402 0 0 0 −5,8641  [m,i] = max(sum(D)) m = 5.4239 i = 1  v = V(:,i) v = −0,5515 −0,7949 −0,2530 Ejercicio 9.  BA = csvread(’hipotecas BA 2003-2008.csv’); (ocultamos la salida)  CA = csvread(’hipotecas CA 2003-2008.csv’);  ES = dlmread(’hipotecas ES 2003-2008.csv’,’;’); (separado por punto y coma)  BA(:,3) = BA(:,2)./BA(:,1);  CA(:,3) = CA(:,2)./CA(:,1);  ES(:,3) = ES(:,2)./ES(:,1);  csvwrite(’hipoBA.csv’,BA)  csvwrite(’hipoCA.csv’,CA)  csvwrite(’hipoES.csv’,ES)


57

´ CAPÍTULO 3. GRAFICOS

59

Cap´ıtulo 3

Gr´ aficos Antes de comenzar con este cap´ıtulo , es importante saber que al ejecutar cualquier comando de tipo gr´ afico, se obtiene una nueva ventana llamada Figure 1 donde se dibuja el objeto o se realiza el cambio requerido. As´ı, después de ejecutar la primera orden gr´ afica, organiza el escritorio de forma que puedas ver a la vez, tanto la ventana de comandos de MATLAB/Octave como la ventana Figure 1. MATLAB tiene m´ as y mejores opciones para las cuestiones gráficas que Octave y esto se observa claramente en las diferencias existentes entre la ventana Figure 1 de MATLAB y la de Octave. Por ello, aunque se indicarán explicitamente las diferencias existentes entre MATLAB y Octave1 , sugerimos que esta sesión se realice con MATLAB.

3.1.

Creaci´ on y manipulaci´ on b´ asica de gr´ aficos

La forma m´ as sencilla de dibujar una función de una variable es mediante el uso de ezplot que tan s´ olo necesita la expresión de la propia función que se desea dibujar.

1

MATLAB indicar´ a que la orden debe ejecutarse solamente en MATLAB y Octave indicar´ a que la orden debe ejecutarse s´ olamente en Octave.


60


Ejemplo 1. Dibujemos la gr´ afica de la función f (x) = x2 + x + 1:  ezplot(’x^2+x+1’) Si el intervalo de definici´ on por defecto no parece el más adecuado, podemos cambiarlo:  ezplot(’x^2+x+1’,[-2,2]) También podemos utilizar una funci´ on anónima como argumento  f= @(x) x^2+x+1;  ezplot(f,[-2,2]) % Dar´ a error Ahora, prueba escribiendo  f= @(x) x.^2+x+1;  ezplot(f,[-2,2]) ¿Sabes por qué dio error? Si no es as´ı, no te preocupes pues lo entenderás en la siguiente secci´ on.

Dibujar una funci´ on de dos variables se puede realizar fácilmente utilizando la función ezmesh.

Ejemplo 2. Dibujemos la gr´ afica de la función 2 2 f (x, y) = x exp(−x − y ):  ezmesh(’x*exp(-x^2-y^2)’) Podemos cambiar el intervalo de definición por defecto por un dominio del tipo [−2, 2] × [−2, 2] escribiendo  ezmesh(’x*exp(-x^2-y^2)’,[-2,2]) o del tipo [−2, 2] × [−5, 7] escribiendo  ezmesh(’x*exp(-x^2-y^2)’,[-2,2,-5,7])



61

Finalmente, el argumento puede ser una función anónima  f= @(x,y) x.*exp(-x.^2-y.^2);  ezmesh(f,[-2,2,-5,7]) Ahora, NO cierres la ventana Figure 1 y contin´ ua leyendo. Es más, no cierres dicha ventana hasta que se indique expl´ıcitamente. Las funciones ezplot y ezmesh no son las u ńicas existentes en MATLAB/ Octave para dibujar gr´ aficas de funciones dadas mediante sus expresiones anal´ıticas como las que acabamos de dibujar. De hecho, existen muchas otras que incluso nos permitir´ an dibujar listas de datos dados como pares/ternas de puntos. Analizaremos algunas de las funciones más importantes de representación gráfica en la siguiente secci´ on. Cuando realizamos una representaci´ on gráfica con MATLAB/Octave lo que necesitamos es obtener una ´ optima visualización de la misma en la ventana Figure 1. Los comandos que introduciremos en lo que resta de sección nos ayudarán a conseguirla independientemente de la funci´ on utilizada para crear la gráfica. El comando axis es de gran utilidad para seleccionar el dominio y la escala de representación adecuados. Adem´ as, complementado con la función grid podemos elegir c´ omo mostrar el objeto representado. Para la representación de gráficas en tres dimensiones también contamos con la ayuda del comando view. Ejemplo 3. Uso de la funci´ on axis. Guardamos el intervalo de definici´ on de la gráfica activa  v=axis La función axis puede utilizarse tanto para modificar el dominio de definición de la funci´ on como para determinar los factores de escala de los ejes.


62


Establecemos el intervalo de definici´ on deseado para los ejes x, y, z  axis([-2,6,-20,10,-1,1]) Imponemos los l´ımites establecidos por defecto de MATLAB/Octave  axis auto Recuperamos los l´ımites originales  axis(v) Establecemos los mismos factores de escala en todos los ejes en MATLAB y en los ejes OX y OY en Octave  axis equal Restablecemos los factores de escala anteriores  axis normal Finalmente, podemos eliminar o no los ejes del dibujo  axis off  axis on

La función grid controla la rejilla que genera cada uno de los ejes y la orden view nos permite especificar el punto de observación del objeto representado en los gráficos en tres dimensiones.

Ejemplo 4.  grid off % Elimina la malla.  grid on % Dibuja una malla en la gr´ afica.  grid minor % Dibuja una malla m´ as tupida.  grid off % Elimina la malla.



63

A continuaci´ on, indicamos que el punto de observación del objeto sea de 28 grados de rotaci´ on horizontal (i.e. en el plano OXY ) con respecto al eje OY y de 90 grados de elevaci´ on vertical con respecto del plano OXY .  view(28,90) Los puntos de observaci´ on por defecto en un espacio de dos y tres dimensiones son respectivamente  view(2) % equivale a view(0,90)  view(3) % equivale a view(-37.5,30) Para conseguir una interpretaci´ on m´ as clara del objeto representado se puede a˜ nadir un t´ıtulo y nombres a los ejes. Para ello, utilizaremos las funciones title, xlabel, ylabel y zlabel. Ejemplo 5.  title (’Realizado con la funci´ on ezmesh’)  xlabel(’Eje de Abcisas’)  ylabel(’Eje de Ordenadas’)  zlabel(’Eje 3D’) Además de los comandos descritos antes para etiquetar gráficas, existe la posibilidad de escribir una anotaci´ on en cualquier otro lugar de la figura mediante el uso de la funci´ on text. Para ello, debemos escribir text(X,Y,Z,{’primera l´ ınea’,’segunda l´ ınea’}) donde X, Y, Z denotan las coordenadas del punto (X, Y, Z) de la gráfica, en este caso en tres dimensiones, donde se realiza la anotación. MATLAB/Octave por defecto interpreta todos los caracteres empezados con una barra invertida “\” como palabras clave en TeX. As´ı, podemos introducir en los gráficos caracteres especiales o letras griegas.


64


Ejemplo 6.  text(0.8,0.05,0.435,{’<- Aqu´ ı se encuentra ... aproximadamente el m´ aximo’}) ¿Te atreves a ´ındicar d´ onde se encuentra aproximadamente el m´ınimo? Pista: Quiz´ as los comandos axis y view sean de utilidad. Como ves, el uso de la funci´ on text no es sencillo. Sin embargo, en el caso dos dimensional y en MATLAB, la orden text(X,Y,{’primera l´ ınea’,’segunda l´ ınea’}), puede sustituirse por otra que utiliza la función gtext y obtener los mismos resultados de una forma m´ as sencilla. Para ello, escribimos, solamente en MATLAB, gtext({’primera l´ ınea’,’segunda l´ ınea’}) en la ventana de comandos y cuando se abra la figura indicamos mediante un clic del ratón el lugar donde ha de ir el texto. En Octave está implementado a partir de la versi´ on 3.2.. Ejemplo 7. MATLAB ezplot(’sin(x)’,[-pi,pi]) MATLAB gtext({’-\pi \leq x \leq \pi’}) Ahora selecciona en la figura el lugar donde quieres leer la anterior información con un clic del rat´ on. Para representar datos gr´ aficamente con mayor calidad necesitamos tener un mayor control sobre las opciones gr´ aficas. Las funciones get y set permiten conocer y manipular todas las opciones posibles sobre el aspecto de un gráfico independientemente de la funci´ on de MATLAB/Octave utilizada para la representación de datos. Adem´ as, su funcionamiento es sencillo. Primero, al mismo tiempo que se ejecuta el comando de dibujo para obtener la ventana



65

habitual debemos asignar un nombre o referencia al propio gráfico. Después, con la función get obtendremos informaci´ on sobre las posibles opciones de manipulación y con la funci´ on set podremos modificarlas a nuestro gusto.

Ejemplo 8. Asignamos un nombre al gr´ afico  p3=ezmesh(’x*exp(-x^2-y^2)’); Obtenemos todas las opciones posibles de control sobre el gráfico  get(p3) Supongamos que queremos cambiar el aspecto y los colores de la superficie. Obtenemos el valor actual de la opci´ on FaceColor y lo guardamos en FC  FC=get(p3,’FaceColor’) Obtenemos todos los posibles valores que puede tomar la opción FaceColor. Esta funci´ on s´ olo est´ a disponible en MATLAB. MATLAB set(p3,’FaceColor’) Establecemos el intervalo de definici´ on de los ejes para obtener una mejor visualizaci´ on de las siguientes ordenes Octave axis([-4,4,-4,4,-0.5,0.5]) A continuaci´ on, probamos los distintos valores hasta elegir el deseado. En nuestro caso recuperaremos el original.  set(p3,’FaceColor’,’none’)  set(p3,’FaceColor’,’flat’)  set(p3,’FaceColor’,’interp’)


66


 set(p3,’FaceColor’,’texturemap’)  set(p3,’FaceColor’,[0 0.5 1])  set(p3,’FaceColor’,[1 2 3]) ¿Qué ha ocurrido?¿Podr´ıas explicarlo?  set(p3,’FaceColor’,FC) Las funciones get y set son realmente muy u ´tiles pues no sólo funcionan sobre el gráfico dibujado sino también sobre cualquier acción que realicemos sobre éste. Tan s´ olo hay que recordar que al mismo tiempo que se ejecuta el comando de dibujo para obtener un resultado sobre la ventana Figure 1 debemos asignar un nombre a dicha acci´ on. As´ı, por ejemplo, si realizamos una anotación sobre el gr´ afico y, después, queremos modificar su posición, estilo de letra, tama˜ no,..... podremos hacerlo.

Ejemplo 9.  t=text(0.8,0.05,0.435,{’<- Aqu´ ı se encuentra ... aproximadamente el m´ aximo’});  get(t) % Obtenemos todas las opciones posibles de control sobre la acci´ on Supongamos que queremos cambiar el tipo de letra.  FN=get(t,’FontName’) % Obtenemos el valor actual de la opci´ on FontName y lo guardamos en FN. MATLAB set(t,’FontName’) % Intentamos obtener informaci´ on sobre los posibles valores que puede tomar la opci´ on FontName. A continuaci´ on, probamos con distintos valores hasta elegir el deseado. En nuestro caso recuperaremos el original.



67

 set(t,’FontName’,’Times’)  set(t,’FontName’,’Arial’)  set(t,’FontName’,’Cambria’)  set(t,’FontName’,FN) Ahora cambia el color y el tama˜ no de la anotación realizada en el ejemplo anterior y, si es necesario, también su posición. Finalmente, vuelve al dibujo original. En determinadas ocasiones es necesario dibujar más de una gráfica en la ventana Figure 1, para ello utilizaremos el comando hold. Además, cuando se dibuja m´ as de una gr´ afica utilizando los mismos ejes de coordenadas, es necesario poder distingir entre las mismas. Utilizando la función legend, podremos dar un nombre a cada una de las gráficas dibujadas. Ejemplo 10.  p1=ezplot(’x*exp(-x^2)’);  hold on  p2=ezplot(’x*exp(-x^2)+0.5’);  axis auto  set(p2,’Color’,[0 1 0])  title(’Uso del comando hold’)  legend(’Funci´ on x*exp(-x^2)’,’Funci´ on x*exp(-x^2)+1’)  set(p1,’Color’,[1 0 0])  hold off Sin embargo, a veces nos ser´ a de m´ as utilidad dibujar gráficas en diferentes ventanas y poder manipularlas independientemente. Para ello, utilizaremos los comandos figure, gcf, shg, clf y close. El comando figure crea una nueva ventana gráfica, la activa y la trae al


68


frente. Además, le asigna un n´ umero n y la denomina Figure n.

Ejemplo 11.  p3=ezmesh(’x*exp(-x^2-y^2)’); Octave axis([-4,4,-4,4,-0.5,0.5])  figure  p4=ezmesh(’x*exp(-x^2-y^2)+3’); Octave axis([-4,4,-4,4,2.5,3.5])  set(p4,’EdgeColor’,[0 0.5 1])  figure(4) %Elegimos el n´ umero de la nueva ventana.  p5=ezmesh(’x*exp(-x^2-y^2)-3’); Octave axis([-4,4,-4,4,-3.5,-2.5])  set(p5,’FaceColor’,[1 0.5 0])  figure  p6=ezmesh(’x*exp(-x^2-y^2)+6’); Octave axis([-4,4,-4,4,5.5,6.5])  set(p6,’EdgeColor’,[0 1 0.5])

Cuando tenemos m´ as de una ventana gráfica es conveniente saber cómo manipularlas. El comando gcf (abreviatura de “get current figure”) nos devuelve el n´ umero de la ventana gr´ afica activa en ese momento. La función shg (abreviatura de “show graph”) trae la ventana gráfica activa al frente de todas (sólamente disponible en MATLAB). Utilizando el comando clf (abreviatura de “clear figure”) podemos borrar el contenido de la ventana gráfica que necesitemos sin cerrarla. Finalmente, utilizando close cerramos las ventanas gráficas que no necesitamos.


´ ´ CAP´ÍTULO GRAFICOS CAP ITULO 3. 3. GR AFICOS

Ejemplo Ejemplo 12. 12.  gcf  gcf

MATLAB MATLAB shg shg Limpiamos la ventana activa

Limpiamos la ventana activa

 clf

 clf

 p7=ezmesh(’x*exp(-x^2-y^2)-6’);

 p7=ezmesh(’x*exp(-x^2-y^2)-6’);

 set(p7,’FaceColor’,[0.5 0 1])

 set(p7,’FaceColor’,[0.5 0 1])

Elegimos la ventana a limpiar

Elegimos la ventana a limpiar

 clf(2)

 clf(2)la ventana sobre la que queremos trabajar Activamos Activamos la ventana sobre la que queremos trabajar  figure(2)  figure(2) p8=ezplot(’x*exp(-x^2)’);  set(p8,’Color’,[1 0.5 0]) p8=ezplot(’x*exp(-x^2)’);   figure set(p8,’Color’,[1 0.5 0])   p9=ezplot(’x*exp(-x^2)’); figure  set(p9,’LineStyle’,’-.’)

 p9=ezplot(’x*exp(-x^2)’);

 figure(500)

 set(p9,’LineStyle’,’-.’)


 figure(500)

 set(p10,’LineWidth’,10)


 gcf

 set(p10,’LineWidth’,10)

Cerramos la ventana activa

 close gcf 

Cerramos la la ventana ventana Figure activa 4 Cerramos close  close(4) MATLAB close(1,3) % Cerramos las ventanas Figure 1 y Cerramos la ventana Figure 4 Figure 3

 close(4)

 close all % Cerramos todas las ventanas gr´ aficas

MATLAB close(1,3) % Cerramos las ventanas Figure 1 y Figure 3  close all % Cerramos todas las ventanas gr´ aficas © Editorial Tébar. Prohibida la reproducción sin la autorización expresa de la editorial

69 69

70


Otra opción interesante de visualizaci´ on gráfica que ofrece MATLAB es la posibilidad de dividir una ventana figura en subventanas gráficas. Para ello, se utiliza el comando subplot que permite dividir la ventana Figure en una matriz m × n de subventanas gr´ aficas. La orden para realizar esto es de la forma subplot(m,n,p) donde los dos primeros argumentos hacen referencia a las dimensiones de la matriz y el u ´ltimo argumento p que toma valores desde 1 hasta m*n hace referencia a la subventana que se quiere activar. Debemos contar las subventanas del modo siguiente: desde la primera fila a la u ´ltima y dentro de cada fila de izquierda a derecha. Ejemplo 13.  subplot(3,2,1); p1=ezmesh(’x*exp(-x^2-y^2)’);  subplot(3,2,2); p2=ezplot(’x*exp(-x^2)’);  subplot(3,2,3); p3=ezmesh(’x*exp(-x^2-y^2)+6’);  set(p3,’EdgeColor’,[0 0 1])  subplot(3,2,4); p4=ezplot(’x*exp(-x^2)’);  set(p4,’LineWidth’,10)  subplot(3,2,5); p5=ezmesh(’x*exp(-x^2-y^2)-6’);  set(p5,’EdgeColor’,[1 0 0])  subplot(3,2,6); p6=ezplot(’x*exp(-x^2)’);  set(p6,’Color’,[1 0 0],’LineWidth’,10)  close Finalmente, veamos c´ omo modificar la apariencia tanto de la ventana Figure como de los ejes de coordenadas. Para ello, debemos ampliar el uso de los comandos get y set mediante las variables gcf (abreviatura de “get current figure”) y gca (abreviatura de “get current axis”) que, como su nom-



71

bre indica, hacen referencia a la ventana Figure y a los ejes de coordenadas activos, respectivamente. Ejemplo 14.  ps=ezplot(’sin(x)’,[pi,5*pi]);  hold on  pc=ezplot(’cos(x)’,[pi,5*pi]);  hold off  set(pc,’Color’,[1 0 0])  legend(’sin’,’cos’) A continuaci´ on cambiaremos propiedades de los ejes y de la propia ventana. En primer lugar, obtenemos todas las opciones posibles de control sobre los ejes.  get(gca) Cambiamos el color de fondo, subdividimos los intervalos del eje x y cambiamos la escala del eje OX de lineal a logar´ıtmica. set(gca,’Color’,[1 1 0],’XMinorTick’,’on’,’XScale’,’log’) Obtenemos ahora todas las opciones posibles de control sobre la ventana gráfica  get(gcf) Cambiamos el color de fondo e identificamos la ventana gráfica con un nombre espec´ıfico.  set(gcf,’Color’,[0.88 1 1],’Name’,’Gr´ aficas del ... seno y del coseno’)

Seguramente, una vez dibujada y modificada a nuestro gusto la gráfica de una función, o de un conjunto de datos, necesitaremos guardarla. Para ello, podemos utilizar el comando print. Su estructura es la siguiente:


72

Aprende MATLAB/OCTAVE en 25 horas print(n,’Extensi´ on’,’Nombre’)

donde n indica el n´ umero de la ventana Figure donde se encuentra el gráfico a guardar y Extensi´ on toma el valor -djpeg o -dpsc dependiendo de si deseamos guardar la imagen como un archivo jpeg o como un archivo ps. El archivo generado se guarda en el directorio de trabajo con el nombre especificado en Nombre. Ejemplo 15.  dir  print(1,’-dpsc’,’SinCos1’) % Genera SinCos1.ps  print(1,’-djpeg’,’SinCos2’) % Genera SinCos2.jpeg  dir % Observa que se han creado dos nuevos archivos.  close Utiliza los comandos introducidos hasta ahora para mostrar la gráfica de la función x2 + x + 1 a t´ u gusto. Finalmente guárdala en formato jpeg.

3.2.

Creaci´ on especializada de gr´ aficos

Con los comandos introducidos a lo largo de la primera sección podemos realizar casi cualquier cambio que nos propongamos sobre una o varias representaciónes gráficas realizadas en una o varias ventanas Figure. Sin embargo, los comandos ezplot y ezmesh no cubren todas las posibles necesidades de reprensentación gr´ afica. En este segundo apartado, estudiaremos diversos comandos de dibujo especializados que cubrirán dicha necesidad. En MATLAB (pero no en Octave) podemos tener una visión de conjunto de todos los posibles comandos para la creación y tratamiento de gráficos.



73

Ejemplo 1. Podemos solicitar informaci´ on de tipo general, sobre gráficos 2D, sobre gráficos 3D o informaci´ on especializada MATLAB help graphics MATLAB help graph2d MATLAB help graph3d MATLAB help specgraph

En lo que resta nos centraremos en conocer los comandos básicos de creaci´ on de gráficos bidimensionales, tridimensionales y estad´ısticos.

3.2.1.

Gr´ aficos 2D

Gr´ aficos elementales

La forma m´ as sencilla de crear un gr´ afico en MATLAB/Octave es unir mediante l´ıneas unos puntos del plano dados. Para ello, los puntos del plano dados se interpretan como una lista de coordenadas (x, y) dadas y se crean dos vectores: el vector X, formado por los valores que toma la coordenada x, y el vector Y, formado por los valores que toma la coordenada y. Finalmente, la función plot une los puntos (x, y) con segmentos lineales.

Ejemplo 2. Conexi´ on de los puntos (1, 3), (4, 9) y (6, 4).  X=[1,4,6]; Y=[3,9,4];  plot(X,Y)


74


Ejemplo 3. Teclea los siguientes comandos:  hold on  px=plot(X); set(px,’Color’,[1 0 0]) ¿Qué puntos acabamos de conectar?  py=plot(Y); set(py,’Color’,[0 1 0]) Y ahora, ¿qué puntos hemos conectado? ¿Cómo funciona la función plot al hacerla actuar sobre un u ńico vector?

Ejemplo 4. Propagaci´ on del SIDA.  meses =[20,21,24,26,29,31,35,37,41,45];  casos=[110,129,220,257,439,514,878,1029,1756,2057];  plot(meses,casos)  xlabel(’Meses transcurridos desde enero de 1980’)  ylabel(’Numero de casos’)  title(’Propagaci´ on del SIDA en EE.UU.’) Notar que también podemos representar las l´ıneas que unen varios puntos para crear pol´ıgonos y adem´ as pintarlos utilizando la función fill. Ejemplo 5. Conexi´ on cerrada de los puntos (1, 3), (4, 9) y (6, 4).  X=[1,4,6,1]; Y=[3,9,4,3];  plot(X,Y) Triángulo de vértices (1, 3), (4, 9) y (6, 4) con relleno de color [1 0 1].  fill(X,Y,[1 0 1])



75

Funciones de la forma y = f (x) Como en el caso de gr´ aficas elementales, para hacer gráficas de funciones de una variable MATLAB/Octave necesita conocer primero la tabla de valores de la variable x para después dibujar la función. En particular, conocido el vector X formado por un n´ umero finito de valores de x se crea el vector imagen Y utilizando la relaci´ on y = f (x). Después, MATLAB interpreta los vectores X e Y como listas de coordenadas (x, y) de puntos sucesivos en el gráfico a representar. Finalmente, une los puntos (x, y) con segmentos lineales. Si el n´ umero de valores contenido en X es suficientemente grande, en el gráfico no apreciaremos los puntos de uni´ on. El comando ezplot crea la tabla de valores de la variable x automáticamente y, por tanto, no tenemos un control total sobre la gráfica pues no podemos elegir qué o cu´ antos puntos representar. Para conseguirlo debemos utilizar la función plot. Ejemplo 6. Dibujamos la gr´ afica de y = cos(x). Primero, creamos una tabla de valores para x con 200 valores entre 0 y 2π.  X=linspace(0,2*pi,200); Segundo, calculamos los valores de y.  Y=cos(X); Finalmente, dibujamos  plot(X,Y) Si la tabla de valores para x no tuviera una cantidad suficiente de valores obtendr´ıamos:  X=linspace(0,2*pi,10);  Y=cos(X);  plot(X,Y)


76

Aprende MATLAB/OCTAVE en 25 horas Ahora, veamos un ejemplo un poco m´ as complicado: 2

Ejemplo 7. Dibujamos la gr´ afica de y = xe−x .  X=-3:.01:3;  Y=X.*exp(-X.^2); ¿Por qué hay que poner los puntos antes de las operaciones?  plot(X,Y) ¿Sabr´ıas responder ahora a la pregunta formulada en el Ejemplo 1 de la primera secci´ on? Algunas funciones gr´ aficas de MATLAB, como la que actualmente estudiamos plot, nos permiten simplificar el procedimiento de realización de cambios sobre el aspecto de la gr´ afica. Es decir, a veces podremos evitar el uso de hold y, para ciertas propiedades podremos evitar el uso de las funciones get y set. El uso del comando help nos permitirá conocer las propiedades de simplificación que ofrece cada funci´ on de dibujo gráfico. Ejemplo 8. Observar como la funci´ on plot ofrece posibilidades, por ejemplo sobre el color o tipo de linea, que no ofrece ezplot.  help plot  help ezplot As´ı, para cambiar el color o tipo de l´ınea con ezplot es necesario el uso de set mientras que su uso puede ser evitado en el caso de dibujar con la función plot. A continuaci´ on, estudiaremos las posibilidades de simplificación más relevantes que ofrece el comando plot con más detalle. A˜ nadiendo un par´ ametro opcional al comando plot podemos controlar el color y el estilo de l´ıneas y puntos. Las opciones predefinidas para los colores son: ’b’ (azul), ’g’ (verde), ’r’ (rojo), ’c’ (cian), ’m’ (magenta), ’y’ (amarillo), ’k’ (negro) y ’w’ (blanco). Las



77

opciones predefinidas de estilo de l´ınea son: ’-’ (l´ınea continua), ’:’ (l´ınea de puntos), ’-.’ (l´ınea que alterna el gui´ on con el punto) y ’- -’ (l´ınea de guiones). Ejemplo 9. Dibujamos la gr´ afica anterior eligiendo el color magenta y el tipo de l´ınea que alterna el gui´ on con el punto.  plot(X,Y,’m-.’) ¿Sabr´ıas colorearla utilizando un color no predefinido?

También podemos elegir un marcador para destacar los puntos (x, y) en la gráfica. Es m´ as, si no nos interesa que los puntos de la gráfica se visualicen conectados por segmentos lineales, tan sólo debemos especificar el tipo de marcador que deseamos para pintar los puntos y no decir nada sobre el tipo de l´ınea. Los tipos de puntos predefinidos son: ’.’ (puntos), ’o’ (c´ırculos), ’x’ (aspas), ’+’ (sumas), ’*’ (estrellas), ’s’ (cuadrados), ’d’ (diamantes), ’p’ (pentágonos), ’h’ (hex´ agonos) y ’v’,’ˆ’,’<’,’>’ (triángulos con orientación distinta). Ejemplo 10. Dibujamos la gr´ afica del coseno eligiendo el color verde y el tipo de punto diamante. Primero, con l´ıneas continuas que conectan los puntos y luego sin ellas.  X=linspace(0,2*pi,20);  Y=cos(X);  plot(X,Y)  plot(X,Y,’gd-’)  plot(X,Y,’gd’) Ahora mostramos una aplicaci´ on de su utilidad:


78


Ejemplo 11. Representaci´ on gr´ afica de la influencia que tienen las embarcaciones a motor en la muerte de los manat´ıes.  matriculas=[13 21 24 16 24 20 15 34 33 33 39 43 50 47];  muertes = [447 460 481 498 513 512 526 559 585 614 ... 645 675 711 719];  plot(matriculas,muertes, ’bo’)  xlabel(’Embarcaciones matriculadas (en miles) desde ... 1977’)  ylabel(’Manat´ ıes hallados muertos’)  title(’Mortalidad de los manat´ ıes en Florida’) Finalmente, una de las propiedades de simplificación más interesantes de plot es el poder evitar el uso del comando hold. Ejemplo 12.  X=0:.5:20;  Y=exp(0.1*X);  Y1=Y.*sin(X); Y2=Y.*cos(X);  plot(X,Y,X,-Y,X,Y1,X,Y2)  plot(X,Y,’kv’,X,-Y,’k^’,X,Y1,’ms:’,X,Y2,’co-’) Sin embargo, cuando utilicemos esta u ´ltima opción debemos llevar cuidado pues el comando plot tiene otras propiedades de simplificación gráfica espec´ıficas sobre el aspecto de la l´ınea. En el caso de MATLAB estas se asignarán a todas las l´ıneas del gr´ afico y en el caso de Octave sólo a algunas de ellas si no se indican correctamente. Estas otras propiedades se especifican mediante parejas de nombre de la propiedad y su valor como se muestra en el siguiente



79

ejemplo.

Ejemplo 13.  plot(X,Y1,’ms:’,X,Y2,’co-’) Ahora, fijamos el grosor de la l´ınea en p´ıxeles.  plot(X,Y1,’ms:’,X,Y2,’co-’,’LineWidth’,2) Además, fijamos el tama˜ no del marcador.  plot(X,Y1,’ms:’,X,Y2,’co-’,’LineWidth’,2, ... ’MarkerSize’, 10) ¡Notar que en MATLAB las propiedades especificadas afectan a todas las l´ıneas mientras que en Octave solo afectan a la gráfica de Y2! ¿Sabes qué hacer para que cada l´ınea tenga su propio grosor y ambas se dibujen sobre los mismos ejes de coordenadas?

En ciertas ocasiones nos puede resultar de utilidad pintar la gráfica de una función o de unos datos dados en escala logar´ıtmica. En el Ejemplo 14 de la primera sección se muestra c´ omo podemos conseguirlo utilizando la función set; Sin embargo, las funciones semilogx, semilogy y loglog, nos permiten obtener la representaci´ on gr´ afica directamente con escala logar´ıtmica en el eje OX, en el eje OY o en ambos, respectivamente.

Ejemplo 14. Descarga el archivo Ejemplo14.m que está en http: //www.editorialtebar.com y gu´ ardalo en el directorio de trabajo. En este archivo se encuentran ciertos datos obtenidos de forma experimental. Cargamos los datos  Ejemplo14


80


Representamos gr´ aficamente los resultados con escala lineal en ambos ejes de coordenadas.  subplot(1,2,1)  plot(diametros,tasas,’ro’)  title(’Grafica dN/dp versus p’)  xlabel(’Diametros p de las particulas’)  ylabel(’Densidad del tama~ no dN/dp’) Representamos gr´ aficamente los resultados con escala lineal en el eje OX y escala logar´ıtmica en el eje OY .  subplot(1,2,2)  semilogy(diametros,tasas,’b*’)  title(’Grafica SemiLog: dN/dp versus p’)  xlabel(’Diametros p de las particulas’)  ylabel(’Logaritmo de la densidad del tama~ no dN/dp’)

Curvas en param´ etricas

Veamos ahora c´ omo se pueden representar curvas en el plano dadas en forma paramétrica, es decir, de la forma r(t) = (x(t), y(t)), t ∈ [a, b]. Para ello, se utilizan principalmente los comandos plot y ezplot. Ejemplo 15. Dibujamos la circunferencia unidad.  subplot(1,2,1)  ezc=ezplot(’cos(t)’,’sin(t)’);



81

 set(ezc,’color’,[1 0 1])  axis square  subplot(1,2,2)  t=0:.01:1;  plot(cos(2*pi*t),sin(2*pi*t),’g’)  axis square Compara ambas gr´ aficas. ¿C´ omo es m´ as sencillo elegir el color de la gráfica? ¿Qué otras diferencias visuales observas? Si lo has realizado con MATLAB, ¿d´ onde resulta m´ as necesario forzar la misma escala en ambos ejes de coordenadas? Notar que hemos utilizado “;” después de la sentencia de ezplot. Sin embargo, la ventana gr´ afica se ha abierto igualmente. El uso de “;” tan sólo suprime el texto de salida en la ventana de comandos. Una particularidad de las curvas en paramétricas es que se pueden dibujar sobre ellas los vectores velocidad, utilizando el comando quiver cuya sintaxis para este caso es quiver(r(t),r’(t)) donde r (t) denota la derivada de r(t) respecto de t. Ejemplo 16.  hold on  s=linspace(0,2*pi,20);  quiver(cos(s),sin(s),-sin(s),cos(s))  axis([-1.5,1.5,-1.5,1.5]) Notar que es mejor elegir un n´ umero reducido de puntos con el comando linspace, pues este ser´ a el n´ umero de vectores que se dibujen.


82


Finalmente, también podemos saber la dirección en que se pinta la curva utilizando el comando comet. Ejemplo 17.  clf MATLAB t=linspace(-5,5,1000); Octave t=linspace(-5,5,100); Ahora, aseg´ urate de ver la ventana Figure a la vez que ejecutas la orden siguiente:  comet((t.*(t.^2-1))./(t.^2+1),(2*(t.^2-1))./(t.^2+1)) Notar que la velocidad de ejecuci´ on depende del n´ umero de puntos que hayamos generado con el comando linspace.

3.2.2.

Gr´ aficos 3D

Curvas en el espacio

El comando m´ as sencillo para dibujar curvas en el espacio es el ezplot3. Sin embargo, al igual que ocurre con los comandos ezplot y plot, si utilizamos el comando plot3 tendremos un mayor control sobre la curva.

Ejemplo 18.  close all  subplot(1,2,1)  ezplot3(’cos(2*pi*t)’,’sin(2*pi*t)’,’t’)



83

 subplot(1,2,2)  T=-2:0.01:2;  plot3(cos(2*pi*T),sin(2*pi*T),T,’r’) ¿Podr´ıas conseguir que ambas gr´ aficas sean exactamente iguales? Pista: Comienza de la forma siguiente  subplot(1,2,1)  ezp=ezplot3(’cos(2*pi*t)’,’sin(2*pi*t)’,’t’,[-2,2]);  set(ezp,’color’,’r’)

Funciones de la forma z = f (x, y)

Para dibujar gr´ aficos de funciones de dos variables, al igual que para funciones de una variable, en primer lugar hay que generar tablas de valores para las variables x e y. en realidad, ahora lo que tenemos que hacer es generar una malla sobre un rect´ angulo del plano OXY . Para ello se utiliza el comando meshgrid. Una vez generada la malla hay dos comandos básicos para dibujar superfices en el espacio tres dimensional: mesh y surf. Su u ńica diferencia es que con mesh obtenemos una superficie mallada transparente y con surf la superficie mallada es opaca. Además, vimos que la forma m´ as sencilla de dibujar una superficie era con la función ezmesh ya que la malla se genera de forma automática. Por supuesto, existe el comando ezsurf cuya u ńica diferencia es que la superficie se dibuja opaca en lugar de transparente.


84


Ejemplo 19. Gr´ afica de la funci´ on z = x2 − y 2 en la región [−2, 2] × [−2, 2]  clf  subplot(2,2,2)  [x,y]=meshgrid(-2:0.5:2);  z=x.^2-y.^2;  mesh(x,y,z)  subplot(2,2,3)  ezsurf(’x^2-y^2’,[-2,2]) Ahora cambiamos la regi´ on por [−2, 2] × [−4, 4]  subplot(2,2,1)  ezmesh(’x^2-y^2’,[-2,2,-4,4])  subplot(2,2,4)  [X,Y]=meshgrid(-2:0.5:2,-4:0.3:4);  Z=X.^2-Y.^2;  surf(X,Y,Z) Nota: Si no observas bien las diferencias maximiza la ventana Figure.

Si f (x, y) es un n´ umero, el procedimiento [X,Y]=meshgrid(x,y); Z=f(X,Y); no define una matriz Z con las mismas dimensiones que X e Y. Esto se resuelve escribiendo Z=f(X,Y)+0.;



85

Por ejemplo, cuando queremos representar la gráfica de un plano perpendicular al eje z Ejemplo 20. Gr´ afica de la funci´ on z = 3 en la región [−5, 5] × [−5, 5]  clf  [X,Y]=meshgrid(-5:0.1:5);  Z=3+0.*X;  surf(X,Y,Z)  close

3.2.3.

Gr´ aficos Estad´ısticos

A pesar de que no se puede decir que MATLAB sea el programa ideal para hacer cálculos relacionados con la Estad´ıstica, dispone de algunos comandos para generar gr´ aficos estad´ısticos como diagrama de sectores o barras e histogramas. Consideremos como ejemplo el n´ umero de estudiantes titulados en ciencias y letras en cinco universidades: Universidades

Titulados en Ciencias

Titulados en Letras

A B C D E

100 60 90 20 40

90 120 30 125 80

Diagramas de sectores El diagrama de sectores permite comparar las diversas opciones en un c´ırculo con sectores cuyo ´ angulo es directamente proporcional al porcentaje de cada opción. Los realizaremos utilizando el comando pie. Para obtener


86


una visión en tres dimensiones en MATLAB (pero no en Octave) utilizamos pie3.

Ejemplo 21.  x=[100 60 90 20 40];  pie(x)  legend(’A’,’B’,’C’,’D’,’E’) Nota: Si la leyenda no aparece en un lugar adecuado en MATLAB (pero no en Octave) podemos moverla utilizando el botón izquierdo del ratón.  title(’Numero de estudiantes titulados en Ciencias’) MATLAB pie3(x) MATLAB legend(’A’,’B’,’C’,’D’,’E’) MATLAB title(’\bf Numero de estudiantes titulados en ... Ciencias’)

Además, tanto para el comando pie, como para el comando pie3 existe la posibilidad de resaltar porciones.

Ejemplo 22.  y=[1 0 0 1 0]; %Indicamos los datos a resaltar.  pie(x,y)  legend(’A’,’B’,’C’,’D’,’E’)  title(’Numero de estudiantes titulados en Ciencias’)



87

MATLAB pie3(x,y) MATLAB legend(’A’,’B’,’C’,’D’,’E’) MATLAB title(’\bf Numero de estudiantes titulados en ... Ciencias’) Representa mediante un diagrama de sectores el n´ umero de estudiantes titulados en Letras por la Universidad XX. Además, compáralos con los titulados en Ciencias utilizando el comando subplot.

Diagramas de barras

El diagrama de barras permite comparar las diversas opciones mediante barras cuya altura es directamente proporcional al porcentaje de cada opción. Existen varias posibilidades para representar diagramas de barras, dependiendo de si queremos que las barras se representen de forma vertical u horizontal y en dos o tres dimensiones. Los comandos que nos permiten realizar estos gráficos son: bar, barh, bar3 y bar3h. Al igual que para representar diagramas de sectores, la visi´ on en tres dimensiones de los diagramas de barras funciona en MATLAB y no en Octave. Ejemplo 23.  x = [’A’ ’B’ ’C’ ’D’ ’E’];  y = [100 60 90 20 40]; MATLAB subplot(2,2,1) Octave subplot(1,2,1)  bar(y)  title(’Barras Verticales-Titulados en Ciencias’)


88


 set(gca,’xtick’,1:length(x));  set(gca,’xticklabel’,x); MATLAB subplot(2,2,2) Octave subplot(1,2,2)  bar(y)  title(’Barras Horizontales-Titulados en Ciencias’)  set(gca,’xtick’,1:length(x));  set(gca,’xticklabel’,x); MATLAB subplot(2,2,3)  bar3(y)  title(’Barras Verticales 3D-Titulados en Ciencias’)  set(gca,’xtick’,1:length(x));  set(gca,’xticklabel’,x); MATLAB subplot(2,2,4)  bar3(y)  title(’Barras Horizontales 3D-Titulados en Ciencias’)  set(gca,’xtick’,1:length(x));  set(gca,’xticklabel’,x); Además, podemos representar datos agrupados de manera que resulte fácil compararlos. Ejemplo 24. Comparaci´ on entre Titulados por Ciencias y por Letras.  clf  y=[100 90;60 120;90 30;20 125;40 80];



89

MATLAB subplot(2,2,1) Octave subplot(1,2,1)  bar(y), title(’Barras Verticales’)  legend(’Titulados en Ciencias’,’Titulados en Letras’) MATLAB subplot(2,2,2) Octave subplot(1,2,2)  barh(y), title(’Barras Horizontales’)  legend(’Titulados en Ciencias’,’Titulados en Letras’) MATLAB subplot(2,2,3) MATLAB bar3(y), title(’Barras Verticales 3D’) MATLAB legend(’Titulados en Ciencias’,’Titulados en ... Letras’) MATLAB subplot(2,2,4) MATLAB bar3h(x,y), title(’Barras Horizontales 3D’) MATLAB legend(’Titulados en Ciencias’,’Titulados en ... Letras’) Averigua cómo modificar el ancho de las barras y cámbialo a tu gusto en los gráficos del ejemplo anterior.

Histograma

El histograma es una representaci´ on gráfica de una distribución de frecuencias por medio de rect´ angulos, cuyas anchuras representan intervalos de valores y sus alturas las correspondientes frecuencias. Es el u ńico gráfico de los estudiados en esta secci´ on que resulta adecuado para representar datos cuantitativos con muchos valores distintos. Se genera con el comando hist.


90


Ejemplo 25. Si generamos 10.000 n´ umeros aleatorios siguiendo una distribución normal est´ andar  x=randn(10000,1) al ejecutar la orden  hist(x) obtenemos un histograma con 10 intervalos. Si aumentamos el n´ umero de intervalos el histograma se parece más a la curva normal  hist(x,100)



3.3.

91

Ejercicios

Antes de comenzar escribe lo siguiente2  clear all  diary capitulo3.txt Ejercicio 1. Dibuja juntas las gr´ aficas de y = x2 − 1 y de y = −x2 + 1 en rojo y azul respectivamente. Ejercicio 2. Dibuja un cuadradado con el borde rojo y rellénalo de verde. Ejercicio 3. Dibuja de forma adecuada un cuadrado inscrito en una circunferencia y preséntalo a t´ u gusto. 2

Ejercicio 4. Representa la gr´ afica de f (x) = 3xe−x como más te guste. Ahora, usando el comando help aprende cómo usa MATLAB el comando ginput. Finalmente, u ´salo para indicar en la gráfica de f dónde se encuentran el máximo y el m´ınimo. Ejercicio 5. Dibuja la curva en el espacio r(t) = (sin(t), cos(t), t). Además, dibuja los vectores velocidad utilizando el comando quiver3. Utiliza la ayuda para descubrir c´ omo funciona y ten en cuenta lo aprendido sobre el comando quiver. NOTA: Aseg´ urate que todos los vectores que defines tienen la misma dimensi´ on. Te puede ayudar el comentario en el u ´ltimo párrafo de la sección 3.2.2. Ejercicio 6. Representa la gr´ afica de la función de dos variables z = sin(x)cos(y) de forma adecuada. Ahora, descubre la función colormap y util´ızala. Ejercicio 7. Elige los comandos gr´ aficos de MATLAB adecuados para estudiar la producci´ on de cereales en Espa˜ na, Alemania y Austria. Los datos son, por supuesto ficticios. los siguientes: 2

Si no deseas hacer un diario de la sesi´ on, no hace falta que lo escribas. En otro caso, el fichero que se genere se guaradar´ a en el directorio en el que te encuentres; recuerda que puedes saber cu´ al es escribiendo pwd.


92

Aprende MATLAB/OCTAVE en 25 horas A˜ no

Espa˜ na

Alemania

Austria

trigo cebada avena arroz ma´ız

100 50 90 20 40

90 120 30 125 80

6 13 2 15 9

Ejercicio 8. Descarga de la direcci´ on http://www.editorialtebar.com los archivos hipoBA, hipoCA y hipoES que generamos en la sesión anterior. Dibuja en una misma gr´ afica la evoluci´ on mensual entre Enero de 2003 y Diciembre de 2008 del precio medio de las hipotecas en Badajoz, Cáceres y Espa˜ na. Ejercicio 9. Dibuja la gr´ afica de la funci´ on  2  si x < 0  x f (x) = 1 si 0 ≤ x ≤ 1   x + 2 si 1 < x

Notar que es muy importante elegir una tabla de valores y una representación adecuadas para que aparezcan los aspectos más representativos de la función. En este caso, las discontinuidades.  diary off



3.4.

Soluciones

Ejercicio 1.  p1=ezplot(’1-x^2’);  hold on  p2=ezplot(’x^2-1’);  axis auto  set(p2,’Color’,[1 0 0])  hold off Ejercicio 2.  X=[-1,1,1,-1,-1];Y=[-1,-1,1,1,-1];  fill(X,Y,[0,1,0])  hold on  plot(X,Y,’r-’)  axis( [-2,2,-2,2])  hold off Ejercicio 3.  P = [1,0;0,1;-1,0;0,-1;1,0];  plot(P(:,1), P(:,2))  hold on  t = linspace(0,2*pi);  plot(sin(t),cos(t))  hold off Ejercicio 4.  ezplot(’3*x*exp(-x^2)’)  axis([-8,8,-2,2])  help ginput  [X,Y] = ginput(2);  text(X(1),Y(1),’<-- M´ aximo’)  text(X(2),Y(2),’<-- M´ ınimo’) Ejercicio 5.  t = linspace(0,2*pi,20);


93

94    

Aprende MATLAB/OCTAVE en 25 horas plot3(sin(t),cos(t),t,’r’) hold on quiver3(sin(t),cos(t),t,cos(t),-sin(t),1+0*t) hold off

Ejercicio 6.  ezmesh(’sin(x)*cos(y)’)  help colormap  colormap(cool)  colormap(summer) Ejercicio 7.  A = [100, 90, 6; 50, 120, 13; 90, 30, 2; 20, 125, 15; 40, 80, 9]  plot(A(:,1),A(:,2))  hold on  plot(A(:,3),’r’)  plot(,A(:,4),’g’)  hold off Otra forma:  subplot(1,3,1); bar(A(:,2))  subplot(1,3,2); bar(A(:,3))  subplot(1,3,3); bar(A(:,4))  close all Ejercicio 8.  BA = csvread(’hipoBA.csv’);  CA = csvread(’hipoCA.csv’);  ES = csvread(’hipoES.csv’);  plot(BA(:,3))  hold on  plot(CA(:,3),’r’)  plot(ES(:,3),’g’)  hold off Ejercicio 9.  ezplot(’x^2’,[-10,0])


´ CAPÍTULO 3. GRAFICOS     

hold on ezplot(’1’,[0,1]) ezplot(’x+2’,[1,10]) axis([-5,5,-1,10]) hold off


95

´ A LA PROGRAMACION ´ CAPÍTULO 4. INTRODUCCION

97

Cap´ıtulo 4

Introducci´ on a la programaci´ on Al igual que los anteriores, este cap´ıtulo está pensado para ser trabajado delante de un ordenador con MATLAB u Octave instalado. Sin embargo a diferencia de los cap´ıtulos anteriores, no veremos comandos, sino cómo crear ficheros que nos permitan automatizar el trabajo que queremos hacer. Para ello utilizaremos el editor de MATLAB (se accede a él a través de “Archivo>Nuevo” o pulsando el icono de nuevo documento). En el caso de Octave, cualquier editor nos servir´ a, por ejemplo gedit o notepad.

4.1.

“Scripts”

Al ser este un curso introductorio, no se asume que el alumno tenga conocimientos de programaci´ on. Para aquellos lectores que los tengan, recomendamos que miren la sintaxis de los distintos comandos de las tres primeras secciones y comiencen por la secci´ on 3.2.3. Hasta ahora hemos visto c´ omo utilizar MATLAB como una calculadora, las funciones que soporta, el manejo de los vectores y matrices y la representación gráfica. En esta sesi´ on aprenderemos c´ omo organizar el trabajo de modo que el ordenador realice toda una serie de tareas sin tener que ir introduciendo los comandos uno a uno. Un archivo de comandos es una serie de órdenes que se escriben en un ar-


98


chivo para usarlas siempre del mismo modo (por ejemplo, si queremos resolver un problema complejo para varios conjuntos diferentes de datos sin tener que escribir cada vez todas las ´ ordenes). MATLAB admite dos tipos de archivos de comandos: “scripts” y funciones. En esta sección veremos los “scripts” y en la siguiente las funciones. Ambos tipos de archivos se pueden crear bien con el editor de MATLAB ( en “Archivo > Nuevo”) o con un editor de texto cualquiera. En todo caso hay que a˜ nadir al final la extensión “.m”. Los archivos “.m” han de guardarse donde el programa pueda encontrarlos. Esto puede ser en el directorio de trabajo o en cualquier directorio del “path” (los directorios también se llaman carpetas). Puedes consultar cuál es el actual directorio con la orden pwd. Los directorios que están actualmente en el path se pueden consultar con path y para a˜ nadir un nuevo directorio se utiliza addpath. También pueden ambos modificarse desde el men´ u “Archivo” de MATLAB. Ejemplo 1. Crea un directorio para trabajar y a˜ nádelo al path.  mkdir(’prueba’)  addpath(’prueba’) Los “scripts” nos permiten agrupar un conjunto de órdenes de MATLAB y ejecutarlas fácilmente. El efecto de teclear el nombre del “script” en la ventana de comandos es ex´ actamente el mismo que si tecleásemos directamente cada una de las l´ıneas del fichero. Ejemplo 2. Abre un nuevo fichero. Escribe en él el siguiente texto: x=0:0.01:2; y=cos(x).*sin(x); plot(x,y) Guárdalo como “ejemplo1.m” en el directorio que has creado antes. Ahora para ver lo que hace, vuelve a la consola de MATLAB y teclea  ejemplo1 Una ventaja de realizar los c´ alculos mediante “scripts” en lugar de directamente en consola, es que se pueden realizar modificaciones (por ejemplo, si



99

nos hemos confundido) y volver a ejecutar. Ejemplo 3. Vamos a modificar el ejemplo de modo que nos represente otra función. Para ello volvemos al editor, abrimos “ejemplo1.m” y lo modificamos para que ahora ponga x=0:0.01:2*pi; y=cos(x).*sin(x); plot(x,y) Ahora para ver lo que hace, vuelve a la consola de MATLAB y teclea  ejemplo1 Los “scripts” de MATLAB comparten variables con la consola. Es decir, el valor de cada variable al comenzar el script es el mismo que tuviese en la consola. Si se modifica una variable (cambia el valor, se crea o se destruye), lo mismo pasa en la consola. Ejemplo 4. Vuelve al editor, teclea las siguientes l´ıneas y guárdalo como “ejemplo2.m”. a=2; d=a+c; b=d+c; Teclea en consola de MATLAB  who  a=5  c=5  b=5  ejemplo2  who  a  b  c=3  ejemplo2  b Esto hace que los ”scripts”sean muy poco prácticos y en general deber´ıan


100


utilizarse funciones. El u ńico caso en el que podr´ıa ser razonable utilizar “scripts” en lugar de funciones es cuando queremos realizar una tarea peque˜ na que no vamos a volver a hacer.

4.2.

Funciones

Como ya hemos visto, los “scripts” no son apropiados para escribir programas largos pues comparten variables con la consola. Para poder tener un mayor control sobre las variables y poder reutilizar el código que escribamos, utilizaremos las funciones. Las funciones son similares a los “scripts” salvo por las siguientes caracter´ısticas: La primera l´ınea de una funci´ on debe comenzar por function y la u ´ltima debe ser end.

Todos los datos que recibe o devuelve una función tienen que ser detallados en la primera l´ınea.

Las variables de una funci´ on son locales. Esto es, sólo existen mientras se está ejecutando la funci´ on. Es m´ as, aunque una variable de la consola tenga el mismo nombre que una de la función, esto no afectará a la función ni a la variable de la consola.

El formato para declarar la funci´ on es el siguiente:

function vs=nombre funcion(ve1,ve2,...,ven)

donde vs es la variable de salida, que es donde se guardará el resultado de la función, ve1, ve2,. . . , ven, son las variables de entrada (los datos que recibirá la función) y nombre f uncion es el nombre del archivo “.m” donde se guardará (sin la extensi´ on).



101

Ejemplo 1. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “sumacuadrados.m”. function c=sumacuadrados(a,b) c=a^2+b^2; end Vuelve a la consola de MATLAB y ejecuta la función que acabas de crear.  cuadrado=sumacuadrados(3,4) También podemos utilizar las funciones directamente, sin asignar el contenido a ninguna variable.  sqrt(sumacuadrados(3,4)) MATLAB también permite funciones que devuelvan más de un valor. La sintaxis de esas funciones es la siguiente: function [vs1,vs2,...,vsm]=nombre funcion(ve1,ve2,...,ven)

Ejemplo 2. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “mediayvarianza.m”. function [media,varianza]=mediayvarianza(v) media=mean(v); varianza=var(v); end Al utilizar esta funci´ on desde la consola, debemos guardar los dos valores que devuelve en dos variables.  [a,b]=mediayvarianza([4,3,5,7]) También se puede utilizar en c´ alculos, aunque en ese caso el valor de la función ser´ a el de la primera variable de salida.  mediayvarianza([4,3,5,7])+1 Nótese que las variables de una funci´ on pueden ser de cualquiera de los tipos que usa MATLAB.


102


Ejemplo 3. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “dedosendos.m”. function v=dedosendos(n,m) v=n:2:m; end Ahora, volvemos a la consola y creamos un vector con los n´ umeros impares de 1 al 30.  dedosendos(1,30)

Todas las l´ıneas de comentario que se escriban a continuación de la primera l´ınea serán la documentaci´ on de la funci´ on, es decir, lo que nos aparecerá cuando utilicemos help seguido del nombre de la función.

Ejemplo 4. Vamos a escribir la documentación de la función del ejemplo anterior, “dedosendos.m”. function v=dedosendos(n,m) % dedosendos(n,m) % Devuelve un vector con los n´ umeros % n, n+2, ... n+2k, ... tales que n+2k<=m. v=n:2:m; end Ahora, volvemos a la consola y pedimos información sobre la función.  help dedosendos

En general MATLAB est´ a pensado para que los argumentos de entrada y salida sean vectores, aunque en ocasiones las funciones se programan pensando que son escalares. Para aplicar una funci´ on programada para recibir escalares a un vector (componente a componente), existe la orden arrayfun.



103

Ejemplo 5. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “minimoconvexa.m”. function y=minimoconvexa(a) % minimoconvexa(a) % Calcula el m´ ınimo (aproximado) de (1-a)sin(t)+a cos(t) % para t entre 0 y 2pi y=min((1-a)*sin(0:0.01:2*pi)+a*cos(0:0.01:2*pi)); end Ahora vamos a calcular el m´ınimo de (1 − a) sin(t) + a cos(t) para t ∈ (0, 2π) y a ∈ [0, 1] (aproximado). Ve a la consola y escribe.  a=0:0.01:1;  min(arrayfun(@minimoconvexa,a))

Las funciones pueden llamar a otras funciones (o incluso a s´ı mismas, como veremos m´ as adelante). Esto permite descomponer los problemas en subproblemas y resolver cada uno de ellos mediante una función.

Ejemplo 6. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “distanciaorigen.m”. function d=distanciaorigen(x0,y0) % distanciaorigen(x0,y0) % Calcula la distancia de un punto (x0,y0) al origen % Si x0,y0 son vectores, devuelve un vector que % en cada posici´ on k tiene la distancia de % (x0(k),y0(k)) al origen. d=sqrt(x0.^2+y0.^2); end


104


Guarda la siguiente funci´ on como “distancia.m”. function d=distancia(x,y,x0,y0) % distancia(x,y,x0,y0) % Calcula la distancia de un punto (x,y) a otro (x0,y0) % Si x,y son vectores, devuelve un vector % que en cada posici´ on k tiene la distancia % de (x(k),y(k)) a (x0,y0). % Si x,y,x0,y0 son vectores, devuelve un vector % que en cada posici´ on k tiene la distancia % de (x(k),y(k)) a (x0(k),y0(k)). d=distanciaorigen(x-x0,y-y0); end Vuelve a la consola de MATLAB y prueba la función que acabas de crear para calcular la distancia de pares de puntos al (1, 1).  x=[0,1,3,5];y=[3,2,1,0]  distancia(x,y,1,1) Todas las funciones aqu´ı mostradas reciben los datos a partir de las variables de entrada y los devuelven en las variables de salida. Nada impide que una función lea datos de un archivo o los pida al usuario o que los imprima por pantalla o represente en una gr´ afica. Sin embargo, para que las funciones sean fácilmente reutilizables y sencillas de depurar lo más conveniente siempre es separar las funciones en tres tipos: las que reciben datos, las que los manipulan (que ser´ıan del tipo que hemos visto) y las que los devuelven. Para enlazar estos tres tipos de funciones se puede hacer en consola, recurrir a los scripts vistos en el primer tema o definir una función que no tenga variables de entrada ni de salida, como veremos en el apartado de entrada y salida.

4.3.

Estructuras de control

Hasta ahora hemos visto c´ omo ejecutar una lista de comandos de MATLAB. En muchas ocasiones nos interesar´ a que la ejecución de los comandos no sea lineal, sino que las instrucciones se ejecuten seg´ un ciertos parámetros, resultados anteriores, etc.



105

Las ordenes que permiten variar la ejecución de los comandos (también llamado flujo del programa) se denominan estructuras de control. Hay dos tipos fundamentales de estructuras de control: condicional y repetición. En el primer caso el siguiente comando a ejecutarse dependerá de si una condición es cierta o falta. En el segundo caso, un mismo grupo de instrucciones se ejecuta repetidamente, bien un n´ umero determinado de veces o bien mientras se cumpla una condici´ on. Comencemos por las estructuras de control condicionales. La principal estructura de control de este tipo es if. Vamos a programar la función valor absoluto utilizando esta estructura de control.

Ejemplo 1. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “abs1.m”. function absx=abs1(x) % abs1(x) % Devuelve el valor absoluto de x absx=x; if x<=0 absx=-absx; end end Ahora, volvemos a la consola y probamos la función.  abs1(3)  abs1(-3)

Como se ha visto en el ejemplo, detr´ as de if debe ir una condición. Si el resultado de evaluar esa condici´ on es verdadero, se ejecutarán las sentencias entre el if y el end. Si el resultado es falso, se continuará detrás del end. Si en el caso en el que la condici´ on del if es falsa queremos que se ejecuten otras l´ıneas adicionales, usaremos la palabra reservada else.


106


Ejemplo 2. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “abs2.m”. function absx=abs2(x) % abs2(x) % Devuelve el valor absoluto de x if x>0 absx=x; else absx=-x; end end Ahora, volvemos a la consola y probamos la función.  abs2(3)  abs2(-3) Cuando necesites usar if en MATLAB, piensa si puedes sustituirlo usando los operadores de comparaci´ on y los l´ ogicos. En general, esto te permitirá escribir programas m´ as compactos e incluso más sencillos de leer. Por ejemplo, para obtener el valor absoluto de una variable x (sin usar abs) se puede escribir (x>0)*x+(~ (x>0))*(-x). Un uso especial de if es en las funciones recursivas. Una función recursiva es aquella que aparece en su definici´ on. El ejemplo más clásico es el factorial de un n´ umero n, definido como n! = n∗(n−1)!, 0! = 1. Para entender mejor lo que hace el siguiente ejemplo puede ser conveniente eliminar los “;” la primera vez que lo ejecutes; de este modo el programa te irá devolviendo los valores seg´ un los vaya calculando. Ejemplo 3. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “factorial1.m”. function nfactorial=factorial1(n) % factorial1(n) % Devuelve n*(n-1)*...*3*2*1



107

if n==0 nfactorial=1; else nfactorial=n*factorial1(n-1); end end Ahora, volvemos a la consola y probamos la función.  factorial1(0)  factorial1(5) Las funciones recursivas son muy frecuentes en programación porque permiten resolver programas muy complejos, pero su estudio y uso exceden los objetivos de este curso. Veamos ahora las estructuras de control que permiten la repetición. Existen dos tipos principales: Repetición de una serie de instrucciones para cada uno de los elementos de un vector. A estos bucles se los denomina bucles for.

Repetición de una serie de instrucciones mientras se cumpla una determinada condici´ on. A estos bucles se los denomina bucles while.

Un bucle for tiene la siguiente sintaxis: for k=v instrucciones a ejecutar end donde v es un vector y k una variable. Las instrucciones se ejecutarán una vez por cada elemento del vector v, comenzando por el primero y avanzando secuencialmente. En cada ciclo, k tendrá el valor del elemento del vector correspondiente.


108


Ejemplo 4. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “sumatorio.m”. function suma=sumatorio(v) % sumatorio(v) % Devuelve la suma de los elementos del vector v suma=0; for k=v suma=suma+k; end end Ahora, volvemos a la consola y probamos la función.  sumatorio([1,3,4,5]) Es muy habitual en un bucle for que el ´ındice recorra las posiciones del vector, en lugar de los elementos, ya que a partir de una posición podemos obtener el valor del elemento en dicha posición, pero no al revés. Ejemplo 5. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “maximo.m”. function imaximo=maximo(v) % maximo(v) % Devuelve la posici´ on del elemento mayor de v imaximo=1; for k=2:length(v) if v(k)>v(imaximo) imaximo=k; end end end Ahora, volvemos a la consola y probamos la función.  v=[1,6,4,5]  maximo(v)  v(ans)



109

Como regla general, cada vez que en MATLAB creas que necesitas usar un bucle for, piensa si puedes definir la operación vectorialmente; será más corto, m´ as rápido y m´ as claro. Por ejemplo, para calcular el n´ umero de elementos positivos de un vector v, se puede usar sum(v>0). Veamos ahora los bucles while. Estos bucles se utilizan cuando no se sabe cuántas veces hay que repetir las instrucciones, sólo que hay que hacerlo mientras se cumpla una condici´ on (o hasta que no se cumpla una condición). Ejemplo 6. Ve al editor, teclea las siguientes l´ıneas y guárdalo como “numeroveces.m”. La funci´ on rand genera una matriz (o un vector en este caso) de n´ umeros aleatorios entre 0 y 1. function nveces=numeroveces(n,m) % numeroveces(n,m) % Devuelve cu´ antas veces ha tenido que obtener n n´ umeros % aleatorios entre 0 y 1 para que su suma sea mayor que % m. if n<=m nveces=0; else v=rand(n,1); nveces=1; while sum(v)
4.4.

Entrada y salida

En esta secci´ on veremos los comandos para que un “script” o función de


110


MATLAB lea o devuelva datos, por teclado/pantalla. Nótese que hasta ahora tenemos un modo de pasar o recibir datos de una función, y es a través de las variables de entrada y salida. Se recomienda que cuando se haga lectura/escritura de datos a través del teclado/pantalla o de un archivo, se utilicen o bien scripts o bien funciones espec´ıficas que sólo hagan la lectura/escritura de datos y la comprobaci´ on de que son correctos. Para pedir que el usuario entre datos por el teclado, usaremos el comando input.

Ejemplo 1. Leemos una distancia (positiva): function d=leedistancia() % d=leedistancia() % Pide al usuario una distancia hasta que sea positiva d=-1; while ˜(d>0) d=input(’Introduce la distancia: ’); end Volvemos a la consola y probamos la función  distancia = leedistancia

En caso de que no se introduzca nada, devuelve una matriz vac´ıa. Se puede comprobar si una variable contiene una matriz vac´ıa usando isempty.

Ejemplo 2. Para indicar que lo que se lee es una cadena de texto, se utiliza la opci´ on ’s’. function v=leevector() % leevector() % Devuelve un vector introducido por el usuario, elemento % a elemento.



111

respuesta=’Si’; v=[]; while respuesta==’Si’ respuesta = input(’Mas datos?\ n [Si]: pulsa enter,\ n [No]: pulsa cualquier tecla y enter\ n ’,’s’); if isempty(respuesta) respuesta = ’Si’; end if respuesta==’Si’ nuevoDato=input(’Nuevo valor:’); v=[v nuevoDato]; end end end Volvemos a la consola y probamos la función  vector = leevector

Para mostrar el uso de la funci´ on display, vamos a crear un peque˜ no script que pida un n´ umero y muestre por pantalla si es primo o no. Copia las siguientes l´ıneas en el editor que estés usando y guarda el archivo como esprimo.m Ejemplo 3. n=input(’Introduce un numero: ’); if isprime(n) disp(’El numero es primo’); else disp(’El numero no es primo’); end El comando disp eval´ ua lo que haya entre los paréntesis y muestra el resultado en pantalla. Es decir, no puede mostrar texto y n´ umeros directamente. Si queremos mostrar texto y n´ umeros en una misma l´ınea, debemos, por ejemplo, transformar los n´ umeros en texto con mat2str y unir ahora los textos con strcat. Vamos a modificar el script esprimo.m para que en la salida repita el n´ ume-


112


ro introducido. Editalo para que ponga lo siguiente: Ejemplo 4. n=input(’Introduce un numero: ’); nstr=mat2str(n); % Convertimos el n´ umero en texto if isprime(n) disp(strcat(’El numero: ’,nstr,’ es primo’)); else disp(strcat(’El numero: ’,nstr,’ no es primo’)); end

4.5.

Manejo de funciones

Las funciones (an´ onimas o no) se pueden utilizar como argumentos para otra función. Ejemplo 1. La siguiente funci´ on, recibe como argumentos una función f y un punto x y devuelve la integral de f entre 0 y x. Sálvala como “integral.m”. function y = integral(f,x) % integral(f,x) % Devuelve la integral de la funci´ on f entre 0 y x % usando cuadratura adaptativa. y=quad(f,0,x); end Volvemos a la consola para probarla  f = @(x)  integral(f,3) En el ejempo anterior hemos usado la función predefinida de MATLAB, quad, que calcula la integral de una funci´ on mediante cuadratura adaptativa. Otro ejemplo de funci´ on de MATLAB que admite funciones como parámetros es ezplot (en general, todas las de representación gráfica). También las funciones de optimizaci´ on o de resoluci´ on de ecuaciones diferenciales admiten funciones como par´ ametros.



Ejemplo 2. Para pasar una funci´ on anónima, simplemente la usamos como argumento.  f = @(x) x;  integral(f,1)

Ejemplo 3. Para pasar una funci´ on no anónima (un archivo “.m”), colocaremos delante del nombre el s´ımbolo @.  integral(@abs1, 1)


113

114

4.6.


Ejercicios

Antes de comenzar crea un directorio y a˜ nádelo al path de MATLAB. Guarda todas las funciones y ejemplos en ese directorio. Ejercicio 1. Crea un “script” que represente la función m´ın(tan2 (x), cos(x)) entre 0 y 2π. Ejercicio 2. Crea una funci´ on que reciba el radio de un c´ırculo r y devuelva 2 el a´rea (πr ). Ejercicio 3. Crea una funci´ on que reciba el radio r de un c´ırculo y devuelva el a´rea y el per´ımetro (2πr). Ejercicio 4. Crea un script que genere un vector aleatorio con 10 n´ umeros entre 0 y 1, v. Usando arrayf un y la función creada en el ejercicio 2, obtén un vector con las ´ areas de los c´ırculos de radio cada uno de los elementos de v. Ejercicio 5. Crea una funci´ on que reciba dos valores a y b y devuelva a/b si b = 0 y a2 + b2 > 1 y devuelva 0 en caso contrario. Ejercicio 6. Crea una funci´ on que reciba un vector y devuelva cuántos elementos positivos tiene. Inténtalo usando un bucle for y si tienes tiempo, sin usarlo. Ejercicio 7. Crea una funci´ on que reciba un valor x e y, repita x = x/2 hasta que el valor sea menor que y y devuelva ese valor. Prueba la función con valores x e y positivos. Ejercicio 8. Crea un “script” que pida al usuario el radio de una circunferencia, llame a la funci´ on creada en el ejercicio 2 y muestre por pantalla el área de la circunferencia. Ejercicio 9. Crea una funci´ on que reciba como parámetros una función f , y dos valores x y h y devuelva la estimaci´ on n´ umerica de la derivada de f en x dada por f (x + h) − f (x − h) f  (x)  2h



115

Ejercicio 10. Crea una funci´ on que reciba como parámetros una función f , dos valores a y b y un n´ umero natural n y devuelva la integral de f entre a y b por el método de los trapecios con n tramos.   b f (a) + 2 n−1 k=1 f (a + kh) + f (b) , h = (b − a)/n. f (x) dx  h 2 a Inténtalo con un bucle for y sin bucle.


116

4.7.


Soluciones

Ejercicio 1. ezplot(’min(tan(x)^2,cos(x))’,[0,2*pi]) Otra forma: x = linspace(0,2*pi); y = min(tan(x).^2,cos(x)); plot(x,y) Ejercicio 2. function A = ejer2(r) A = pi*r^2; end Ejercicio 3. function [A,P] = ejer3(r) A = pi*r^2; P = 2*pi*r; end Ejercicio 4. v = rand(1,10); arrayfun(@ejer2,v) Ejercicio 5. function c = ejer5(a,b) if (a^2+b^2) >1 && b >0 c = a/b; else c = 0; end end Ejercicio 6. Usando un bucle for: function w = ejer6(v) n = length(v);


´ A LA PROGRAMACION ´ CAPÍTULO 4. INTRODUCCION w = 0; for i = 1:n if v(i)>0 w = w+1; end end end Otra forma: function w = ejer6b(v) n = length(v); w = 0; for i = 1:n w = w+(v(i)>0); end end Sin usar bucle for: function w = ejer6c(v) w=sum(v>0); end Ejercicio 7. function y = ejer7(x,y) while (x>=y) x = x/2; end y = x; end Ejercicio 8. r = input(’Introduce el radio: ’); ejer2(r); Ejercicio 9. function D = ejer9(f,x,h) D = (f(x+h)-f(x-h))/2*h; end


117

118


Ejercicio 10. function A = ejer10(f,a,b,n) h = (b-a)/n; k = 1:n-1; A = (f(a) + 2*sum(f(a+k*h)) + f(b))*h/2; end


CAPÍTULO 5. AJUSTE DE CURVAS A DATOS

119

Cap´ıtulo 5

Ajuste de curvas a datos Mostramos en un par de ejemplos c´ omo se utiliza MATLAB para calcular el polinomio de un grado determinado que mejor se ajusta (en el sentido de m´ınimos cuadrados) a un conjunto de pares de datos.

5.1.

La funci´ on polyfit

Dado que solamente vamos a mostrar un par de ejemplos del uso de la función polyfit, es recomendable echar un vistazo a la ayuda de MATLAB si se desea conocer esta funci´ on con mayor profundidad.  help polyfit

En nuestro primer ejemplo calcularemos y dibujaremos la recta que mejor se ajusta a un conjunto de pares de datos. Ejemplo 1. En un experimento realizado entre 1987 y 1989, se determinó el ´ındice de grosor (TI, por sus siglas en inglés, thickness index) de la concha de mejillones recogidos en distintos lugares de la costa de Maine (EE.UU.), as´ı como la concentraci´ on de tributilesta˜ no (TBT) en sus tejidos celulares. El ´ındice de espesor TI se define como el cociente de 100 veces el grosor de la concha divido por su longitud y la concentración de TBT se midió en micro-


120


gramos (mg) por gramo (g) de tejido. El resumen de los resultados obtenidos vienen dados en la tabla siguiente: TBT(mg/g) TI

0,04 0,11 0,50 0,60 0,73 0,77 1,75 4,31 0,40 0,42 0,43 0,43 0,46 0,46 0,48 0,53

TBT(mg/g) TI

1,66 3,56 0,27 0,45 0,74 1,75 0,04 0,04 0,50 0,56 0,48 0,45 0,44 0,48 0,42 0,39

Estos datos los puedes descargar del fichero TI TBT.mat que está en http: //www.editorialtebar.com. Desc´ argalo donde MATLAB pueda encontrarlo y escribe:  load TI TBT Recuerda que una vez cargado en MATLAB, puedes ver los nombres de las variables usando who. Pintemos ahora los pares de puntos (TI, TBT) con c´ırculos azules (o como t´ u prefieras) y de paso pongamos t´ıtulo a la gráfica y a los ejes.     

plot(TI,TBT,’bo’) xlabel(’Indice de espesor (TI)’) ylabel(’Concentracion de TBT’) title(’Contaminacion marina’) hold on

Observa que hemos activado la opci´ on hold para poder seguir dibujando sobre la misma gr´ afica. A continuaci´ on, calculemos la recta que mejor se ajusta a nuestros datos:  p = polyfit(TI,TBT,1)



121

El par de valores obtenidos son la pendiente y el término independiente, respectivamente, de la ecuaci´ on de una recta en la forma y = bx + a. El u ´ltimo argumento (el uno) indica precisamente que queremos una recta (grado 1). Si, por ejemplo, quisiéramos hallar la par´ abola que mejor se ajusta a los datos, bastar´ıa escribir polyfit(TI,TBT,2) en este caso obtendremos un polinomio de MATLAB de grado 2, es decir, un vector con tres coordenadas a, b y c de tal forma que nuestra par´ abola ser´ıa y = ax2 + bx + c. Finalmente, dibujamos en nuestra gr´ afica la recta obtenida:  x = linspace(0.35,0.65);  y = polyval(p,x);  plot(x,y,’r’); Si quieres, puedes calcular el polinomio de grado dos o tres que mejor se ajusta a nuestros datos y pintarlos con otro color en la figura activa. Veamos, en el siguiente ejemplo, c´ omo se puede calcular la exponencial que mejor se ajusta a un conjunto de puntos dados. Antes de comenzar cerremos todas las figuras y borremos todas las asignaciones a variables.  close all  clear all

Ejemplo 2. Con el fin de determinar la vida media de una sustancia radiactiva mediante la observaci´ on directa de los datos, la sustancia debe desintegrarse con la suficiente rapidez para que se pueda medir realmente su cambio. No es sorprendente, pues, que las primeras sustancias radiactivas estudiadas fueron aquellas de desintegraci´ on r´ apida. Los datos de la siguiente tabla provienen de un art´ıculo publicado por Meyer y von Schweidler en 19051 . Las entradas de 1

St. Meyer and E. v. Schweidler, Untersuchungen u ¨ber radioaktive Substanzen. V. Mit¨ ¨ teilung. Uber Radioblei und Radium-Restaktivit¨ aten, Sitzungsberichte der Osterreichischen Akademie der Wissenschaften, mathem.-naturw. Kl., Abt. IIa, 114 (1905), 1195–1220.


122


la segunda columna no corresponden a cantidad de materia sino a mediciones de actividad relativa. Dado que, seg´ un la ley de desintegración radiactiva, la tasa de emisión radiactiva es proporcional a la cantidad presente de materia, se pueden medir cantidades de materia (relativas a la cantidades iniciales) observando en su lugar las tasas de emisi´ on (relativas a la tasas de partida). Tiempo (d´ıas) Nivel de act. rel.

0,2 2,2 4,0 5,0 6,0 8,0 11,0 35,0 25,0 22,1 17,9 16,8 13,7 12,4

Tiempo (d´ıas) Nivel de act. rel.

12,0 15,0 18,4 26,0 33,0 39,0 45,0 10,3 7,5 4,9 4,0 2,4 1,4 1,1

Estos datos los puedes descargar del fichero radiactivo.mat que está en http://www.editorialtebar.com. Desc´ argalo donde MATLAB pueda encontrarlo y escribe:  load radiactivo Recuerda que una vez cargado en MATLAB, puedes ver los nombres de las variables usando who. Pintemos ahora los pares de puntos (tiempo, nivel) con cuadrados magenta (o como t´ u prefieras) y también podemos aprovechar para poner nombres y t´ıtulos a nuestra gr´ afica.     

figure(1) plot(tiempo,nivel,’ms’) xlabel(’Tiempo (en dias)’) ylabel(’Nivel de actividad relativa’) title(’Datos de Meyer-v.Schweidler’)

Veamos qué ocurre cuando dibujamos (con cruces rojas, por ejemplo) nuestros pares de puntos con escala logar´ıtmica en el eje OY .



123

 figure(2)  semilogy(tiempo,nivel,’r+’)  hold on  xlabel(’Tiempo (en d´ ıas)’)  ylabel(’Log (nivel de actividad relativa)’)  title(’Datos de Meyer-v.Schweidler (en escala ... semilogaritmica)’)

Observamos que ahora los puntos est´ an prácticamente alienados. Calculemos la recta que mejor se ajusta a estos puntos.  p = polyfit(tiempo,log(nivel),1) Nótese que hemos tomado logaritmo en las ordenadas.  x = linspace(0,45);  y = exp(polyval(p,x));  semilogy(x,y,’b’) Obsérvese con atenci´ on el uso de la exponencial en la segunda l´ınea de la anterior lista de comandos. Dado que nuestro segundo dibujo ha sido realizado en escala semilogar´ıtmica, hemos tenido que hacer uso de log y exp para ajustar los datos al dibujo. Veamos finalmente el resultado de nuestros cálculos en la primera figura.  figure(1); hold on  plot(x,y,’b’)


CAPÍTULO 6. AUTOVALORES Y AUTOVECTORES

125

Cap´ıtulo 6

Autovalores y autovectores Resumimos las ordenes elementales de MATLAB para calcular los autovalores y autovectores de una matriz cuadrada (diagonalizable) y las ilustramos con algunos ejemplos elementales. Incluimos asimismo un ejemplo de aplicaci´ on del cálculo de autovalores y autovectores.

6.1.

La funci´ on eig

Comencemos recordando qué se entiende por matriz diagonalizable: diremos que una matriz cuadrada A ∈ Cn×n es diagonalizable si existe una matriz invertible P ∈ Cn×n tal que P −1 AP = D := diag(d1 , . . . , dn )

(6.1)

donde diag(d1 , . . . , dn ) denota la matriz diagonal cuya entrada (i, i)-ésima es di ∈ C, i = 1, . . . , n. Los elementos de la diagonal de D se llaman autovalores de A y son las ra´ıces del polinomio det(A − xIn ), donde In es la matriz identidad de orden n, este polinomio se llama polinomio caracter´ıstico de A. Para cada autovalor λ de A, el sistema de ecuaciones lineales (A−λIn )x = 0 es compatible indeterminado (es decir, tiene infinitas soluciones), sus soluciones no nulas se llaman autovectores de A asociados a λ. Se puede comprobar f´ acilmente que si pi es la columna i−ésima de la matriz P de la ecuación (6.1), entonces pi es un autovector de A asociado a di .


126


La función eig de MATLAB calcula una matriz P y una matriz diagonal D tales que P −1 AP = D, cuando A es una matriz diagonalizable.  A = [4,-1;6,-1]  [P,D] = eig(A) Por supuesto podemos cambiar los nombres de las variables que contendrán la matriz de paso y la matriz diagonal:  [Q,E] = eig(A) Por defecto, la funci´ on eig normaliza las columnas de P. En MATLAB (no en Octave), podemos evitarlo a˜ nadiendo la opción ’nobalance’. MATLAB>> [P,D] = eig(A,’nobalance’)

Es importante tener presente que por razones de eficiencia computacional, MATLAB no resuelve el problema de calcular P invertible y D diagonal tales que P −1 AP = D, sino que resuelve el siguiente problema: calcular P y D tales que AP = P D. A simple vista podr´ıa parece equivalente al nuestro, pero si A no es diagonalizable la cosa no funcionar´ a como esperamos.  A = [0,1;0,0]  [P,D] = eig(A)  inv(P) % P no es invertible. Si queremos calcular el polinomio caracter´ıstico de una matriz A (no necesariamente diagonalizable), podemos usar la función poly. As´ı, los autovalores de A se pueden calcular usando la funci´ on roots aplicada al polinomio caracter´ıstico, aunque también podemos usar la función eig con ninguna o con una sola variable de salida.



    

127

B = [1,-1,0;0,2,-3;0,0,3] p = poly(B) roots(p) % Autovalores de B eig(B) % Sin variable de salida d = eig(B) % Con una variable de salida

La función roots no es muy precisa cuando el polinomio tiene ra´ıces m´ ultiples. Por lo que es recomendable usar la función eig para calcular los autovalores de una matriz. Veamos un divertido ejemplo. Ejemplo 1. Consideremos la matriz diagonal A de orden 7 cuya diagonal está formada por menos unos, es decir, la opuesta de la matriz identidad de orden 7.  A = -eye(7) Es claro que A tiene un s´ olo autovalor, −1, repetido siete veces.  eig(A) Veamos qué ocurre si primero calculamos el polinomio caracter´ıstico de A y luego sus ra´ıces.  p = poly(A)  roots(p) ¡No coincide ni uno!

6.2.

Un ejemplo “cl´ asico”

Consideremos el problema de la movilidad social que involucra la transición entre distintas clases sociales a través de las generaciones sucesivas de una familia. Supongamos que cada individuo es clasificado socialmente, seg´ un su ocupación, como clase alta, media o baja, que etiquetamos como estados 1, 2 y 3, respectivamente, y que la matriz de transición que relaciona la clase de un hijo con la de su padre es   0,45 0,05 0,05   P =  0,45 0,70 0,50  . 0,10 0,25 0,45


128


De tal forma que, por ejemplo, la probabilidad de que un hijo sea de clase alta, media o baja cuando su padre es de clase baja viene dada por la u ´ltima columna de P. As´ı, si p0 es el vector de las proporciones iniciales de hombres de cada clase social, la proporci´ on de hombres en cada clase social tras n generaciones ser´ a: pn = P pn−1 = P (P pn−2 ) = . . . = P n p0 . Un simple an´ alisis de los autovalores y autovectores de P revela que P tiene un u ńico autovector p∞ cuyas coordenadas son positivas y suman 1. Veamos cómo podemos calcularlo:

>> P = [0.45, 0.05, 0.05; 0.45, 0.70, 0.50; 0.10, ... 0.25, 0.45] >> [U,D] = eig(P); [a,i] = max(eig(P)); >> u = -U(:,i); pinf = u/sum(u)

El vector p∞ predice la distribuci´ on por edades tras una cantidad considerable de generaciones. Nuestra predicci´ on es, por tanto, que en el futuro la población masculina consistir´ a en un 8.3 % de clase alta, un 62 % de clase media y un 29.7 % de clase baja. Comprobemos experimentalmente que la proporción entre las entradas de p∞ rige la distribuci´ on por clases cuando el tiempo es suficientemente grande. Comencemos fijando un vector que recoja la proporción inicial de población masculina de cada clase, por ejemplo, 

 0   p0 =  0,25  , 0,75 y veamos en una gr´ afica que efectivamente la proporción de población en cada clase se estabiliza con el tiempo.



>> >> >> >> >> >> >> >> >>

129

p0 = [0;0.25;0.75]; g = 10; % Numero de generaciones. X = zeros(3,g); X(:,1)=p0; for t=2:g X(:,t) = P*X(:,t-1); end plot(X’) legend(’Clase alta’,’Clase media’,’Clase baja’)

Si quieres puedes cambiar el n´ umero de generaciones y observar que la distribución por clases se estabiliza igualmente. Veamos ahora que nuestra predicci´ on de equilibrio se ajusta a la realidad. Observemos, por ejemplo, qué ocurre tras 10 y 100 generaciones: >> >> >> >>

p10 = P^10*p0 err10 = norm(p10 - pinf)*100/norm(p10) P100 = P^100*p0 err100 = norm(p100 - pinf)*100/norm(p100)

Observemos finalmente que el resultado es el mismo para cualquier dato inicial: >> p0 = rand(3,1);p0 = p0/sum(p0) >> p100 = P*p0 >> err100 = norm(p100 - pinf)*100/norm(p100)


CAPÍTULO 7. ANIMACIONES CON MATLAB

131

Cap´ıtulo 7

Animaciones con MATLAB Aqu´ı se explican los comandos b´ asicos para realizar una animación y crear una pel´ıcula con ella. Desafortunadamente, a nivel básico éstas sólo pueden ser realizadas con el programa MATLAB1 . Además, al igual que en el cap´ıtulo sobre Gr´ aficos es importante recordar que al ejecutar cualquier comando de tipo gr´ afico, se obtiene una nueva ventana llamada Figure 1 donde se dibuja el objeto o se realiza el cambio requerido. As´ı, después de ejecutar la primera orden gráfica, organiza el escritorio de forma que puedas ver a la vez, tanto la ventana de comandos de MATLAB como la ventana Figure 1.

7.1.

Creaci´ on de animaciones

La forma más sencilla de producir una animación es utilizando el comando comet introducido en la sesi´ on Gr´ aficos, donde pod´ıamos ver cómo se trazaba una curva en el plano dada en forma paramétrica. Ahora, para producir animaciones m´ as complicadas necesitamos introducir el comando getframe que nos permitir´ a guardar la imagen (o fotograma) que se acaba de representar en la ventana Figure o ventana gráfica activa. En el siguiente ejemplo crearemos una funci´ on que devuelva un vector formado por 100 imágenes. Adem´ as, al ejecutar dicha función obtendremos una animaci´ on. 1

Para realizar animaciones utilizando Octave se recomienda seguir las indicaciones de Karol Krizka dadas en su web http://www.krizka.net/2009/11/06/ creating-animations-with-octave/


132


Es decir, veremos c´ omo la ventana Figure va cambiando como consecuencia de los comandos que se van ejecutando.

Ejemplo 1. Abrimos el editor de programas de MATLAB y escribimos: function M=string % Animaci´ on para ver como vibra un cord´ on. X=0:0.01:1; n=100; % N´ umero de fotogramas. for j=1:n plot(X,sin(j*pi/5)*sin(pi*X)) axis([0,1,-2,2]) M(j)=getframe(gcf); % Guardamos el fotograma j. end end Ahora, guardamos el archivo creado con el nombre string.m en el directorio de trabajo y volvemos a la ventana de comandos de MATLAB.  figure % Activamos una ventana gr´ afica. IMPORTANTE: Ahora nos aseguramos que la ventana Figure no se cubre con ninguna otra ventana ANTES de teclear la orden siguiente. Vemos la animaci´ on y guardamos los fotogramas en el vector M.  M=string;

7.2.

Pel´ıculas o MOVIES

Una vez creada una animaci´ on nos gustará poderla reproducir tantas veces como deseemos. Para ello utilizaremos el comando movie.



133

Ejemplo 1. Reproducimos 4 veces la animación M a 25 fotogramas por segundo.  movie(M,4,25) Como podéis observar la pel´ıcula se reproduce en la ventana Figure. Sin embargo, MATLAB posee una ventana especial con opciones adecuadas para la reproducción de pel´ıculas. El comando movieview reproducirá la pel´ıcula en ésta ventana. Ejemplo 2. Reproducimos la animaci´ on M en la ventana adecuada y mostramos el t´ıtulo de la pel´ıcula.  movieview(M,’Cuerda vibrando’) Finalmente, utilizando la funci´ on movie2avi podremos guardar la pel´ıcula generada a partir de una animaci´ on en un fichero AVI para poder visualizarla en cualquier otro programa de reproducci´ on.2 Ejemplo 3. Creamos la pel´ıcula ’string.avi’ a 25 fotogramas por segundo a partir de la animaci´ on M. WINDOWS movie2avi(M,’string.avi’,’fps’,25, ... ’compression’, ’Indeo3’) LINUX movie2avi(M,’string.avi’,’fps’,25, ... ’compression’, ’None’)

7.3.

Ejemplos de pel´ıculas o MOVIES

Recreando una cama el´ astica

2

WINDOWS indicar´ a que la orden debe ejecutarse solamente si el sistema operativo es WINDOWS y LINUX indicar´ a que la orden debe ejecutarse solamente si el sistema operativo es LINUX.


134

Aprende MATLAB/OCTAVE en 25 horas Abrimos el editor de programas de MATLAB y escribimos:

function C=cama % Animaci´ on para recrear una cama el´ astica. [x,y]=meshgrid(-1:.1:1); n=200; % N´ umero de fotogramas. for j=1:n t=(2*pi/19)*(j-1); z=2*sin(t)*exp(-x.^2-y.^2); surf(x,y,z) axis([-1,1,-1,1,-2,2]) C(j)=getframe(gcf); % Guardamos el fotograma j. end end Ahora, guardamos el archivo creado con el nombre cama.m en el directorio de trabajo y volvemos a la ventana de comandos de MATLAB.

 figure % Activamos una ventana gr´ afica. IMPORTANTE: Ahora nos aseguramos que la ventana Figure no se cubre con ninguna otra ventana ANTES de teclear la orden siguiente. Vemos la animaci´ on y guardamos los fotogramas en el vector C.  C=cama; Finalmente, reproducimos la animaci´ on C. Primero, lo hacemos a 25 fotogramas por segundo, segundo la reproducimos en la ventana adecuada y finalmente creamos la pel´ıcula cama.avi.



 movie(C,1,25)  movieview(C,’Cama El´ astica’) WINDOWS movie2avi(C,’cama.avi’,’fps’,25, ... ’compression’, ’Indeo3’) LINUX movie2avi(C,’cama.avi’,’fps’,25, ... ’compression’, ’None’)

Creando un reloj

Abrimos el editor de programas de MATLAB y escribimos:

function R=reloj % Animaci´ on para ver girar las agujas de un reloj. n=100; % N´ umero de fotogramas. for j=1:n t=0:((2*pi)/1000):(2*pi); plot(cos(t),sin(t)) % Dibujamos el reloj axis square hold on horas=0:12; plot(.9*cos(horas*2*pi/12),.9*sin(horas*2*pi/12),’k*’) hor=pi/2-(j-1)*2*pi/(n-1); % horaria plot([0 .5*cos(hor)],[0 .5*sin(hor)]) min=pi/2-(j-1)*12*2*pi/(n-1); % minutera


135

136


plot([0 .8*cos(min)],[0 .8*sin(min)]) hold off R(j)=getframe(gcf); % Guardamos el fotograma j. end end Ahora, guardamos el archivo creado con el nombre reloj.m en el directorio de trabajo y volvemos a la ventana de comandos de MATLAB.  figure % Activamos una ventana gr´ afica. IMPORTANTE: Ahora nos aseguramos que la ventana Figure no se cubre con ninguna otra ventana ANTES de teclear la orden siguiente.  R=reloj; % Vemos la animaci´ on y guardamos los fotogramas en el vector R. Finalmente, reproducimos la animaci´ on R. Primero, lo hacemos a 5 fotogramas por segundo, segundo la reproducimos en la ventana adecuada y finalmente creamos la pel´ıcula reloj.avi.  movie(R,1,5)  movieview(R,’Reloj’) WINDOWS movie2avi(R,’reloj.avi’,’fps’,5,’compression’, ’Indeo3’) UNIX movie2avi(R,’reloj.avi’,’fps’,5,’compression’, ’None’)


CAPÍTULO 8. ECUACIONES DIFERENCIALES

137

Cap´ıtulo 8

Ecuaciones diferenciales En este cap´ıtulo veremos como resolver numéricamente ecuaciones diferenciales. Por razones exclusivamente metodológicas, nos centraremos en MATLAB, que tiene diferentes funciones para los distintos métodos de resolución de ecuaciones diferenciales, mientr´ as que Octave emplea esencialmente una sóla función.

8.1.

Ecuaciones Diferenciales Ordinarias

Las diversas funciones que tiene MATLAB para calcular de modo numérico las soluciones de una ecuaci´ on diferencial ordinaria son ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb. Vamos a utilizar ode45, esto es, el método Runge-Kutta (4,5), por ser el que primero se deber´ıa probar. En caso de que el método devuelva alg´ un error o sea demasiado lento, se deberá cambiar a otro método de resoluci´ on (la ayuda de MATLAB proporciona una muy buena información sobre qué método utilizar en cada caso). Para pasar una ecuaci´ on diferencial a MATLAB, debemos escribirla como una ecuación diferencial de primer orden, es decir, de la forma x (t) = f (t, x(t)),

t ∈ R, x(t) ∈ Cn , f : R × Cn → Cn .

La función ode45 tiene la siguiente sintaxis: [t,x]=ode45(f,T,x0),


138


donde f es la funci´ on que define el lado derecho de la ecuación, T es un vector cuyos extremos son los instantes inicial y final de tiempo y x0 es el valor de x(t) en el instante inicial. Consideremos la ecuaci´ on log´ıstica: x = x(x − 1). Podemos calcular la soluci´ on que comienza en t = 0, x = 1/2, entre t = 0 y t = 1, con el comando  f=@(t,x) x*(x-1);  [T,X]=ode45(f,[0,1],[0.5]); Y representarla  plot(T,X) Puedes probar pintando ahora una gr´ afica con las soluciones que comienzan en t = 0 y x = −0.5, x = 0, x = 0.5, x = 1. Para trabajar con ecuaciones de grado superior, debemos transformarlas en un sistema equivalente. Veamos un ejemplo con el oscilador de Van der Pol: x (t) − µ(1 − x2 (t))x (t) + x(t) = 0. Mediante el cambio y(t) = x (t), tenemos x (t) = y(t),

  

 y  (t) = µ(1 − x2 )y − x. 

Escribimos el lado derecho como una función que recibe dos variables, siendo la segunda un vector de longitud dos y devuelve un vector columna de longitud dos.  mu=2;  f=@(t,x) [x(2),mu*(1-x(1)^ 2)*x(2)-x(1)]’; Obtenemos la soluci´ on determinada por la condición inicial (x(0), y(0)) = (0, 2):  [t,x]=ode45(f,[0,1000],[0;2]);


CAPÍTULO 8. ECUACIONES DIFERENCIALES

139

Y podemos representar bien la evoluci´ on de las coordenadas respecto del tiempo, bien la gr´ afica de (x(t), y(t)).  plot(t,x(:,1),t,x(:,2))  plot(x(:,1),x(:,2)) Prueba a dibujar varias soluciones o la solución con la misma condición inicial pero distintos valores del par´ ametro µ.

8.2.

Ecuaciones en derivadas parciales

El comando para resolver numéricamente ecuaciones en derivadas parciales es pdepe, que resuelve sistemas de ecuaciones parabólicas o el´ıpticas en una variable espacial x y otra temporal t, de la forma:  ∂  m c(x, t, u, ux )ut = x−m x f (x, t, u, ux ) + s(x, t, u, ux ), ∂x

donde t0 ≤ t ≤ tf , a ≤ x ≤ b y m = 0, 1, 2. Además, hay que proporcionar una condición inicial para x cuando t = t0 : u(t0 , x) = v(x), y una condición en la frontera (x = a ´ o x = b) de la forma:  pa(t, u) + qa(t)f (a, t, u, ux ) = 0    pb(t, u) + qb(t)f (b, t, u, ux ) = 0. 

Con las anteriores definiciones, la solución se calcular´ıa como: sol = pdepe(m,cfs,v,pqab,x,t); donde cfs es una funci´ on que devuelve un vector resultado de evaluar las funciones c,f,s, pqab es una funci´ on que devuelve un vector resultado de evaluar pa,qa,pb,qb y v es la funci´ on dada por las condiciones fronteras (todas estas funciones habr´ıa que crearlas, bien como archivos “.m”, bien como funciones anónimas) y x, t ser´ıan los puntos en los que se van a evaluar las variables x y t. Veámoslo mejor con un ejemplo, la ecuación del calor en una dimensión.


140


Supongamos que tenemos la siguiente ecuación del calor con condiciones frontera: ut = uxx ,

u(t, 0) = 0, u(t, 1) = 0,

u(0, x) =

2x . 1 + x2

Comparando la ecuaci´ on anterior con la forma general que necesitamos para MATLAB, tenemos: c(t, x, u, ux ) = 1, f (t, x, u, ux ) = ux , s(t, x, u, ux ) = 0, m = 0. y para las condiciones frontera: v(x) =

2x , pa(t, u) = u, qa(t) = 0, pb(t, u) = u − 1, qb(t) = 0. 1 + x2

Debemos definir ahora las funciones auxiliares en archivos de función “.m”. Comencemos con las que definen la ecuaci´ on: function [c,f,s]=cfs(x,t,u,ux) c=1; f=ux; s=0; end La condición de frontera para t = t0 se puede definir también como función anónima:  v=@(x) 2*x/(1+x^ 2); Las condiciones frontera en los bordes de x también han de definirse en un archivo “.m”: function [pa,qa,pb,qb]=pqab(xa,ua,xb,ub,t) pa=ua; qa=0; pb=ub-1; qb=0; end Finalmente, definimos la malla de puntos donde vamos a calcular la solución:


CAPÍTULO 8. ECUACIONES DIFERENCIALES  x = linspace(0,1,20);  t = linspace(0,2,10); Ahora podemos calcular la soluci´ on:  u = pdepe(0,@cfs,v,@pqab,x,t); y representarla:    

surf(x,t,u); title(’Superficie solucion.’); xlabel(’Distancia x’); ylabel(’Tiempo t’);


141

BIBLIOGRAFÍA

143

Bibliograf´ıa 1. G. Borrell i Nogueras, Introducci´ on Informal a Matlab y Octave. http://iimyo.forja.rediris.es/ 2. R. Echevarr´ıa L´ıbano, Breves apuntes para comenzar con MATLAB. http://personal.us.es/echevarria/documentos/APUNTESMATLAB.pdf 3. J. Garc´ıa de Jal´ on, J.I. Rodr´ıguez y J. Vidal, Aprenda Matlab como si estuviera en primero. http://mat21.etsii.upm.es/ayudainf/ aprendainf/Matlab70/matlab70primero.pdf 4. P. Howard, Partial Differential Equations in MATLAB 7.0. Spring 2005 5. B. Hunt, R. Lipsman, J. Rosenberg, K. Coombes, J. Osborn y G. Stuck, A guide to MATLAB for beginners and experienced users. Cambridge University Press, Cambridge, 2006 6. R. Johnson, MATLAB Programming Style Guidelines. www.datatool. com/downloads/matlab_style_guidelines.pdf 7. I. Ojeda Mart´ınez de Castilla y J. Gago Vargas, Métodos matem´ aticos para estad´ıstica. Manuales UEx, 58, 2008 8. A. Quarteroni y F. Saleri, C´ alculo cient´ıfico con MATLAB y Octave. Springer-Verlag Italia, 2006 9. G. W. Recktenwald, Numerical Methods with Matlab: Implementations and Applications. Prentice-Hall 2001


144


10. The Connected Curriculum Project: Interactive Learning Materials for Mathematics and Its Applications. Department of Mathematics, Duke University, Durham, USA. http://www.math.duke.edu/education/ ccp/.


MatlabOctave Aprender MatlabOctave en 25 Horas

Recommend Documents