MATLAB para novatos
1.
2
3
4 5
6
7
8
1 Tabla de contenido
Introducción.........................................................................................................3 1.1 MATLAB MATLAB por primera primera vez...........................................................................4 vez ...........................................................................4 1.2 Uso elemental elemental de MATLAB.................................................... MATLAB......................................................................... ..................... 5 1.3 La ayuda en MATLAB ................................................................................7 1.4 Definición Definición de Variables................................................................................ Variables................................................................................ 8 1.5 Exhibición Exhibición de resultados resultados ..............................................................................9 1.6 Guardar el trabajo trabajo ...................................................................................... 10 1.7 Ejercicios...................................................................................................10 Arreglos Arreglos y matrices............................................................................................ matrices............................................................................................ 11 2.1 Construir Construir arreglos de números.................................................................... 12 2.2 Acceder Acceder a los números números en los arreglos arreglos ........................................................ 18 2.3 Operaciones Operaciones con matrices........................................................................... matrices........................................................................... 21 2.4 Operaciones Operaciones con arreglos arreglos ........................................................................... 24 2.5 Ejercicios...................................................................................................27 Programas Programas (scripts) (scripts) y funciones funciones ......................................................................... 29 3.1 Scripts ........................................................................................................29 3.2 Functions Functions ................................................................................................... 30 Trabajar con y mejorar scripts 3.3 scripts ....................................................................31 3.4 Funciones Funciones de funciones funciones .............................................................................. 33 3.5 Ejercicios...................................................................................................36 Control de flujo flujo de programa programa ............................................................................. 37 4.1 Ejemplos....................................................................................................37 4.2 Ejercicios...................................................................................................38 Métodos numéricos............................................................................................ numéricos............................................................................................ 43 5.1 Mínimos Mínimos cuadrados, modelo modelo lineal. lineal. ............................................................ 43 5.2 Promedio Promedio y desviación desviación estándar activo activo ....................................................... 47 5.3 Algoritmo Algoritmo simplex simplex ..................................................................................... 48 5.4 Método indirecto indirecto de solución solución de ecuaciones ecuaciones simultáneas............................ simultáneas............................ 51 5.5 Raíces Raíces de polinomios polinomios (Método de Laguerre) Laguerre) .............................................. ..............................................53 53 5.6 Runge Kutta............................................................................................... 56 5.7 Descomposici Descomposición ón de LU .............................................................................. 59 5.8 Interpolaci Interpolación ón suave (splines) (splines) ...................................................................... 62 Gráficos.............................................................................................................65 6.1 Dos dimensi dimensiones ones ........................................................................................ 65 6.2 Tres dimensiones dimensiones ....................................................................................... 68 6.3 Guardar imágenes imágenes ...................................................................................... 69 6.4 Ejercicios...................................................................................................69 Prácticas Prácticas de buena programación programación ....................................................................... 71 7.1 Favorecer Favorecer vectores vectores en lugar lugar de bucles......................................................... bucles......................................................... 71 7.2 Evitar Evi tar que los arreglos arreglos crezcan mientras mientras se corre el programa. programa. ...................72 ................... 72 7.3 Escribir Escribir MEX-File MEX-File cuando los bucles bucles son inevitables inevitables ................................. 72 7.4 Otras consideracione consideraciones................................................................................. s................................................................................. 73 GUIs..................................................................................................................76 8.1 Ejemplo de script para FFT........................................................................76 8.2 Ejemplo de function para FFT .................................................................... 77
MATLAB para novatos
2
8.3 Ejemplo Ejemplo de GUI para FFT.......................................................................... 77 8.4 Pasarlo Pasarlo a ejecutabl ejecutablee .................................................................................... 80 9 Conexión Conexión con C ................................................................................................. 81 10 Simulink Simulink ........................................................................................................ 85 10.1 Simulink Simulink para resolver ecuaciones ecuaciones algebraicas. algebraicas. .......................................... 85 10.2 Simulink Simulink para resolver ecuaciones ecuaciones diferenciales diferenciales......................................... ......................................... 87 10.3 Simulink Simulink para sistemas sistemas de ecuaciones ecuaciones ........................................................ 88 11 SimMechanics SimMechanics ............................................................................................... 91
MATLAB para novatos
2
8.3 Ejemplo Ejemplo de GUI para FFT.......................................................................... 77 8.4 Pasarlo Pasarlo a ejecutabl ejecutablee .................................................................................... 80 9 Conexión Conexión con C ................................................................................................. 81 10 Simulink Simulink ........................................................................................................ 85 10.1 Simulink Simulink para resolver ecuaciones ecuaciones algebraicas. algebraicas. .......................................... 85 10.2 Simulink Simulink para resolver ecuaciones ecuaciones diferenciales diferenciales......................................... ......................................... 87 10.3 Simulink Simulink para sistemas sistemas de ecuaciones ecuaciones ........................................................ 88 11 SimMechanics SimMechanics ............................................................................................... 91
MATLAB para novatos
3 MATLAB para novatos
Rafael Quintero Torres* 2010 Como resultado de que el alfabetismo en MATLAB en las universidades continúa siendo muy bajo y de que la necesidad de emplear las buenas herramientas de cómputo matemático de manera mas extendida esta presente, es que emprendemos la construcción de este material. El material aquí contenido esta pensado para estudiantes que desconocen MATLAB, pero que les gustaría introducirse en el de manera autodidáctica. También puede ser de utilidad para estudiantes de computación y métodos numéricos ya que los ejemplos y el interés principal es en el control de las acciones de la maquina por medio de la sintaxis de MATLAB. Suponemos que existe un genuino interés o una necesidad en introducirse en esta temática, que se tiene un conocimiento básico de algebra y que se aceptan los retos con singular alegría. 1. Introducción MATLAB es un software para computación en ingeniería y ciencias. Nos ofrece un poderoso lenguaje de programación, posibilidad de realizar y manipular gráficas con facilidad, además de reunir conocimiento acumulado en estas áreas de manera estructurada y útil. út il.
MATLAB se centra en el cómputo y no en las matemáticas, así que las expresiones y la manipulación simbólica no son su fuerte, a excepción de las empleadas con la herramienta symbolic o con el lenguaje MuPAD. Las soluciones son numéricas e inexactas, asociadas a los artificios del redondeo. Esto puede verse como un problema y como una ventaja al mismo tiempo. Comparado con otros lenguajes de programación numéricos como C++ o FORTRAN, MATLAB es mucho más fácil de usar y esta integrado con una biblioteca formada de una gran, gran cantidad de información, módulos, programas, subrutinas, relacionadas con el cómputo numérico. La única comparación donde MATLAB queda en desventaja es en el rubro de velocidad de ejecución, pero puede ser menos dramática con buenas prácticas de programación como veremos aquí. MATLAB no es la herramienta de elección para cómputo de alto desempeño (aunque probablemente si lo sea en el futuro cercano), pero es una plataforma inmejorable para proyectos desde elementales hasta intermedios y en muchos casos para *
UNAM-Querétaro, UNAM-Querétaro, CFATA. rquintero@fata rq
[email protected] .unam.mx
MATLAB para novatos
4
probar lo que en otros lenguajes seria código muy elaborado o que requiere de demasiado tiempo. MATLAB permite obtener una respuesta numérica a una expresión algebraica, de manera similar a una calculadora, también permite obtener una aproximación numéricas a problemas donde la substitución no es posible. Lo que consideramos como la solución enfrenta dos dificultades; que MATLAB entienda una cosa diferente de lo que estamos pensando en virtud de la falta de dominio de la herramienta y que el procedimiento emprendido sea adecuado para alcanzar la solución. Esta escrito, espero, ayude a aminorar el primer problema, el segundo se logra con el dominio de los métodos numéricos y de la experiencia. No se espera que MATLAB pueda resolver un problema que no se pueda resolver con papel, lápiz y tiempo. En general es más exitoso si antes de expresar el problema en MATLAB se ha expresado en papel y lápiz. 1.1
MATLAB por primera vez
Al iniciar MATLAB aparece el escritorio con diversas áreas de trabajo, Multipaneled desktop. Que se puede modificar para ser útil a cada usuario, acción que no es necesaria en un principio. La figura 1.1 ilustra las áreas de trabajo, Command Window (elipse amarilla), donde se le escribe a MATLAB las acciones a realizar, inmediatamente frente al prompt >> (rectángulo azul), y donde aparecen las respuestas numéricas
prompt
Start
Figura 1.1 Multipaneled desktop en MATLAB.
MATLAB es un lenguaje de programación que interpreta las instrucciones línea por línea, a diferencia de otros lenguajes de programación que compilan el programa total. Esta es una afirmación correcta en parte, ya que como veremos existen las funciones que no se interpretan o incluso programas completos en MATLAB que pueden ser ejecutables. Esta parte se aclarara al ver MEX files y GUIs.
MATLAB para novatos
5
(elipse verde), muestra las variables definidas e información al respecto; en un principio esta vacía, al definir variables o realizar operaciones se acumulan y se pueden editar. Doble clic sobre la variable, inicia un editor que permite su visualización y edición, siempre que sea menor a 65 536 renglones o columnas. En esta misma área se puede visualizar la carpeta en la que se esta trabajando al seleccionar la pestaña correspondiente a Current Directory, MATLAB solo ve el contenido de la carpeta que esta seleccionada. Para cambiar el directorio emplear el acceso rápido (circulo azul claro, arriba al centro). Workspace
(elipse rojo), lista las acciones solicitadas en command window, de la más reciente a la más antigua y se pueden pasar al Workspace navegando por medio de la flecha en el teclado (arriba y abajo) o haciendo un doble clic en el comando deseado. Command History
(rectángulo magenta), es el puerto de entrada a muchas de las funciones de MATLAB, herramientas de trabajo, ayudas y mas, usualmente debe estar preparado para pasar horas de sano entretenimiento en los ejemplos y guías de uso. Cuando se le ha pedido a Matlab que realice una actividad que toma tiempo, desaparece el prompt y a la derecha de Start aparece la leyenda Busy, lo que indica la necesidad de tener paciencia, si se desea cancelar el proceso se puede hacer presionando al mismo tiempo las teclas ctrl-c lo que pone disponible el prompt y aparece la leyenda Ready a la derecha de Start . Start
En cálculo numérico elemental, un camino trivialmente simplista es usar una calculadora cuando la tarea es muy simple, después Excel si se requiere de gráficas mas elaboradas y finalmente C++ cuando se tiene en mente acciones no desarrollas por un programa particular. El salto de una herramienta a la otra puede ser intimidante e incluso limitante. Si quisiéramos que MATLAB estuviera presente en esta evolución, igualmente puede ser una alternativa a cualquiera de las herramientas anteriores. Una de las limitaciones para iniciar a usar MATLAB como calculadora (inicialmente y después el cielo es el límite), es que no se parece a una. Superando esta diferencia, se puede fácilmente superar las ventajas de Excel como hoja de cálculo y de ahí pasar a usar MATLAB con más provecho. 1.2
Uso elemental de MATLAB
Si se escribe una expresión valida en el área de command window, la respuesta aparece de manera inmediata. >> 5*3 ans = 15 >>
La multiplicación de dos números se realiza y el resultado se guarda en la variable ans. El prompt indica que esta listo para nuevas instrucciones. >> 3^3
MATLAB para novatos
6
ans = 27 >>
La potencia es igualmente evidente, el resultado se escribe en la variable ans, substituyendo cualquier valor o contenido que tuviera anteriormente. >> 1/0 Warning: Divide by zero. ans = Inf >>
Aquí Inf es la manera de representar ∞, que indica un número más grande que cualquiera que puede ser escrito y que en este caso resulto de intentar dividir el número uno entre cero. >> cos(pi/4) ans = 0.7071 >>
Aquí la función trigonométrica está en radianes, pi representa a ; número de veces que el diámetro del círculo se requiere para formar el perímetro del círculo. La misma función trigonométrica en grados donde el círculo tiene 360 grados se escribe cosd. >> exp(i*pi/4) ans = 0.7071 + 0.7071i >>
Aquí i (también pude ser 1i, “uno-i” que evita confusiones con la variable i) representa a 1 ; unidad de los números imaginarios. i puede ser también una variable, únicamente si no ha sido empleado como variable, entonces toma la representación de los números imaginarios. En cuanto a la notación exp(x) es la representación de e x , donde e es el número de Euler que es aproximadamente 2.718281828. >> 0/0 Warning: Divide by zero. ans = NaN >>
MATLAB para novatos
7
Aquí NaN representa a un resultado indefinido, se formó de Not a Number , no es un número y aquí surge de intentar dividir cero entre cero. >> (3-exp(log(3)))/3 ans = 1.4803e-016 >>
Aquí la respuesta algebraica esperada es cero, la aritmética de punto flotante realiza por aproximaciones los cálculos y en su lugar aparece un número muy pequeño, 1.4803 x 1016 (1.4803e-016=1.4803E-016) ; Los cálculos en la computadora producen errores que con la ventaja del punto flotante (notación exponencial) es muy pequeño. eps representa al valor del error máximo aritmético en mi máquina, que es 2.2204 x 10 -16; cualquier numero mas pequeño que eps es aceptable. Dicho de otra manera el siguiente número mayor que uno y distinguible de uno en mi máquina es 1+eps. >> eps ans = 2.2204e-016 >>
Algunas reservadas es logaritmo natural deenxMATLAB log(x) palabras
es logaritmo natural de 10 log(x) x de x es el logaritmo base log10(x) log10(x) logaritmo base 10cuya de x tangente es x en radianes atan(x) eseselelángulo en radianes atan(x) seno hiperbólico decuya sinh(x) es el ángulo x tangente es x seno hiperbólico sinh(x) raízelcuadrada de x de x sqrt(x) es sqrt(x) raíz cuadrada de x 1.3
La ayuda en MATLAB
MATLAB contienen un enorme material, es imposible dominarlo todo y menos aún recordarlo todo, por lo que es importante es aprenda a aprender y recordar de manera eficiente. No es accidente que al iniciar el programa, la primera línea invite a explorar el material de ayuda que contiene: To get started, select MATLAB Help or Demos from the Help menu. >>
El texto en azul es una liga directa a la ayuda en MATLAB, que abre una ventana como la mostrada en la figura 1.2. En la figura 1.2 se ve la ventana de ayuda; en el panel izquierdo esta el navegador, donde se puede explorar por contenido, por índice, hacer búsquedas o explorar las
MATLAB para novatos
8
demostraciones. Aquí se muestra el navegador por contenido (elipse azul), el libro de MATLAB (elipse roja), las funciones por categoría (elipse verde) y las funciones matemáticas elementales. En el panel derecho se ve la lista de funciones exponenciales, doble clic en cualquier nombre en azul para ver la ayuda preparada en cuanto a su definición y uso.
Figura 1.2 Ventana de ayuda de MATLAB.
Muchas veces no se recuerda o conoce el nombre de la función que hace lo que buscamos, la búsqueda alfabética de funciones puede ser útil. Si recordamos el nombre de la función y solo no recordamos el uso, podemos emplear Command Window y teclear help y el comando del que requiere ayuda, por ejemplo help log y aparece en el Command Window la ayuda al respecto, usualmente simplificada. 1.4
Definición de Variables
Se puede definir las variables usadas >> a=exp(1) a = 2.7183 >>
Aquí a tiene el valor de la constante de Euler. >> c=a*b ??? Undefined function or variable 'b'.
MATLAB para novatos
9
>>
Aquí intentamos guardar en la variable c el resultado de multiplicar el número de Euler (a) por el contenido de b, como no habíamos definido previamente el valor de b nos lo indica y no realiza nada, queda en espera una instrucción valida. La multiplicación a*b tiene la notación correcta, (a)(b) no es aceptable como multiplicación, como tampoco lo es ab ni a.b, estos dos últimos casos corresponden a otras variables, ab y la componente b del objeto a. A continuación están algunas convenciones útiles de tomar en cuenta; comentarios que ignora el intérprete de MATLAB, asignar caracteres alfanuméricos a las variables y continuación de comandos en varias líneas. Los caracteres que aparecen en verde no son interpretados por MATLAB y están en verde cuando van después del símbolo de porcentaje %. >> % este es un comentario, no hay ejecución >> x=a^2; % ; al final de la instrucción evita que aparezca resultado en Command Window >> y='letreros' y='letreros'; ;% las comas definen los caracteres alfanuméricos >> z=1+2+3+... 4+5+6 % ... es una continuación de línea
el
z = 21 >>
1.5
Exhibición de resultados
Los resultados numéricos aparecen en forma de 4 decimales, también se pueden exhibir en forma de 15 decimales (long ) o en forma fraccional aproximada (rat ). ). La palabra format regresa a la forma predeterminada. >> P= pi P = 3.1416 >> format long >> P P = 3.141592653589793 >> format rat >> P P = 355/113 >> format >>
MATLAB para novatos 1.6
10
Guardar el trabajo
File-Save workspace as y se guarda con Si se desea guardar las variables se procede por File-Save el nombre deseado. Para guardar la historia de comandos se procede a seleccionarlos en el panel de Command History y con el botón derecho del Mouse se selecciona Create M File, aparece otra ventana con la historia de los comandos y ahí se puede guardar su contenido en un archivo. Las figuras, mas adelante se vera como generarlas, se pueden guardar en formatos muy diversos, *.fig es el formato nativo de MATLAB, es modificable y el mas conveniente ya que guarda los valores que la generaron, *.eps, es postscript y se pueden anexar a documentos Word y la calidad no cambia al cambiar el tamaño, *.jpg es el mas común pero la calidad es muy limitada y disminuye rápidamente al cambiar de tamaño.
1.7
Ejercicios
1. Evaluar las siguientes expresiones en MATLAB. a) A e
ee
1 b) B cos 1 3 c) C log e 7 d) determinar que número es primo {0,1,2,3,4,5,6,7,8,9,10,97,99} Respuestas: a) A=3.8143e+006 b) B=1.9106 [rad] c) C =0.1718 =0.1718 e) [2, 3, 5, 7, 97]
MATLAB para novatos
11
2 Arreglos y matrices En el centro de MATLAB se posicionan las matrices, el acrónimo del nombre se desprende de MATrix LABoratory, y se espera un uso intensivo de matrices y vectores. Un array es un grupo de números que en MATLAB se especifica con índices, cada elemento del grupo se identifica con uno o varios índices, los índices siempre son enteros que empiezan con el número uno. La cantidad de índices necesarios para especificar a un elemento forma las dimensiones del array.
Para el caso de dos dimensiones, la función del índice es la de representar la posición del renglón y la columna. La estructura de cada índice esta definida como (i,j), donde i denota la posición del renglón y j la posición de la columna:
(i, j ) Renglón
Columna
Por ejemplo, la matriz es un array de números de dos dimensiones, dado como: Columna R e n g l ó n
a11 a 21 a31 am1
a12
a13
a22
a23
a32
a33
am 2
am 3
a1n
a3n amn a2n
son los renglones, líneas horizontales y las columns son las columnas, líneas verticales. Por otro lado, los vectores son arreglos numéricos de una sola dimensión, y pueden ser vectores renglón o vectores columna. Conocer la aritmética de matrices es fundamental para aprovechar las ventajas de MATLAB. rows
MATLAB para novatos
12
Vector columna R e n g l ó n
a11 a 21 a 31 a m1
Vector renglón Columna
a11
a12
a13
a1n
MATLAB no hace distinción entre un escalar; una matriz con un renglón y una columna, y un arreglo de números de múltiples índices. En la manera de escribir las expresiones es que se respeta la manera de hacer el algebra de matrices o el algebra de arreglos. Una matriz de dos dimensiones representa los coeficientes de sistema de ecuaciones lineales, un arreglo de dos dimensiones puede representar la manera de agrupar información por conveniencia, MATLAB no sabe si el conjunto de números es una matriz o un arreglo, la aritmética de matrices esta bien definida, la aritmética de arreglos es arbitraria dependiendo de cómo se arreglaron los números y que se quiera hacer con ellos; se debe tener cuidado en lo que le pedimos hacer a MATLAB. 2.1
Construir arreglos de números
Un arreglo de dos dimensiones se define al incluir los números entre paréntesis cuadrados, las columnas se separan por espacio o por comas y los renglones se separan por punto y coma o por nueva línea. Por ejemplo, una matriz de tres renglones y tres columnas está dada como: >> A = [1 2 3; 4 5 6; 7 8 9] A = 1 4 7
2 5 8
3 6 9
Es importante que los renglones contengan el mismo número de columnas de lo contrario Matlab desplegará un mensaje de error de dimensiones. La manera de definir en Matlab un vector columna está dado como: >> B=[1 ;4 ; 7 ] B = 1 4 7
Un vector renglón puede ser definido como:
MATLAB para novatos
13
>> C=[1 2 3] C = 1
2
3
Un escalar o constante puede ser definido como: >> D=5 D = 5
Después de haber definido esas variables, el usuario puede revisar ciertas propiedades directamente en la ventana de Work Space. Sin embargo, existen diferentes funciones en Matlab que le permiten al usuario obtener información del grupo de datos que se está manejando. Algunas veces, esta información es importante para separar un cierto grupo de datos o tener un control en el número de ciclos de operaciones matemáticas recursivas basado en la dimensiones del array. Por ejemplo, el comando size proporciona el número de renglones y columnas que contiene una variable en particular. Estos parámetros pueden ser almacenados con variables alfanuméricas y tenerlos presentes en la memoria de la sesión. La función está dada como: >> [n,m]=size (A) n = 3 m = 3
Donde, n representa el número de renglones y m representa el número de columnas. En este caso particular se trata de una matriz cuadrada de 3x3. Otra función, relacionada con el número de elementos que almacena un array, es la función length. Esta función proporciona la cantidad de elementos más grande que contiene el array. Por ejemplo si se trata de un vector columna, como aquel definido en la variable B, la función esta dada como: >> length(B)
3
El número 3 representa la dimensión mas grande en la variable B y que corresponde a tres renglones. No es fácil darse cuenta si el comando proporciona el número de elementos en renglón o columna. Pero este comando es utilizado en ciclos de iteración que corresponden al máximo número de elementos en un arreglo. Así se limita un ciclo
MATLAB para novatos
14
de operación en un proceso recursivo matemático sin la necesidad de declarar, por separado, el número de veces que un ciclo se tiene que repetir, ya que esta basado únicamente en las dimensiones del array. Matlab tiene la capacidad de poder representar datos en dos y tres dimensiones. La función dims representa la dimensión en la cual fue definido el array. >> ndims(A) 2
El número dos indica que el arreglo de datos esta representado en dos dimensiones. Cuando se trata de un arreglo de tres dimensiones el resultado numérico es tres. Una de las ventajas de Matlab es que no hay necesidad de realizar programas para poder extraer o identificar cierto grupo de datos numéricos en un array, basta usar adecuadamente las funciones definidas en su biblioteca. Esto le permite al usuario identificar y separar propiedades en un array para después poderlos manejar por separado. Por ejemplo, la función find nos permite encontrar datos numéricos con ciertas propiedades matemáticas. Supongamos que por alguna razón estamos interesados en identificar aquellos datos numéricos que son mayores que 7 en el arreglo definido en la matriz A. La función debe ser declarada como: >> [row,col] = find(A>7) row = 3 3
col = 2 3
El resultado desplegará la posición en el arreglo mediante los índices de los renglones con sus respectivas columnas. En este ejemplo, el resultado indica que hay dos números en el array mayores que siete. Uno de ellos se encuentra en la posición (3,2), es decir, en el renglón tres y columna dos, y el otro se encuentra en la posición (3,3), es decir, en el renglón tres y columna tres. También se pueden emplear las relaciones lógicas para identificar elementos en un array.
MATLAB para novatos
15 idéntico a diferente de menor que mayor que menor o igual que mayor o igual que
== ˜= < > <= >=
En resumen podemos describir las anteriores funciones mediante el siguiente recuadro con la función que desempeña cada una de ellas: size da los elementos en cada dimensión length da los elementos en la dimensión mas grande ndims da el numero de dimensiones find encuentra los índices de elementos diferentes de cero. [row,col] = find(A>7) encuentra en A los electos mayores
a7
Matlab le permite al usuario construir nuevas matrices basándose en los resultados o variables declaradas. El requisito es que todo renglón debe contener el mismo número de datos numéricos para hacerlos compatibles. Por ejemplo, la tarea es agrupar los datos almacenados en la variable A y en el variable B. El resultado tiene que ser almacenado en una nueva variable llamada E. La única forma de hacer compatibles ambos array es anexando el vector columna de la variable B; ya sea a la izquierda o la derecha de la matriz A. De aquí que la función está dada como: >> E=[A B] E = 1 4 7
2 5 8
3 6 9
1 4 7
El resultado de la variable E será una matriz que contiene tres renglones y cuatro columnas. Note que la operación requiere de agrupar ambas variables separadas por un espacio entre ellas. Esta misma acción se logra, como se comento al inicio de esta sección, cuando las variables son separadas por una coma. Como regla, el usuario tiene que memorizar que la función de espacio y coma son equivalentes en un array. Ahora, supongamos que deseamos agrupar el vector renglón de la variable C con la matriz de datos A. La única posibilidad de hacer esto es anexar la variable C es en el primer renglón o el último renglón. La función, entonces, esta dada como: >> F=[A;C] F = 1 4
2 5
3 6
MATLAB para novatos 7 1
8 2
16
9 3
El resultado es una matriz almacenada en una variable F con dimensiones de cuatro renglones y tres columnas. Note que la operación esta caracterizada por punto y coma porque se requiere de definir un nuevo renglón. Después de haber entendido como declarar nuevas variables usando matrices, pasemos a redefinir elementos dentro de un array. De esta manera, es posible modificar o sustituir con nuevos datos ciertas posiciones dentro del array. Suponga que en el array F se desea sustituir el número ocho por un cero. La función está dada como: >> F(3,2)=0 F = 1 4 7 1
2 5 0 2
3 6 9 3
Note que en la estructura se declara primero el array a ser modificado, en este caso el array esta representado con la variable F. Después hay que indicar la posición del índice a ser sustituido y que en este caso corresponde al tercer renglón y segunda columna (3,2). Finalmente usando el símbolo igual se relaciona con el nuevo dato a ser almacenado dentro del array y que en este caso particular es un cero. Cuando se desee eliminar todo un renglón en el array, la función está dada como: >> F(:,2)=[] F = 1 4 7 1
3 6 9 3
La estructura de esta operación modificó el array de la variable F. El símbolo de dos puntos (:) en el índice le indica a Matlab que se tomarán todos los renglones de la segunda columna (2) para ser eliminados por []. De aquí que el símbolo de paréntesis cuadrado denota vacío. Esta misma operación puede realizarse si se cambia la notación del índice por (2,:). Por ejemplo si se declara en Comman Windows: >> H=F, H(2,:)=[] H = 1 3 7 9 1 3
La operación indica que se tomarán todas las columnas del segundo renglón para ser eliminadas []. Se puede incluir mas de una instrucción en un renglón, separado por comas
MATLAB para novatos
17
(en caso de querer ver el resultado en el Comman Windows) o punto y coma (en caso de no desearlo). El usuario puede darse cuenta que el lenguaje y la estructura de Matlab esta basado en una estructura lógica y fácil de entender. Por otro lado, existen diferentes operaciones que son útiles cuando se manejan arreglos de matrices. Una operación básica cuando se trabaja con arreglos matriciales consiste en obtener la transpuesta de una matriz, la cual se construye escribiendo los renglones en las columnas. Si los números de la matriz son reales, la transpuesta de A se escribe como A’. Esta operación es muy útil para acceder a los números de la matriz, por ejemplo para construir matrices a partir de otras. Otra manera de concatenar arreglos es la función cat(dim, A, B) donde los arreglos A y B se unen en la dimensión señalada por dim, por supuesto las dimensiones a pegar tienen que ser compatibles. >> V= cat (1,E,F') V = 1 4 7 1 3
2 5 8 4 6
3 6 9 7 9
1 4 7 1 3
7 8 9 7
1 4 7 1
>> H=cat (2,E',F) H = 1 2 3 1
4 5 6 4
3 6 9 3
>>
La construcción de matrices grandes se puede realizar con la ayuda de los comandos mostrados en el siguiente recuadro. Se deja al lector, como ejercicio, leer la ayuda en línea y definir matrices con ciertas dimensiones. Esto le ayudará para resolver los problemas mostrados al final de la sección. eye(m,n) matriz de m x n con unos en la diagonal y ceros en el resto. matriz de m x n con elementos aleatorios entre 0 y 1. rand(m,n) zeros(m,n) matriz de m x n con todos los elementos ceros ones(m,n) matriz de m x n con todos los elementos unos linspace(a,b,n) produce un renglón de n puntos linealmente espaciados repmat(A,m,n) crea una matriz de m x n copias de A
de a hasta b
Otra función relacionada con el operador dos puntos (:), y muy útil en MATLAB, es que puede generar secuencias numéricas como: >> T=1:10
MATLAB para novatos
18
T = 1
2
3
4
5
6
7
8
9
10
>>
En este caso se genera un vector renglón desde 1 hasta 10. El tamaño de paso entre cada digito es por norma 1. Pero si se desea modificar el tamaño de paso la función debe ser dada como: >>TT= -3:-.5:-6 TT = -3.0000
-3.5000
-4.0000
-4.5000
-5.0000
-5.5000
-6.0000
>>
En este caso el vector renglón formado empieza en -3 y termina en -6. El tamaño de paso de -0.5 esta representado entre los operadores dos puntos (:). De aquí que el formato general para definir un vector está dado como: S=número inicial : paso entre números : número final
el resultado produce un vector renglón o vacío si el número final es menor que el número inicial. Esta función es útil cuando se desean declarar nuevos vectores para relacionarlos con otros valores de una función matemática previamente calculada. Al final, ambos vectores son concatenados para graficarlos. La función linspace realiza una función similar y se recomienda al lector tenerla presente. 2.2
Acceder a los números en los arreglos
Matlab permite extraer datos numéricos del arreglo. Por ejemplo si se desea extraer el segundo renglón del arreglo y almacenarlo en una variable G. La estructura a usar está dada como: >> G=E(2,:) G = 4
5
6
4
La estructura indica que primero se debe declarar la nueva variable, que en este caso es G. El símbolo igual representa que la operación será guardad en esta nueva variable y el grupo de datos se tomará del arreglo original representado por la variable E. Después, se tiene que indicar mediante un índice para referirse a todas las columnas del segundo renglón. En este caso, el número dos en el índice usado indica que nos referiremos al segundo renglón y el símbolo dos puntos (:) se refiere a todas las columnas. La estructura para extraer toda una columna es similar a la anterior con la diferencia de que ahora el índice puede ser declarado como (:,2). En este caso se extraerán todos los renglones de la segunda columna y almacenarlo en la misma variable G.
MATLAB para novatos
19
>> G=E(:,2) G = 2 5 8
El acceso a los números de un arreglo puede involucrar ciertos sectores dentro del arreglo. En otras palabras, no requerimos de extraer todo un renglón o una columna. Para ejemplificar este caso, definamos una nueva matriz usando la función rand . Cabe resaltar que cada vez que se ejecuta esta función, Matlab proporcionará diferentes números porque la variable rand esta relacionada con números aleatorios entre 0 y 1. >> A=rand(2,5) A = 0.9501 0.2311
0.6068 0.4860
0.8913 0.7621
0.4565 0.0185
0.8214 0.4447
El arreglo se compone de dos renglones y cinco columnas. Ahora bien, si se desea extraer un valor almacenado en el arreglo, esto puede ser expresado de dos maneras. La primera estructura está dada como: >> ele1=A(1,3) ele1 = 0.8913
La segunda estructura está dada como: >> ele2=A(5) ele2 = 0.8913 >>
En el primer formato A(1,3) representa al primer elemento de la tercera columna. Que también se puede llamar A(5), quinto elemento del listado, la numeración es consecutiva agotando los elementos de la primer columna se sigue con el resto de las columnas. Por otro lado es posible extraer una parte del arreglo de la siguiente manera: >> mat1=A(1:2,[3 5]) mat1 = 0.8913 0.7621
0.8214 0.4447
MATLAB para novatos
20
>>
Los elementos que se seleccionan para formar mat1 son los primeros dos elementos de las columnas 3 y 5. También se puede obtener una matriz nueva a partir de los elementos de otra matriz, por ejemplo, tomar todos los elementos de la columna y formar una nueva matriz con las columnas 2 1 2 4 5 y 5. >> mat2=A(1:end,[2 1 2 4 5 5]) mat2 = 0.6068 0.4860
0.9501 0.2311
0.6068 0.4860
0.4565 0.0185
0.8214 0.4447
0.8214 0.4447
>>
Aquí end indica el ultimo elemento disponible, así que 1:end es una manera compacta de decir que se tiene que considerar desde el primer elemento hasta el ultimo, equivalente a solo escribir dos puntos (:). Por otro lado cuando se requiere agregar un nuevo elemento al array. Se puede hacer de la siguiente manera: >> mat2(3,1)=8 mat2 = 0.6068 0.4860 8.0000
0.9501 0.2311 0
0.6068 0.4860 0
0.4565 0.0185 0
0.8214 0.4447 0
0.8214 0.4447 0
>>
El elemento que se agrega requiere de tener un renglón extra, lo que no es problema en principio, pero puede crear problemas cuando se requiere que el tamaño del arreglo este cambiando de tamaño constantemente, ya que puede hacer el cómputo mucho más lento. >> F=find(mat2>.6) F = 1 3 4 7 13 16
En este caso F representa al índice único que identifica a los elementos de la matriz mat2 que son mayores que 0.6. La siguiente instrucción identifica estos elementos de la matriz. También se podría escribir [a b] = find(mat2>.6) para identificar los dos índices que etiquetan a los números referidos. >> mat3=mat2(F) mat3 =
MATLAB para novatos
21
0.6068 8.0000 0.9501 0.6068 0.8214 0.8214 >>
2.3
Operaciones con matrices
Aquí solo se enfatiza la información presentada en help para las matrices, las operaciones entre elementos de matrices se especifica en la siguiente sección. Dos nuevos vectores columna son generados conteniendo los siguientes números. >> A=[1;2;3], B=[4;5;6] A = 1 2 3
B = 4 5 6 >>
A+B suma A y B. A y B deben tener el mismo tamaño, a menos que una de las matrices sea un escalar, un escalar se puede sumar a una matriz de cualquier tamaño. C ij Aij Bij + Adición de matrices.
>> C=A+B C = 5 7 9 >>
- Sustracción de matrices .
A-B resta B de A. A y B deben tener el mismo tamaño, a menos que una de las matrices sea un escalar, un escalar se puede sumar a una matriz de cualquier tamaño. C ij Aij Bij >> C=A-B C =
MATLAB para novatos
22
-3 -3 -3 >>
* Multiplicación de matrices .
C = A*B es el producto algebraico lineal de las matrices
A y B. Mas precisamente,
n
C ij
Aik Bkj k 1
Para matrices no escalares A y B, el número de columnas de A debe ser igual al número de renglones de B. Un escalar puede multiplicar matrices de cualquier tamaño. >> C=A*B ??? Error using ==> mtimes Inner matrix dimensions must agree.
1 4 2 5 ? 3 6 >> C=A*B' C = 4 8 12
5 10 15
6 12 18
1 4 5 6 24 5 6 8 10 12 3 12 15 18 >> C=A'*B C = 32 >>
4 1 2 35 32 6 / Diagonal o división derecha de matrices .
B/A es aproximadamente lo mismo que B*inv(A). De manera mas precisa es, B/A = (A'\B')'. Esta operación no tienen razón de existir, se tiene por completes pero no representa ninguna operación o concepto fundamental. \ Diagonal invertida o división izquierda de matrices .
Si A es una matriz cuadrada, A\B es aproximadamente lo mismo que inv(A)*B, excepto que se calcula de manera diferente. Si A es una matriz de n x n y B es un vector columna con n componentes, o una
MATLAB para novatos
23
matriz con varias columnas de ese tipo, entonces X = A\B es la solución de la ecuación AX = B calculada por eliminación de Gauss. Una advertencia aparece si A esta mal escalada o es casi singular. Es muy común encontrar un sistema de ecuaciones de la siguiente manera: a11 a12 a13 x1 b1 a a a x b elementos del sistema variables fuerzas 21 22 23 2 2 a31 a 32 a33 x3 b3 En donde la solución o valores para el vector X se puede determinar de la siguiente forma: X= inv(A)*B=A\B. >> A=[1 2 3;1 5 6;1 3 9] A = 1 1 1
2 5 3
3 6 9
>> B=[32;65;73] B = 32 65 73 >> X=A\B X = 4 5 6 >>
^ Potencia de matrices .
X^p es X elevado a la potencia p, si p es un escalar. Si p es un entero, la potencia se calcula por cuadrados repetidos. Si el entero es negativo, X se invierte primero. Para otros valores de p, el cálculo se realiza con la ayuda de los eigenvalores y eigenvectores, tales que [V,D] = eig(X), determina los valores de V y D, donde V es la matriz de columnas de eigenvectores y D es la matriz con diagonal de eigenvalores, entonces X^p = V*D.^p/V. Aquí AX X , X es un eigenvector y es un eigenvalor. F=3*ones(3) F = 3
3
3
MATLAB para novatos 3 3 >> F^2
3 3
3 3
27 27 27
27 27 27
24
ans = 27 27 27 >>
Si x es un escalar y P es una matriz, x^P es x elevado a la potencia matriz P usando eigenvalores y eigenvectores. X^P, donde X y P son ambas matrices se indica por un error. ' Matriz transpuesta . A' es la transpuesta algebraica de A. Para matrices complejas, es la
transpuesta del complejo conjugado. >> G=A+i*A G = 1.0000 + 1.0000i 4.0000 + 4.0000i 7.0000 + 7.0000i
2.0000 + 2.0000i 5.0000 + 5.0000i 8.0000 + 8.0000i
3.0000 + 3.0000i 6.0000 + 6.0000i 9.0000 + 9.0000i
4.0000 - 4.0000i 5.0000 - 5.0000i 6.0000 - 6.0000i
7.0000 - 7.0000i 8.0000 - 8.0000i 9.0000 - 9.0000i
>> G' ans = 1.0000 - 1.0000i 2.0000 - 2.0000i 3.0000 - 3.0000i >>
expm(A). Exponencial de una matriz .
Debemos usar expm(A) cuando deseamos obtener el exponencial de una matriz, o sea, cuando queremos calcular: I + A + A2/2 + A3/6 + · . 2.4
Operaciones con arreglos
Las operaciones entre los elementos se realizan fácilmente incluyendo un punto antes del operador, siguiendo la secuencia empleada en la sección anterior. >> x=[1;2;3], y=[4;5;6] x = 1 2 3 y =
MATLAB para novatos
25
4 5 6 >>
.+ Adición de arreglos .
A+B suma A y B. La suma de los elementos de las matrices es idéntica a la suma de matrices. .- Sustracción de arreglos .
A-B resta B de A. La resta de los elementos de las matrices es idéntica a la resta de matrices. .* Multiplicación de arreglos .
A.*B es el producto de elemento por elemento de las matrices y deben ser del mismo tamaño, a menos que alguna de las dos matrices sea un escalar. >> c=x.*y c = 4 10 18 >> x.*y' ??? Error using ==> times Matrix dimensions must agree. >> x'.*y ??? Error using ==> times Matrix dimensions must agree. >>
./ División derecha de arreglos .
A./B forma la matriz cuyos elementos se forman de A(i,j)/B(i,j). A y B deben tener el mismo tamaño, a menos que alguno sea un escalar A = 1 4 7
2 5 8
3 6 9
>> E=2*ones(3) E = 2 2 2
2 2 2
2 2 2
>> P=A./E P = 0.5000
1.0000
1.5000
MATLAB para novatos 2.0000 3.5000
2.5000 4.0000
26 3.0000 4.5000
>>
.\ División izquierda de arreglos .
A.\B forma la matriz cuyos elementos se forman de B(i,j)/A(i,j). A y B deben tener el mismo tamaño, a menos que alguno sea un escalar. Esta operación no tienen razón de existir, se tiene por completes pero no representa ninguna operación o concepto fundamental. .^ Potencia de arreglos .
A.^B es la matriz con elementos A(i,j) elevado a la potencia B(i,j). A y B deben ser del mismo tamaño a menos que una sea un escalar. >> F=3*ones(3) F = 3 3 3
3 3 3
3 3 3
9 9 9
9 9 9
>> T=F.^2 T = 9 9 9 >>
.' Transpuesta de un arreglo .
A.' es una matriz de la transpuesta de A. Para matrices complejas, esta operación no hace el complejo conjugado como lo hace la transpuesta de una matriz. >> G=A+i*A G = 1.0000 + 1.0000i 4.0000 + 4.0000i 7.0000 + 7.0000i
2.0000 + 2.0000i 5.0000 + 5.0000i 8.0000 + 8.0000i
3.0000 + 3.0000i 6.0000 + 6.0000i 9.0000 + 9.0000i
4.0000 + 4.0000i 5.0000 + 5.0000i 6.0000 + 6.0000i
7.0000 + 7.0000i 8.0000 + 8.0000i 9.0000 + 9.0000i
>> GG=G.' GG = 1.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 3.0000i >>
MATLAB para novatos
27
exp(A). Exponencial de un arreglo .
elemento se obtiene de exp(aij). 2.5
Esta función produce una matriz donde cada
Ejercicios
1. Construir el siguiente arreglo: 8 2 0 0 0 0 2 8 2 0 0 0 0 2 8 2 0 0 0 0 2 8 2 0 0 0 0 2 8 2 A 0 0 0 0 2 8 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 - 1 0 0 0 0 0
0 0 0 0 0 2 8 2 0 0
0 0 0 0 0 0 2 8 2 0
0 0 0 0 0 0 0 2 8 2
- 1 0 0 0 0 0 0 0 2 8
a) Consulte help para emplear únicamente diag en su construcción. b) Consulte help para emplear únicamente toeplitz en su construcción. c) Construir los siguientes arreglos de la manera más económica.
D
1 2 3 0 1 2 0 0 1 0 0 0 0 0 0
4 3 2 1 0
5 11 1 4 2 3 , E 13 1 2 4 15 1
1 2 1 1 1 2 1 3 1 4
1 3 1 2 1 1 1 2 1 3
1 4 1 3 1 2 1 1 1 2
1 5 1 4 1 3 1 2 1 1
.
2. Si A es una matriz de números aleatorios de 8x8. a) Determinar el máximo valor en cada renglón. b) Determinar el mínimo valor en cada columna c) Determinar el máximo número del arreglo d) Cambiar por cero los elementos que son mayores que 0.5 3. Si A es una matriz mágica (magic) de 4x4, construir una instrucción que le permita suma los elementos de la diagonal principal y de la antidiagonal. 4. Si se tienen un arreglo de 52 números, determine alguna manera de hacer un listado aleatorio de ellos. Respuestas: 1. A=diag(ones(10,1)*8,0)+diag(ones(9,1)*2,1)+ …
MATLAB para novatos diag(ones(9,1)*2,-1)+diag(-1,9)+diag(-1,-9) note que eye(10)*8= diag(ones(10,1)*8,0) B=toeplitz([8 2 zeros(7,1)' -1],[8 2 zeros(7,1)' -1]) B=toeplitz([8 2 0 0 0 0 0 0 0 -1]) D=toeplitz(1:5,zeros(5,1))' E=toeplitz(1./(1:5),1./(1:5))
3.
sum(diag(magic(4))) A=(magic(4)), sum(diag(A(4:-1:1,:)))
4. sortrows([rand(52,1),(1:52)'])
28
MATLAB para novatos
29
3 Programas (scripts) y funciones Los M-File son archivos de texto que contienen instrucciones para MATLAB, y se guardan como *.m Existen dos tipos de archivos de este tipo, los scripts o secuencia de instrucciones para MATLAB, y las functions que realizan una secuencia de acciones con argumentos de entrada y de salida. MATLAB tiene integrado un editor de texto que es muy conveniente ya que reconoce las acciones que se solicitan y aporta ayuda visual como colores diferentes para comentarios, palabras reservadas por MATLAB y texto regular, además de permitir identificar errores de sintaxis y ayuda con indentaciones. El editor de archivos se puede llamar por New M-file y que se puede acceder como acceso directo bajo el menú de File, es la hoja en blanco, véase la figura 1.1.
Los archivos se guardan de inicio en la carpeta \work † o en la carpeta \MATLAB, en mi caso la ruta completa es (C:\Users\bigcheese\Documents\MATLAB), se pueden guardar en otra carpeta pero para ser ejecutados por MATLAB deben estar en una ruta reconocida, editpath o el menú de File se emplea para modificar la ruta. Los archivos no son compilados, solo se escribe el nombre frente al prompt o por el editor y las acciones se ejecutan. 3.1
Scripts
Si abro el editor de archivos y escribo: % Este script me permite practicar la manera de guardar instrucciones % y aumentar mi alfabetismo en MATLAB x= rand(1) % genera un número aleatorio. y=sin(x) % calcula el seno del número aleatorio x
y lo guardo con el nombre tarea1, que se guardara como tarea1.m. Entonces podría ejecutarlo desde el editor de archivos presionado el acceso directo Save and run o podría ejecutarlo escribiendo tarea1 frente al prompt. >> help tarea1 Este script me permite practicar la manera de guardar instrucciones y aumentar mi alfabetismo en MATLAB >> tarea1 x = 0.8913 y = 0.7779 >>
†
C:\Program Files\MATLAB\R2006a\work
MATLAB para novatos
30
Al escribir help tarea1 siempre se escriben los comentarios que se agregaron al principio del archivo. La razón principal para escribir scripts es tener una secuencia de instrucciones secuenciales en MATLAB que al ejecutarse facilitan nuestro trabajo, mas útil, conforme mas elaborada es la lista de instrucciones que contiene. Otra importancia de los scripts es que permite revisar y editar la secuencia de instrucciones en cualquier momento y por terceros. Aquí la advertencia de buenas prácticas de escritura y de comentarios nunca se enfatizara lo suficiente. 3.2
Functions
Las funciones aumentan las posibilidades de MATLAB, se construyen de manera similar a los scripts, solo que la primera instrucción es: function [sal1, sal2, sal3]=nombre (ent1, ent2, ent3)
las variables sal1, etc son las salidas de la función y las variables ent1, etc son las entradas de la función. Se pueden emplear tantas variables como se necesiten, incluso ninguna variable, y pueden tener el nombre que les guste. El nombre que emplee en la función debe ser el mismo del nombre empleado para guardar el archivo m. Las líneas que continúan se guardaron en un archivo llamado leycoseno.m. % % % % %
Esta function calcula la ley de los cosenos, requiere los valores de dos lados del triangulo (B, C) y el ángulo entre ellos (a) en radianes La respuesta es la longitud del otro lado del triangulo (A) el formato es: [respuesta]=leycoseno(B,C,a)
function [A]=leycoseno (B, C, a); A=sqrt(B^2+C^2-2*B*C*cos(a));
Como en el caso de los scripts, escribir help y el nombre en el prompt despliega los comentarios iniciales que escribimos. >> help leycoseno Esta function calcula la ley de los cosenos, requiere los valores de dos lados del triangulo (B, C)y el ángulo entre ellos (a), en radianes La respuesta es la longitud del otro lado del triangulo (A) el formato es: [respuesta]=leycoseno(B,C,a)
Para emplear la función se pueden usar cualquier nombre de variables, a continuación la respuesta del valor del lado del triangulo, cuando los otros dos lados tienen un valor de uno y el ángulo entre ellos es 90 grados se guarda en ans. Recordemos que ans lo genera MATLAB cuando no se asigna una variable para la respuesta y solo existe un ans. >> leycoseno(1,1,pi/2) ans =
MATLAB para novatos
31
1.4142 >>
Una propiedad interesante de las funciones es que disponen de su espacio privado para variables, de tal manera que no puede ver las variables del espacio de trabajo de donde se llamo la función, y libera todo el espacio de memoria que empleó cuando termina, además cada función emplea su propio espacio privado. Para el caso de funciones, es posible agregar sub-funciones dentro del mismo archivo, esto no es posible en los scripts. La función principal puede llamar tantas sub-funciones como necesite. Un comentario final en relación a los archivos-m, la mayoría de las funciones en MATLAB son archivos-m, que se pueden leer y se pueden copiar. Lo que representa una oportunidad para aprender como programan los profesionales en MATLAB y copiar algunos de sus secretos. Un tipo de función que se puede usar es la función anónima, es un tipo de función para el que no se requiere de escribir un m-file, la sintaxis es: f = @(arglist) expression. Como ejemplo veamos la construcción de la raíz cuadrática mayor, donde los argumentos son tres a, b, y c. y la raíz mayor es: b b 2 4ac x max 2a >> f=@(a,b,c) (-b+sqrt(b^2-4*a*c))/2/a f = @(a,b,c) (-b+sqrt(b^2-4*a*c))/2/a >> x=f(1,-(1e8+1e-8),1) x = 100000000
3.3
Trabajar con y mejorar scripts
El editor de MATLAB tiene accesos directos para correr el script , en la figura 3.1 esta dentro del círculo rojo. En el circulo azul esta el setpoint , que permite marcar a la derecha del número de línea puntos de referencia que se pueden usar para correr el script de punto a punto. Esto es ventajoso para probar solo una sección del script , estos marcadores se pueden eliminar con el accedo directo que esta a la derecha. Al momento de detener el
MATLAB para novatos
32
programa se pueden revisar las variables en el Workspace para juzgar la fuente de error, y se puede continuar después de revisar. El acceso directo que esta indicado por el círculo amarillo se activa al correr el script y puede emplearse para avanzar un paso a la vez la ejecución del programa.
Figura 3.1 Apariencia del editor de MATLAB y accesos comunes. En la figura 3.1 y dentro del circulo rosa se encuentra un indicador de sintaxis, el color en la figura es que se tiene una advertencia, en este caso un error que esta mas adelante del script , si se corrige, el color del recuadro cambia a verde. La posición del cursor esta señalada en la esquina inferior derecha, rectángulo gris, esta referencia es conveniente cuando aparecen errores de ejecución e indican en que lugar esta el error. La flecha verde indica las pestañas de los scripts que están abiertos en el editor y cualquiera de ellos puede editarse, solo seleccionándolo. Una capacidad muy provechosa para mejorar el script es profiling, con ella se puede averiguar que tiempo se pasa la computadora en cada función del script y juzgar si es posible mejorar o substituir alguna sección de líneas por otra alternativa, como por ejemplo crear una función. profile on , inicia esta posibilidad, al terminar el script se dice profile report , esto nos da el reporte y cuando no se requiere se deshabilita con profile off .
MATLAB para novatos 3.4
33
Funciones de funciones
MATLAB cuenta con funciones que requieren de definir la función sobre la que operará. La lista incompleta de estas funciones es:
Encontrar las raíces de una función de una variable. Encontrar el valor mínimo en un intervalo de una función de una variable. Encontrar el valor mínimo de un intervalo de una función multidimensional. Encontrar la integral de una función en un intervalo definido. Resolver ecuaciones diferenciales ordinarias con condiciones iniciales.
La manera de definir la función que se usará depende del grado de complejidad de la función; funciones anónimas son las mas simples cuando la función es relativamente simple y un archivo de función cuando la función es muy compleja o cuando se tienen que definir el sistema de ecuaciones diferenciales. Para encontrar las raíces de la función f x 5 3 x 4 13 x 3 91 x 2 196 x 120 , se puede definir la función f por medio de una funciona anónima y evaluarla con fzero como en el siguiente ejemplo, (este polinomio tiene raíces -3, -1, 5, -2-i2 y -2+i2). >> f=@(x) x^5+3*x^4-13*x^3-91*x^2-196*x-120 >> z = fzero(f,6) z = 5
Para determinar el valor de la variable x donde se tiene el valor mínimo de la función f definida anteriormente en el intervalo [0,6] se emplea la función fminbnd como en el siguiente ejemplo. >> [x,fenx] = fminbnd(f, 0, 6) x = 3.5529 fenx = -1.5039e+003
Para determinar el mínimo de una función de dos dimensiones se puede usar la función anónima y fminsearch como en el siguiente ejemplo tomado de la ayuda de MATLAB. z 100( y x 2 ) 2
(1 x) 2
MATLAB para novatos
34
x = -2:.01:2; y = -2:.01:4; [X,Y] = meshgrid(x,y); Z=100*(Y-X.^2).^2+(1-X).^2; A=find(Z>100); Z(A)=100; surf(X,Y,Z)
Figura 3.2 Script para graficar la función limitada en z a 100, (véase sección 5.2), sin importar que en esta escala se percibe una región con posibles mínimos, si se limita el crecimiento de z fácilmente se ve que el mínimo esta en (1,1) como lo indica fminsearch . >> banana = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2; >> [x,fval] = fminsearch(banana,[-1.2, 1]) x = 1.0000
1.0000
fval = 8.1777e-010
Para calcular la integral reprobabilidad de Gauss de 0 a 1 se puede usar una función anónima.
( x ) 2 2 exp f 2 2 >> f = @(x) 2/sqrt(2*pi)*exp(-x.^2/2); >> Q = quad(f,0,1) Q = 0.6827
El caso de ecuaciones diferenciales en general requiere de definir una función con las ecuaciones a resolver, (Simulink permite un mejor control de la simulación y lo recomiendo)
2
d 2 f 2
dt
1
df dt
5 f 3 .05t 3 5 cos(1.5t )
Que podemos reescribir con las funciones y(1)= f y y(2)= f’ que es la manera habitual de trabajar con ecuaciones diferenciales, escribiendo un sistema de ecuaciones diferenciales de primer orden.
MATLAB para novatos y1
f
y 2
df dt
35
dy1 dt
3 .05t 3 5 cos(1.5t ) y 2 5 y1 2 dt
dy 2
Se crea el archivo con nombre derivada1.m function dy = derivada1(t,y) dy=zeros(2,1); dy(1)= y(2); dy(2)= (3+0.05*t^3+5*cos(1.5*t)-y(2)-5*(y(1)))/2;
y se resuelve el sistema de ecuaciones con ode45. La única línea indispensable es la tercera, donde las respuesta se guarda en la columna t, y en dos columnas y una para f y otra para su derivada, si no se agrega el comando options, la simulación procede con los parámetros típicos, [0 10] es el intervalo de la simulación de 0 a 10 segundos y las dos condiciones iniciales son [1.3 8]. clear all options = odeset('RelTol',1e-7,'AbsTol',[1e-7 1e-7]); [t,y]=ode45(@derivada1,[0 10],[1.3 8], options); figure(3) subplot 131 plot(t,y(:,1)); title('1');xlabel('Time [s]'); ylabel('f [a.u.]'); subplot 132 plot(t,y(:,2)); title('2');xlabel('Time [s]'); ylabel('df/dt [a.u.]'); subplot 133 plot(t(1:536),diff(y(:,2))./diff(t)); title('3');xlabel('Time [s]'); ylabel('d(df/dt)/dt [a.u.]');
Figura 3.3 Gráficas de las soluciones de la ecuación diferencial, la sección 5 explica la manera de hacer las graficas, el problema es resuelto una vez más en la sección de Simulink.
MATLAB para novatos
3.5
36
Ejercicios
1. Encontrar las primeras seis raíces de la función de Bessel de orden 1. 2. Encontrar la solución de la ecuación diferencial con condiciones iniciales y(0)=0, dy dz 0
1 en el intervalo de z (0,20] y comparar con besselj(1,z). 2
z
d 2 y 2
dz
z
dy dz
( z 2 1) y 0
Respuesta 1.
2.
3.8317, 7.0156, 10.1735, 13.3237, 16.4706, 19.6159; f=@(x) besselj(1,x); z = fzero(f,4)
Crear el archivo derivada2.m
function dy = derivada2(t,y) dy=zeros(2,1); dy(1)= y(2); dy(2)= (-(t^2-1)*y(1)-t*y(2))/t^2; y resolver con ode45 [t,y]=ode45(@derivada2,[0.0001 20],[0 1])
La comparación entre ambas debe ser igual.
MATLAB para novatos
37
4 Control de flujo de programa La escritura de programas requiere que se pueda repetir alguna instrucción o tomar decisiones respecto a la dirección del programa. MATLAB dispone de estas herramientas. Las herramientas de control de programa se puede dividir en:
Decisiones de control if, else elseif switch, case, otherwise
Bucles for, while, continue, break
Acciones cuando hay error try, catch
Finalizar el programa return 4.1
Ejemplos
Del número x, indicar si es un número negativo, en caso de serlo; si es real no entero dividirlo entre 2 y si es entero obtener su factorial. if x < 0 disp('es un numero negativo' ) elseif rem(x,1) == 0 A = factorial(x) else A = x/2 end
De la variable x, si tiene el valor -1 imprime “uno negativo’, si tiene el valor 0 imprime ‘cero’, si tiene el valor 1, imprime ‘uno positivo’, cualquier otro valor imprime ‘otra cosa’. switch x case -1 disp('uno negativo'); case 0 disp('cero'); case 1 disp('uno positivo'); otherwise disp('otra cosa'); end
MATLAB para novatos
38
La siguiente instrucción se realiza 5 veces, la primera, cuando n=2, produce x(2)=2*1 y sucesivamente toma valores que son el doble del anterior, generando x=[1 2 4 8 16 32]. x=1; for n = 2:6 x(n) = 2 * x(n - 1); end
En el siguiente ejemplo A 1 2 2 3 2 5 2 6 2 7 2 , while debería ejecutarse siempre que x < 10; continue evita que se sume 42 y break termina el bucle cuando x=7. x=x+1 se incluye al principio para evitar que x deje de crecer por continue . Las instrucciones while e if se ejecutan siempre que sean verdaderos lógicos, x==4 es verdadero solo cuando x tienen el valor cuatro. x=1; A=1; while x<10 x=x+1 if x==4 continue end A=A+x.^2 if x==7 break end end
En el siguiente ejemplo se ilustra el posible manejo de errores de interacción con el programa. Aquí se espera que los datos introducidos sean: una constante o una matriz de 2xn. De no ser el caso, pasa el control a catch, que advierte del error. while nos da una nueva oportunidad para que se satisfaga la condición y al logarlo se sale del while. A=[1 2] done = false; while ~done B = input('introduzca datos congruentes: ' ); try C=A*B; done = true; catch disp('Los datos no son congruentes!' ) end end
4.2
Ejercicios
1.- Hacer un script que genere la siguiente expresión: Los primeros 12 números de la serie de Fibonacci son: 0 1 1 2 3 5 8 13 21 34 55 89
2.- Recorra una serie de 200 números aleatorios y graficar las dos columnas, una con los números y la otra con el promedio de los números hasta ese momento.
MATLAB para novatos
39
3.- Hacer el código que permita graficar la siguiente función para x entre -12 y 12: 0 x 10 sin( x) 10 x 0 f ( x) 2 exp( ( x 8) ) 0 x 10 0 x 10 4.- ¡Dibujar un helecho! Este se puede realizar haciendo una transformación repetida de un punto en el plano (x1, x2) por medio de la transformación que toma al punto xi y lo transporta a un nuevo punto x f Axi b , si x es un vector columna de dos elementos, las coordenadas en el plano, A es una matriz de 2x2 y b es un vector columna de dos elementos. a) Tomar un punto arbitrario xi =[0.5;0.5] b) Obtener un numero aleatorio R entre 0 y 1 Si el numero aleatorio esta entre 0 y 0.85 obtener un nuevo punto x f por medio de la transformación x f Axi b , donde A1=[.85 .04;-.04 .85];b1=[0;1.6]; Si el numero aleatorio esta entre 0.85 y 0.92 obtener un nuevo punto x f por medio de la transformación x f Axi b , donde A2=[.20 -.26;.23 .22];b2=[0;1.6]; Si el numero aleatorio esta entre .92 y .99 obtener un nuevo punto x f por medio de la transformación x f Axi b , donde A3=[-.15 .28;.26 .24];b3=[0;.44]; Si el numero aleatorio esta entre .99 y 1 obtener un nuevo punto x f por medio de la transformación x f Axi , donde A4=[0 0;0 .16]; c) graficar las coordenadas del punto sin borrar las coordinas de punto anterior d) regresar al inciso “b” y repetir 10000 veces.
5. Modificar el programa para producir un triangulo de Sierpinski.
MATLAB para novatos a) Iniciar con el punto en el plano s=[0;0] 12 0 b) para cada trasnformacin A 1 0 2
14 12 0 c) se elige b , b , b 3 con la misma probabilidad 0 0 4
Posibles soluciones: 1. N=input('numero de elementos de la serie a mostrar: ' ); x=zeros(N,1); x(1)=0; x(2)=1; for ii=3:N x(ii)=x(ii-2)+x(ii-1); end fprintf('Los primeros %g números de la serie de Fibonacci son:\n' ,N) disp(int2str(x'))
2. a=zeros (200,2); a(1,1)=rand(1); a(1,2)=a(1,1) for ii=2:200 a(ii,1)=rand(1); a(ii,2)=(a(ii-1,2)*(ii-1)+a(ii,1))/ii; end plot(a)
3.
x = linspace(-12,12,1000)'; f=zeros(length (x),1); ind1=find(x>=-10 & x<0);
40
MATLAB para novatos
41
f(ind1)=sin(x(ind1)); ind2=find(x>0 & x<=10); f(ind2)=exp(-(x(ind2)-8).^2); plot(x,f)
3bis.
x=linspace(-12,12,1000); N=length (x) for j=1:N if x(j)>=-10 && x(j)<=0 f(j)=sin(x(j)); elseif x(j)>0 && x(j)<=10 f(j)= exp(-(x(j)-8)^2); else f(j)=0; end end plot (x,f)
4. function []=fern2() tic x=[0.5;0.5]; h=plot (x(1),x(2),'.' (x(1),x(2), '.'); ); set(h,'markersize' set(h,'markersize',1, ,1,'color' 'color', ,'green' 'green', ,'erasemode' 'erasemode', ,'none' 'none'); ); set(gcf,'color' set(gcf, 'color', ,'white' 'white', ,'menubar' 'menubar', ,'none' 'none', ,'numbertitle' 'numbertitle', ,'off' 'off', ,'name' 'name', ,'He lecho') lecho' ) axis([-3 3 0 10]) axis off drawnow p=[.85 .92 .99 1.00]; A1=[.85 .04;-.04 .85];b1=[0;1.6]; A2=[.20 -.26;.23 .22];b2=[0;1.6]; A3=[-.15 .28;.26 .24];b3=[0;.44]; A4=[0 0;0 .16]; for m=1:10000 r=rand; if r
5.
tic a=0; x=[0;0]; plot (x(1),x(2),'.' (x(1),x(2), '.') ) hold on for m=1:10000
MATLAB para novatos p=[.33 .66 1]; A1=[.5 0;0 .5];b1=[0;0]; A2=[.5 0;0 .5];b2=[.5;0]; A3=[.5 0;0 .5];b3=[1/4;sqrt(3)/4]; A4=[.5 0;0 .5]; r=rand; if r
42
MATLAB para novatos
43
5 Métodos numéricos En esta sección se presentarán algunos ejemplos comunes a los cursos de métodos numéricos, se puede saltar esta sección en la primera lectura o se pueden ver los ejemplos de manera arbitraria, la intención es comentar el código y considerarlo como una manera de ver una posible solución a un problema simple. MATLAB tiene soluciones para todos estos problemas de manera mucho mas eficiente y elegante, aquí solo se intenta emplear a MATLAB como el medio ejecutor de las soluciones, controlar a MATLAB a que haga lo que deseamos como medio de aumentar nuestro control de las instrucciones básicas de MATLAB. 5.1
Mínimos cuadrados, modelo lineal.
Tenemos una colección de puntos P( x xi, y yi) que se encuentra dispersos en el plano XY. Podemos decir que cada punto representa el resultado de un experimento, encuesta, censo, etc. Queremos encontrar una representación lineal que mejor se aproxime a nuestros puntos obtenidos. Para eso queremos que la distancia D entre el modelo y el punto experimental sea muy pequeña. Esta distancia la podemos definir como: D yi y *
Donde yi es el valor real y y* es el valor modelado.
Figura 5.1 Error D a minimizar en este método. El modelo lineal para y* la podemos definir como una función del tipo: y* mxi b
MATLAB para novatos
44 D yi
Entonces:
mxi b
Pero queremos que esta distancia D sea pesada de igual manera sin importar si es positiva o negativa así que queremos que esto se cumpla para cada uno de los puntos por lo que debemos agregar la suma de todos los puntos P( xi ,yi). R
D 2
n
( yi mxi b) 2 i 1
Ahora para obtener los parámetros del modelo m y b que mejor se ajustan a los puntos, queremos que el error sea minimizado, para lo que obtenemos las primeras derivadas para igualarlas a cero, de esta forma minimizamos la distancia de los puntos a la función f y determinar nuestras constantes m y b para el modelo lineal. n R 2 ( yi mxi b) xi m i 1 n
y i xi
i 1
n
n
i 1
i 1
m xi 2 b xi
n R 2 ( yi mxi b) b i 1 n
y
n
i
i 1
m xi nb i 1
Donde n es el número de puntos que tenemos. A continuación resolvemos el sistema para encontrar m y b a partir de los puntos tabulados.
n n x i 1 i
b n xi yi i 1 ni 1 n xi2 m xi yi i1 i 1 n
Despejando b de una de las ecuaciones y sustituyendo en otra obtenemos lo siguiente igual para el caso de substituir m.
MATLAB para novatos
45 n
n
n
( y x ) x y i
i
i
i 1
m
i 1
n
n
x
i
n
n
i
i 1
n xi i 1
2
i 1
n
2
n
x y x ( y x ) 2
i
b
n
i 1
i
i 1 n
n
x
2
i
i 1
i
i 1
i
i
i 1
n xi i1
2
Así nos queda una función con coeficientes m y b de la forma y*=mxi+b. La calidad del ajuste se puede estimar por medio del coeficiente de correlación: 2 n 1 n n ( yi xi ) xi yi n i 1 i 1 i 1 2 r n 2 1 n 2 n 2 1 n 2 xi xi yi yi i 1 n i 1 i 1 n i 1 A continuación se muestra el algoritmo realizado en MATLAB para este método de mínimos cuadrados. %Ajuste de puntos por el método de mínimos cuadrados, %ajustando puntos a un modelo lineal y=m*x+b %creamos un arreglo de números de -10 al 10 que conste de 200 números %igualmente espaciados. clear all x=linspace(-10,10,200); n=length(x); %creamos nuestro arreglo de números experimentales, formados por una %correlación definida, y=3x+5 %además se agrega una componente aleatoria % se grafican los puntos. y=(3*x+5)+(randn(1,200))*10; plot(x,y,'*'); %Sección de mínimos cuadrados %Ya que para encontrar m y b el denominador es el mismo lo calculamos y lo %guardamos en den. den = (n*sum(x.^2))-(sum(x))^2; %calculamos las constantes m y b de acuerdo con la formula. m = (n*sum(x.*y)-sum(x)*sum(y))/den; b= (sum(x.^2)*sum(y)-sum(x)*sum(x.*y))/den; %calculamos el coeficiente de correlación, 1 es excelente, 0 es malo.
MATLAB para novatos r=(sum(x.*y)-sum(x)*sum(y)/n)^2/(sum(x.^2)-sum(x)^2/n)/(sum(y.^2)sum(y)^2/n); n,r,m,b %graficamos esta función lineal en la misma ventana para comprarla %con los puntos. f=m*x+b; hold on plot(x,f,'red')
46
MATLAB para novatos
5.2
47
Promedio y desviación estándar activo
El objetivo en este caso es de calcular el promedio y la desviación estándar de un conjunto de números con condición de solo usar el valor actual y el anterior de los cálculos. En este caso xn es el número n de la serie, xn es el promedio de los números del 1 al n y n es la desviación estándar de los números del 1 al n. La regla es que para cada cálculo solo se accede el número y los cálculos inmediatos anteriores. 1 n El promedio se define de la siguiente manera: x xi n
i 1
La desviación estándar se define de la siguiente manera:
1
n
( x x ) n
2
i
i 1
Si se calcula para el evento n y el evento n+1 se tiene la probabilidad de conocer el promedio a partir del número recientemente agregado y del promedio conocido. 1 n xn xi n
i 1
1 n1 1 n 1 n xn xn1 xn 1 xi xi xn 1 1 n 1 i 1 n 1 i 1 n Y para la desviación estándar, se puede obtener la nueva desviación estándar únicamente con el valor del número agregado, del promedio anterior y la desviación estándar anterior. n n 1 n 2 2 2 ; de donde se tiene n x x ( ) ( ) x x n i n xi 2 nxn 2 i n n n
i 1
i 1
i 1
1 n 1 1 n 2 2 2 n 1 x x ( ) xi xn1 (n 1)( xn1 ) 2 i n 1 n 1 i 1 n 1 i1 Substituyendo la suma de los cuadrados con el valor de σ n. 1 n n2 n xn2 xn21 (n 1)( xn1 ) 2 n 1 n 1 %Cálculo del promedio y la desviación estándar después de agregar un %nuevo número y empleando solo los valores del promedio y desviación %estándar conocidos. %creamos un arreglo de 100 números consecutivos en x %creamos un arreglo de 100 números aleatorios en xx %creamos dos arreglos de 100 ceros para guardar los promedios %creamos dos arreglos de 100 ceros para guardar las desviaciones % clear all x=1:100; %xx=x; xx=rand(100,1); xpro=zeros(100,1);
MATLAB para novatos
48
xdev=zeros(100,1); xproa=zeros(100,1); xdeva=zeros(100,1); % %se inicializan los promedios y las desviaciones xpro(1)=xx(1); xdev(1)=0; xproa(1)=xx(1); xdeva(1)=0; % %de acuerdo a las ecuaciones del texto, se calcula el promedio y la %desviación por los dos métodos, el que requiere todos los números y %que emplea solo el último. for m=2:100 xpro(m)=sum(xx(1:m))/m; xdev(m)=sqrt(1/m*sum((xx(1:m)-xpro(m)).^2)); xproa(m)=((m-1)*xproa(m-1)+xx(m))/m; xdeva(m)=sqrt(1/(m)*((m-1)*(xdeva(m-1)^2+xproa(m-1)^2)+xx(m)^2(m)*xproa(m)^2)); end % %gráfica simple de las cuatro curvas, deben ser idénticas por pares. plot(x,xpro,x,xdev) hold on plot (x,xproa,x,xdeva,'r')
5.3
Algoritmo simplex
Una compañía manufactura dos productos (A, B) que requieren de pasar por dos departamentos (α, β). Producto A requiere de 2 horas en el departamento α y 3 horas en el departamento β. Producto B requiere de 4 horas en el departamento α y 3 horas en el departamento β. El tiempo disponible para fabricar ambos modelos es de 100 horas a la semana en el departamento α y 90 horas a la semana en el departamento β. La ganancia esperada del producto A es de 5 pesos y 7 pesos para el producto B. Determinar en número de productos A y B a manufacturar que maximice la ganancia suponiendo que todo artículo fabricado será vendido. X es el número de productos A. Y es el número de productos B. F 5 X 7Y max ; % función objetivo, máxima ganancia en pesos. 2 X 4Y 100 ; % restricción, un máximo de 100 horas en departamento α. 3 X 3Y 90 ; % restricción, un máximo de 90 horas en departamento β. X >=0; % restricción, el número de artículos de A es positivo. Y >=0; % restricción, el número de artículos de B es positivo.
Arreglar la función objetivo para obtener el máximo número para F y modificar las desigualdades agregando dos nuevas variables S 1 y S 2.
MATLAB para novatos F 5 X 7Y 0 2 X 4Y S 1 100 3 X 3Y S 2 90
X, Y, S1 y S2 >=0
49 F 1 5 7 0 0 X 0 0 2 4 1 0 Y 100 0 3 3 0 1 S 1 90 S 2
1. Identifica la columna pivote . Encontrar la columna con el valor más negativo en el renglón de la función objetivo. 2. Identifica el renglón pivote . Encontrar el renglón de restricciones que produzca el número menor al dividir el número de la restricción de los diferentes renglones entre el valor de la columna pivote. 3. Hacer el elemento pivote unitario . La columna y renglón pivote forman el elemento pivote; hacer el elemento pivote unitario, dividiendo el renglón pivote entre el valor del elemento pivote. 4. Eliminar los elementos de la columna pivote . En esa columna solo queda el elemento pivote igual a uno y los otros elementos son cero. 5. Si algún elemento de la función objetivo es negativo repetir el paso uno. 6. Al cumplir la condición anterior se tiene la optimización S 1=S 2=0 y solo falta evaluar X , Y y F .
F 1 1.5 0 1.75 0 X 175 0 0.5 1 0.25 0 Y 25 0 1.5 0 0.75 1 S 1 15 S 2
F 1 X 190 1 0 0 1 0 0 1 0.5 1 / 3 Y 20 0 1 0 0.5 2 / 3 S 1 10 S 2
clear all A=[1 -5 -7 0 0 0;0 2 6.5 1 0 100;0 3 3 0 1 90]; m=length(A(:,1)); falta =0; while sum(A(1,:)<0)>0 [val,ind]=min (A(1,:));% columna pivote [val2,ind2]=min (A(:,ind)/val); % valor pivote & renglón pivote A(ind2,:)=A(ind2,:)/A(ind2,ind); % renglón pivote con elemento pivote unitario for k=1:m% columna pivote con un solo uno if k==ind2 continue end A(k,:)=A(k,:)-A(ind2,:)*A(k,ind); end
MATLAB para novatos falta =falta +1; end if falta ==1 ind=5-ind;% columna ind2=5-ind2;% renglón A(ind2,:)=A(ind2,:)/A(ind2,ind); % renglón pivote con elemento pivote unitario for k=1:m% columna pivote con un solo uno. if k==ind2 continue end A(k,:)=A(k,:)-A(ind2,:)*A(k,ind); End
50
MATLAB para novatos
5.4
51
Método indirecto de solución de ecuaciones simultáneas
El producto punto puede tener la notación matricial y escalar siguiente: a b c n n 1 índice igual a b c c ij ai b j , con ij i 1 j 1 0 índice diferente a b c El producto crus puede tener la notación matricial y escalar siguiente: ab c 1 ijk cíclico n n ci ijk a j bk , con ijk 1 a b c ijk cíclico j 1 k 1 0 índice repetido ab c Un sistema de ecuaciones algebraicas simultáneas tiene la ecuación matricial A x b Y también tiene la notación escalar n
A x
ij j
bi para n 1,2,, n
j 1
Que al sacar xi de la suma, se pude usar para calcular a partir de una aproximación anterior x j, una aproximación posterior de xi. n
Aii xi
Aij x j bi , n 1,2,, n j 1 j i
De donde se obtiene la nueva aproximación de xi n 1 xi bi Aij x j , n 1,2,, n Aii j 1 j i 1. 2. 3. 4.
(1)
Conociendo Aij y bi y una aproximación para x j Calcular el valor de xi a partir de la ecuación (1). Calcular la diferencia entre el nuevo vector xi y el viejo vector x j . Si la diferencia es aceptablemente pequeña ya se encontró la solución, de no ser así proponer que x j es xi y repetir desde el paso 2. 5. Traducir este algoritmo en lenguaje de MATLAB y probarlo para obtener la solución xi de la siguiente ecuación.
MATLAB para novatos
52
4 1 1 x1 12 1 4 2 x 1 2 1 2 4 x3 5 function [x,m]=gaussSeidel1(x,maxiter,epsilon) %Soluciona Ax=b por el método de Gauss-Seidel. % %USO: [x,m]=gaussSeidel1(x,maxiter,epsilon) %Ejemplo: [A,m]=gaussSeidel1([0;0;0],100,1e-7) %INPUT: %x=valor inicial arbitrario del vector solución %maxiter=numero máximo de iteraciones (si no se escribe es 500) %epsilon=tolerancia en el error (si no se escribe es 1e-7) %OUTPUT: %x=vector solución %m=número iteraciones para alcanzar la solución if nargin<3; epsilon=1e-7; end if nargin<2; maxiter=500; end for m=1:maxiter xOld=x; x=func(x); dx=norm (x-xOld); if dx
MATLAB para novatos
5.5
53
Raíces de polinomios (Método de Laguerre)
Si se tiene un polinomio de grado n 1 P n ( x ) a1 x n a 2 x n 1 an x an 1 Un método económico para evaluar el polinomio es: 1a P 0 ( x ) a1 ; P i ( x) ai1 xP i 1 ( x); i 1,2, n De la misma manera se puede demostrar que se obtienen las expresiones para las primeras dos derivadas: 2 P 0 ( x) 0 ; P i ( x) P i 1 ( x) xP i 1 ( x); i 1,2, n 3 P ( x ) 0 ; P ( x) 2 P ( x) xP ( x); i 1,2, n 0
i 1
i
i 1
Encontrar la primera raíz y eliminar la posibilidad de repetirla P n ( x) ( x r ) P n1 ( x ) Si expresamos el nuevo polinomio como sigue: P n 1 ( x ) b1 x n 1 b2 x n 2 bn 1 x bn Al igualar los dos polinomios se obtiene b1 a1 ; bn an rbn 1 i 2,3 n
4 5 6
Método de Laguerre Propone un polinomio especial, que de manera interactiva, puede aproximar a cualquier polinomio (algebraicamente es un caso particular, la computadora lo hace un caso general). 7 P n ( x) ( x r )( x q ) n 1 Derivando el polinomio dos veces P n ( x) n 1 1 8 P n ( x) x r x q 2 1 P n ( x ) P n ( x) n 1 9 2 ( ( P n ( x) P n ( x) x r ) x q) 2 De donde se puede definir P n ( x ) P n ( x) G( x )2 y H ( x) 10 G ( x) P n ( x) P n ( x) Y se obtiene la formula de Laguerre, despejando ( x-q) de 8 y substituyendo en 9 y resolviendo para ( x-r ).
x r
n G ( x)
(n 1)nH ( x) G 2 ( x)
El procedimiento es: 1.- Proponer un valor de prueba x como primera raíz del polinomio
11
MATLAB para novatos
54
2.- Evaluar el polinomio, su derivada y su segunda derivada en x. Como se propone aquí o como puedan. 3.- Calcular G y H, con eq. 10, solo se necesita la función y sus dos derivadas. 4.- Determinar r con la ecuación 11. Eligiendo la magnitud más grande en el denominador. 5.- r es x y repetir desde el punto 2 hasta que el error sea mínimo. Esto determina la primera raíz. 6.- con las ecuaciones 6, empezar el procedimiento para la siguiente raíz, con un polinomio de grado menor. Después de esto, queda un polinomio de grado menor y con las ecuaciones: a(1)=b(1) y b(n)=a(n)+r*b(n-1) (siendo a(n) los coeficientes del polinomio de orden n y b(n) los de orden n-1) obtenemos sus coeficientes. Repetimos el procedimiento anterior (de 1 a 5) hasta obtener la segunda raíz. Después de esto quedará un polinomio de orden n-2 y volvemos a calcular sus coeficientes. Repetimos todo el procedimiento hasta encontrar todas las raíces. Ejemplo
Encontrar las raíces del polinomio x3 2 x 2 3 x 4 0 function raiz = raices_polinomio(a,tol) % Regresa todas las raíces del polinomio de la forma % a(1)*xˆn + a(2)*xˆ(n-1) + ... + a(n+1). %Si quieres un error mínimo usa como valor de tolerancia de entrada % tol = error de tolerancia (1.0e4*eps) n=length(a)-1; raiz=zeros(n,1); for i=1:n r=laguerre(a,tol); if abs(imag(r))
MATLAB para novatos if abs(G+F)>=abs(G -F); dr=n/(G+F); else dr=n/(G-F); end r=r-dr; if abs(dr)
55
MATLAB para novatos
5.6
56
Runge Kutta
Primer orden
Una ecuación diferencial es una regla para construir funciones a partir de sus derivadas. Para obtener la función por este método numérico se plantea la ecuación diferencial de la siguiente forma: y f ( x, y ) Si la ecuación diferencial es de la forma y ( n) f ( x, y, y, y,, y ( n 1) Transformarla usando la notación y3 y ; yn y ( n 1) y1 y ; y2 y ; Para obtener y y2 y1 1 1 y 2 y3 y2 2 ; con condiciones iniciales y3 3 en x= xini y4 y3 f ( x, y , y ,, y ) yn n 1 2 n yn En algunas situaciones, como el caso de sistemas de ecuaciones diferenciales de orden mayor a uno, es necesaria hacer algo parecido para cada ecuación diferencial. Finalmente se tienen tantas ecuaciones diferenciales como la suma de los órdenes de las ecuaciones de que se parte. La construcción de la función ( x, y) se hace partiendo de la condición inicial y por medio de la definición de la derivada o Taylor a primer orden dy y ( x h) y ( x) lim dx
h0
h
Que en la aproximación de h pequeño o Taylor truncado da y( x h) y ( x) h y( x)
RK-primer orden
Método R-K de cuarto orden K 1 hF ( x, y1 , y2 ,, yn ) K 2
h
K
h
K
hF( x , y 1 ) 2 2
hF( x , y 2 ) 2 2 K 4 hF ( x h, y K 3 ) K 3
y ( x h) y ( x) Ejemplo
1 (K 2K 2 2K 3 K 4 ) 6 1
RK-cuarto orden
MATLAB para novatos
57
Un ejemplo para ilustrar este método se puede obtener del siguiente caso: El movimiento newtoniano en coordenadas polares es: 2r ) r ˆ a (r r 2 )eˆr (r
r
Vo=velocidad inicial
θ
H=altura inicial
r r 2
GM r 2
2r r
R=Radio de la tierra Figura 1. Esquema de un objeto de masa mucho menor que la tierra que se desplaza por efecto de la atracción gravitatoria. De donde se propone el siguiente conjunto de ecuaciones: y1 r y r 2 y3 y4 Y se obtiene el conjunto de ecuaciones diferenciales siguientes: y2 y 1 y1 y42 GM y2 y12 y4 y3 2 y y y 2 4 4 y1 Para las que podemos emplear las siguientes condiciones iniciales: Posición radial inicial, y1(0)=R+H=(6378.14 + 772) × 103= 7. 15014 × 106 m Velocidad radial inicial, y2(0)=0 m/s Posición angular inicial, y3(0)=0 rad Velocidad angular inicial, y4(0)= Vo/r(0) = 6700 /(7.1501 × 106 ) = 0.937045 × 10−3 rad/s Este archivo se llama “ruku.m” %Método de Runge Kutta de cuarto orden %para resolver ecuaciones diferenciales %con condiciones iniciales. %Esta función tienen que ser modificada %y agregar las condiciones iniciales para %t(ini) y para la función y(t=ini). %La función y se obtiene de un arreglo de %ecuaciones diferenciales de primer orden. % [Y’]=[F(t,y)]
MATLAB para novatos
function [B]=ruku() [t,Y,h]=CI(); y=Y(1,:); for m= 2:length(t); K1=h*yprima(t(m-1),y); K2=h*yprima(t(m-1)+h/2,y+K1/2); K3=h*yprima(t(m-1)+h/2,y+K2/2); K4=h*yprima(t(m-1)+h,y+K3); y=y+(K1+2*K2+2*K3+K4)/6; Y(m,:)=y; end B=[t' Y];
function [t,y,h]=CI() tini=0;%CAMBIAR t INICIAL tfin=1200;%CAMBIAR t FINAL h=50;%CAMBIAR PASO t=tini:h:tfin; y=zeros(length(t),4); %CAMBIAR COLUMNAS CON VARIABLES y(1,:)=[7.1501e6 0 0 9.3704e-004]; %CAMBIAR Y INICIAL
function [yp]=yprima(t,y) yp=zeros(1,4); yp(1)=y(2); yp(2)=y(1)*y(4)^2-3.9860e14/y(1)^2; yp(3)=y(4); yp(4)=-2*y(2)*y(4)/y(1);
58
MATLAB para novatos
5.7
59
Descomposición de LU
Es un método muy utilizado en programación para resolver sistemas de ecuaciones lineales, tiene la gran ventaja de ser un método fácil y rápido. La descomposición LU se basa en el proceso de eliminación de gauss. Es conveniente usar descomposición de LU cuando se requiere resolver un sistema de ecuaciones lineales de la forma Ax=b que requiere varias iteraciones. Supongamos que la matriz A se puede descomponer como el producto de dos matrices triangulares, una llamada L (lower triangular) y otra llamada U (upper triangular). Existen varias posibilidades para este producto, aquí se hablara de la descomposición de Doolittle
a11 a 21 a31 am1
a12
a13
a22
a23
a32
a33
am 2
am 3
a1n x1
a3n amn a2n
b1 x b 2 2 x3 b3 xn bn
El cual representan el sistema lineal Ax=b, y que se puede escribir LUx=b, ya que podemos representar a A como un producto de una cierta matriz triangular inferior (L) y otra matriz triangular superior (U), de la forma: [ L][U]=[A] Donde 1 0 0 0 u11 u12 u13 u1n l 1 0 0 0 u u u 21 22 23 2n L l 31 l 32 1 0 U 0 0 u33 u3n l m1 l m 2 l m3 1 0 0 0 umn Observe que los elementos de la diagonal principal de L siempre son iguales a 1 (si la matriz A es cuadrada), y siempre tiene elementos iguales a cero en la parte triangular superior de la matriz. U siempre tiene elementos iguales a cero en su parte triangular inferior. Los demás valores de los elementos de las matrices L y U se obtienen mediante la eliminación de gauss. Al aplicar el método de la eliminación de gauss y escalonar la matriz A, esta pasa a ser la matriz triangular superior U. En la matriz L se extraen de manera ordenada los coeficientes que se utilizaron para escalonar la matriz, por ejemplo, nos interesa que los
MATLAB para novatos
60
elementos de las posiciones a21, a31, de la matriz A, sean ceros, esto lo podemos lograr sustrayendo un múltiplo adecuado de las filas de A renglon2 = renglon2 - m21 * renglon1, donde m21 = a21/a11 renglon3 = renglon3 - m31 * renglon1, donde m31 = a31/a11 Estas operaciones hacen cero los elementos a21 y a31 respectivamente. Por lo tanto en la matriz triangular inferior L, asignamos en la posición L21 el valor del coeficiente m21 que se utilizó para hacer cero el elemento de la posición a21, y en L31 asignamos de igual manera el valor de m31,
1 0 0 m 21 1 0 L m31 l 32 1 l m1 l m 2 l m3
0 0 0 1
Y así sucesivamente para los demás elementos de L, hasta escalonar A y obtener todos los valores de la matriz triangular inferior L. Una vez que obtenemos LU, descomponemos A en el sistema de ecuaciones [L][U][x] = [b] Proponemos [U][x]=[z] sustituimos [L][z]=[b] Donde conocemos los valores de [L] y del vector [b], y se obtienen los valores de [z] haciendo una sustitución hacia adelante.
1 0 0 l 1 0 21 l 31 l 32 1 l m1 l m2 l m3
0 0 0 1
z 1 b1 z b 2 2 z 3 b3 z n bn
Posteriormente podemos conocer el valor de [x], por medio de una sustitución hacia atrás, de [U][x]=[z] ya que se conoce a [z] y a [U].
MATLAB para novatos
61
u11 u12 u13 0 u u 22 23 0 0 u33 0 0 0
u1n x1
u3n u mn u 2n
z 1 x z 2 2 x3 z 3 xn z n
A continuación se presenta un algoritmo de programación en Matlab, sobre una manera de obtener la descomposición LU de una matriz cuadrada A=round(rand(5)*10); %Definimos una
matriz arbitraria cuadrada de números enteros aleatorios del 1 al 10 n=length(A(:,1)); %Se asigna a la variable n la longitud de las %columnas de A L=zeros(n);%Se reserva el espacio para L y U, para que no %estén creciendo dentro de los ciclos U=zeros(n); for k = 1:n %Este es nuestro ciclo principal L(k, k) = 1; %Elementos de la diagonal principal de L %iguales a uno for i = k+1:n L(i,k) = A(i, k) / A(k, k); %Se guardan de manera ordenada %en L los coeficientes con los que vamos a escalonar la %matriz A for j = k+1:n A(i, j)= A(i, j) - L(i, k)*A(k, j); %Se escalona la matriz A end end for j = k:n U(k, j) = A(k, j); %Al terminar la matriz A, pasa a ser la matriz U end end
MATLAB para novatos
5.8
62
Interpolación suave (splines)
Se tienen n pares de puntos. De donde se tiene que son n-1 segmentos Asignando una cúbica en cada segmento, con 4 incógnitas por segmento Produce 4(n-1) incógnitas totales De cada segmento conocemos sus extremos 2(n-1) incógnitas conocidas La pendiente es continua en cada unión de segmentos (n-3) segmento se pegan por ambos lados y 2 segmentos (el primero y ultimo) se pegan en un lado. 2(n-3)+2 incógnitas conocidas La curvatura del primer punto y el último se fijan a cero. Así se reúnen todas las incógnitas: 2(n-1)+2(n-3) +2 +2= 4n-4 Función de partida, donde las k son las curvaturas de los puntos. k ( x xi1 ) k i1 ( x xi ) f i,i 1 i ( xi xi1 ) ( xi 1 xi ) 1 Integrar dos veces y evaluando en los extremos k i 1 ( x xi ) 3 k i ( x xi 1 ) 3 ( )( ) f i ,i 1 ( x xi 1 )( xi xi 1 ) x x x x ( x x ) i i i 1 ... 6 ( xi xi 1 ) 6 i i 1 y ( x xi 1 ) y i 1 ( x xi ) i ( xi xi 1 ) 2 Para el segmento anterior k i ( x xi 1 ) 3 k i 1 ( x xi ) 3 ( )( ) ( )( ) f i 1, i x x x x x x x x ( x x ) i i 1 i i 1 i 1 i ... 6 ( xi 1 xi ) 6 i 1 i y ( x xi ) yi ( x xi 1 ) i 1 ( xi 1 xi ) Derivando una vez y conectando segmentos k i 1 3( x xi ) 2 yi yi 1 k i 3( x xi 1 ) 2 ( ) f i ,i 1 ( xi xi 1 ) x x ( x x ) i i 1 6 ( xi xi 1 ) 6 i i 1 ( xi xi 1 )
k i 3( x xi 1 ) 2 yi 1 yi k i 1 3( x xi ) 2 ( ) ( ) f i 1, i x x x x ( x x ) i 1 i i 1 i 6 ( xi 1 xi ) 6 1 i i ( xi 1 xi ) yi 1 yi yi yi 1 ( xi 1 xi ) ( xi xi1 )
k i 1 ( xi 1 xi ) 2k i ( xi 1 xi1 ) k i 1 ( xi xi 1 ) 6
3
MATLAB para novatos
63
Y como conocemos que k 0 0 y k n 0 Entonces la solución para las k son: 0 0 2( x1 x3 ) ( x2 x3 ) k 2 6 y1 y2 y2 y3 k 1 ( x1 x2 ) x1 x2 x2 x3 ( ) 2 ( ) ( ) 0 x x x x x x k 2 3 3 2 4 3 4 y y y y 6 2 3 3 4 0 x x x x 2 3 3 4 ( x3 x4 ) 2( x3 x5 ) 0 k 4 y y y y 6 3 4 4 5 x x x x 4 5 3 4 y yn 1 yn 1 y n 0 0 0 2( xn 2 xn ) k n1 6 n 2 ( ) k x x xn2 xn 1 xn1 xn n n1 n lo que se simplifica si los puntos están igualmente espaciados. De este último sistema de ecuaciones se obtienen las k . De la ecuación 2 se obtiene f Y con esta se obtiene el valor interpolado. function z=interpol % x=[-5 -4 -3 -2 -1 0 1 2 3 4 5]; % y=[0 .01 .05 .25 .7 1 .7 .25 .05 .01 0]; x=[1 2.5 3 5 8]; y=[1 2 2 0 0]; % %Plantear la matriz tri diagonal [a,b,c,d]=coeficientes(x,y); a, b, c, d % %Resolver por LU la matriz tri diagonal [k]=seconder(a,b,c,d); k=[0;k;0]; % % gráfica de la solución m=length(x); mmm=0; fx=zeros((x(m)-x(1))/.1,1); f=zeros((x(m)-x(1))/.1,1); for n=1:m-1 xx=x(n):.1:x(n+1); mm=length(xx); fx(mmm+1:mmm+mm)=xx; f(mmm+1:mmm+mm)=k(n)/6*((xx-x(n+1)).^3/(x(n)-x(n+1))-(xxx(n+1))*(x(n)-x(n+1))) ... -k(n+1)/6*((xx-x(n)).^3/(x(n)-x(n+1))-(xx-x(n))*(x(n)x(n+1)))... +(y(n)*(xx-x(n+1))-y(n+1)*(xx-x(n)))/(x(n)-x(n+1)); mmm=mmm+mm; end
MATLAB para novatos
64
plot (x,y,'*') hold on plot (fx,f)
%Obtener los coeficientes de la matriz tri diagonal. Esta manera es muy %económica en cálculo y muy común en los problemas de recurrencia. function [a,b,c,d]=coeficientes(x,y) m=length(x); a=zeros(m-3,1); b=zeros(m-2,1); c=zeros(m-3,1); d=zeros(m-2,1); for n=2:m-2 a(n-1)=x(n)-x(n+1); b(n-1)=2*(x(n-1)-x(n+1)); d(n-1)=6*(y(n-1)-y(n))/(x(n-1)-x(n))-6*(y(n)-y(n+1))/(x(n)-x(n+1)); end b(m-2)=2*(x(m-2)-x(m)); d(m-2)=6*(y(m-2)-y(m-1))/(x(m-2)-x(m-1))-6*(y(m-1)-y(m))/(x(m-1)-x(m)); c=a;
% solución del sistema de ecuaciones empleando LU function [k]=seconder(a,b,c,d) m=length(b); for n=2:m lambda=a(n-1)/b(n-1); b(n)=b(n)-lambda*c(n-1); a(n-1)=lambda; end for n=2:m%Forwardsubstitution d(n)=d(n)-a(n-1)*d(n-1); end d(m)=d(m)/b(m);%Backsubstitution for n=m-1:-1:1 d(n)=(d(n)-c(n)*d(n+1))/b(n); end k=d;
MATLAB para novatos
65
6 Gráficos Las capacidades de trabajar los gráficos es una de las razones por lo que muchas personas se acercan a MATLAB, los aspectos mas básicos son muy simples, pero además se pueden hacer cosas espectaculares. En general no se requiere de programas adicionales, y lleva los resultados hasta archivos de la mayor calidad. A estas alturas del partido espero que no se tenga dudad de las ventajas de MATLAB, gráficos es un premio adicional. 6.1
Dos dimensiones
Las ilustraciones mas simples y comunes surgen de tener coordenadas formadas por un vector de números como x y otro vector de números como y.
Figura 6.1 Ilustración simple con caracteres variados. La figura 6.1 se construyó con las siguientes instrucciones >>figure1=figure % Definir tamaño de letra en el margen y rejillas axes('FontSize',20,'XGrid','on','YGrid','on','Parent',figure1) box('on'); hold('all'); % Grafica la función de gauss, con ejes x de 0 a 6 y y de 0 a 0.41 x=0:.01:6; plot(x,1/(1*sqrt(2*pi))*exp(-(x-3).^2./2)) axis([0 6 0 0.41]); % Definir leyenda en x con tamaño 20 xlabel('x','FontSize',20); % Definir leyenda en y con tamaño 20 ylabel('f(x)','FontSize',20); % Definir el titulo title('Gauss'); % Definir caja de texto para el símbolo en x=.7 y y=.75 de uno annotation(figure1 ,'textbox','Position',[.7 .75 .1 .1],'FontSize', 20,'String',{'\mu=3 \sigma=1'}); % Crear caja de texto para la función f(x) en latex text('Interpreter','latex',... 'String','$$f(x)=\frac{1}{\sigma\sqrt{2 \pi}}exp(-{\frac{(x\mu)^2}{2\sigma^2}})$$', ... 'Position',[1 .05],... 'FontSize',16)
MATLAB para novatos
66
% Crear caja de texto para la integral en latex text('Interpreter','latex', ... 'String','$$\int_{-\infty}^{\infty} f(x) dx=1$$', ... 'Position',[2 .15],... 'FontSize',16)
No desesperar, no se requiere de recordar nada de esto para que nos salga bien la figura. Si se escribe las coordenada x y y se pide que se grafiquen con plot aparece la figura 6.2. >>x=0:.01:6; plot(x,1/(1*sqrt(2*pi))*exp(-(x-3).^2./2))
Figura 6.2 grafico simple sin edición adicional. En la figura 6 se indica dentro del círculo rojo el acceso directo para Show Plot Tools, que es la entrada al paraíso de las graficas. El círculo azul muestra el acceso directo a los datos de la grafica, se puede ver las coordenadas de los puntos que forman la grafica.
MATLAB para novatos
67
Figura 6.3. Diversa opciones de edición para los gráficos, en la parte inferior de las figuras se muestran las posibilidades. A) edición de los ejes, leyendas, límites y escalas; además del titulo y de las rejillas. B) edición de la apariencia del trazo, etiquetas, líneas marcadores y colores. C) edición de una caja de texto, en este caso nombre del eje vertical, color, tamaño y formato. La figura 6 muestra ejemplos de las posibilidades de edición para modificar la apariencia de las figuras. Observe que las marcas en la figura (pequeños cuadros) coinciden con las posibilidades de edición. Además de las opciones evidentes, se pueden acceder a mas detalles de control con el acceso directo More Properties… que se encuentra en el ángulo inferior derecho. También en los menús superiores se puede tener acceso cambios adicionales o a agregar características adicionales a la figura. Otra maravilla de la edición es que en las cajas de texto, nombres de ejes, leyendas, etc se pueden mezclar caracteres, por ejemplo secuencia \alpha \beta \gamma \delta \delta \pi \lambda
símbolo
secuencia \Omega \Phi \Psi \Delta \Sigma \sim \0
símbolo
MATLAB para novatos \surd
68
√
\nabla
Más aún, se pude escribir ecuaciones complejas como las ilustradas en la figura 6.1. text('Interpreter','latex','String','$$f(x)=\frac{1}{\sigma\sqrt{2 \pi}}exp(-{\frac{(x-\mu)^2}{2\sigma^2}})$$')
Si se desea que un trazo permanezca al hacer un nuevo trazo la instrucción es hold on. >> >> >> >> >>
t = (-pi:pi/100:pi); plot(t,sin(t),'r') hold on plot(t,cos(t),'k')
Si desea tener un panel de figuras, 2 por ejemplo, en un renglón y dos columnas: >> subplot (1,2,1); plot(t,sin(t),'r') >> subplot (1,2,2); plot(t,cos(t),'k')
figure semilogx semilogy loglog
6.2
abre una nueva ventana de figura grafica con la escala en x logarítmica grafica con la escala en y logarítmica grafica con la escala en x-y logarítmica
Tres dimensiones
En virtud de que las graficas en MATLAB se construyen uniendo los puntos lo primero que se debe construir es la malla x-y para evaluar f(x,y). >> >> >> >>
x = 0:pi/50:pi; y = 2*x; [X,Y] = meshgrid(x,y); surf(X,Y,sin(X.^2+Y))
Si solo se desea ver la malla se puede ver con: >> figure; plot(X(:),Y(:),'ko')
Cuando las curvas son paramétricas, la secuencia de puntos construye una curva, se pude usar plot3. La cuarta instrucción mantiene el aspecto de cubo. >> >> >> >>
t = 0:pi/50:10*pi; plot3(sin(t),cos(t),t) grid on axis square
Instrucciones convenientes para explorar funciones son:
MATLAB para novatos >> >> >> >> >> >>
69
ezplot('x^2-y^4') ezsurf('real(atan(x+i*y))') ezcontour('sqrt(x^2 + y^2)') ezmesh(' x.*exp(-x.^2-y.^2)') ezplot3('s/2','2*s^2','3*s^3') ezpolar('1+cos(t)')
6.3
Guardar imágenes
La manera mas conveniente de guardar las imágenes es en *.fig ya que pueden editarse. Cuando se esta satisfecho de la apariencia se puede guardar en varios formatos *.jpg es el menos recomendable, compacto y fácil de leer pero de baja calidad y no es modificable. *.eps es mas aceptado Word no es un buen interprete pero aun así es muy bueno. *.emf es otra alternativa de buena calidad. Yo recomendaría guardar tres versiones, tener la versión de *.fig , la *emf para trabajar y *.eps para la presentación final. 6.4
Ejercicios
1. Graficar un orbital atómico tipo f para el átomo de hidrogeno 2. Graficar 4 funciones en una sola hoja con dinámica de interferencia Soluciones: 1. t = 0:pi/50:pi;%definición del libro f = 0:pi/80:2*pi;%def del libro [T,F] = meshgrid(t,f); %construcción de puntos completos Y=(5*cos(T).^2-1).*sin(T).*exp(-1i*F); %L=3,M=1 [x,y,z] = sph2cart(F,pi/2-T,abs(Y.^2)); % recordar que matlab y el libro definen los ángulos de manera diferente, T=F y F=90-T mesh(x,y,z)
2.
%iNTERFEROMETRO DE LUZ BLANCA % %Un pulso es fijo y el orto de mueve de 4 unidades de tiempo a -4 unidades %de tiempo %la fase dentro del pulso se mantiene % %tic clear all uni=8;%unidades de tiempo num=100;%numero de puntosa calcular 100 es rapido, 1000 es preciso pul=1;%1/ancho del pulso, .000001 es sin pulso casi for delta=1:num; t=linspace(0,10,1000); e1=1*exp(-pul*(t-5).^2).*1.*cos(10*(1*t-5)); e2=1*exp(-pul*(t-9+(delta-1)*uni/num).^2).*1.*cos(10*(1*t-9+(delta1)*uni/num)); Int=(e1+e2).^2;% exponente 2, 4 o 6 para primero, segundo y tercer orden Int2(delta)= norm(Int,1);
MATLAB para novatos
%end %toc figure (1) subplot (2,2,1) plot (t,e1) subplot (2,2,2) plot (t,e2) subplot (2,2,3) plot (t,Int) subplot (2,2,4) plot (Int2) end plot (linspace(0,uni,num),Int2)
70
MATLAB para novatos
71
7 Prácticas de buena programación La eficiencia de la programación en MATLAB se puede veneficiar empleando criterios que emplean extensivamente las fortalezas del leguaje, en particular el uso extensivo de matrices.
Las líneas más importantes para mantener la eficiencia se pueden agrupar en: 1. Favorecer vectores en lugar de bucles. 2. Evitar que los arreglos crezcan mientras se corre el programa. 3. Escribir MEX-File cuando los bucles son inevitables 4. Evitar el uso innecesario de scripts y cambios innecesarios de variables 5. Si solo se hace calculo real, evitar usar funciones para números complejos cuando existen alternativas para números reales. 6. Usar los operadores lógicos correctamente 7.1
Favorecer vectores en lugar de bucles.
Para las personas que aprendieron a programar con ciclos repetitivos (bucles) es natural buscar como implementarlos en MATLAB y es muy simple pero en general puede ser ineficiente. En el siguiente ejemplo se generan 100 índices m sobre los que se construye las variables t, y. for m = 1:100 t(m) = 3*pi*(m-1)/100; y(m) = sin(t(m)); end
la mejor manera de hacer lo mismo es por medio de vectores t = 0:pi/100:3*pi; y = sin(t);
La primera línea genera el vector renglón t que contiene valores desde 0 hasta 3π, en incrementos de π/100. La segunda línea genera para cada valor de t un valor de y, esto es un vector renglón con los mismos elementos que t. En mi computadora la versión vectorial es 19 veces más rápida que la del bucle. Para ver un ejemplo de la implementación de vectorización en las funciones de MATLAB se pueden ver la funcione repmat en el directorio de MATLAB; en mi caso en C:\Program Files\MATLAB\R2006a\toolbox\matlab\elmat.
MATLAB para novatos 7.2
72
Evitar que los arreglos crezcan mientras se corre el programa.
Al permitir que un arreglo crezca dentro de un bucle reduce el tiempo de ejecución, en el siguiente ejemplo la variable x aumenta dentro del bucle de tener un elemento a 10000 elementos. tic x = 0; for k = 2:10000 x(k) = x(k-1) + 5; end toc
En esta nueva versión se define de antemano el tamaño de x como de 10000 elementos, todos cero y se cambian por le valor deseado dentro del bucle. tic x = zeros(1, 10000); for k = 2:10000 x(k) = x(k-1) + 5; end toc
Este simple cambio permite que la segunda versión sea 730 veces más rápida. 7.3
Escribir MEX-File cuando los bucles son inevitables
Este es un comentario que tiene que ser juzgado por el usuario de MATLAB, es cierto que MEX-Files (véase sección 9) no son interpretados por MATLAB ya que son equivalentes a las funciones escritas en C, el problema pudiera ser que no se dispone del tiempo para implementar el código apropiado en cuyo caso la promesa de velocidad se ve limitada por la lentitud de crear el programa. El la sección 9 daremos algunos puntos para acelerar los programas que pueden incorporarse a la base del conocimiento en niveles posteriores, cuando MATLAB es mas familiar y se desea explotarlo mas eficientemente. 7.4
Comentarios de graficas activas
En la sección de control reflujo (4.1), se pide en un ejercicio que se dibuje un helecho y la manera de presentar el gráfico es un poco extraña, qui se explica la razón del código al comparar la velocidad con que se grafica datos en dos caminos alternativos. Esta primera alternativa solo recuerda el ultimo número, y sin menosprecio a que guarde los números, es muy rápida, si se trata de editar, sólo tiene el último número graficado function []=espiral() tic x=[0;1]; %guarda el manejador de la gráfica en h
MATLAB para novatos
73
h=plot (x(1),x(2),'.'); %define en la grafica el color a graficar y no borra los datos set(h,'markersize',1,'color','red','erasemode','none'); %define el color del fondo, el titulo de la garfica set(gcf,'color','white','menubar','none','numbertitle','off','name','es piral') axis([-3 3 -2 2]) %borra los ejes de la figura axis off drawnow %espiral simple de 1000 puntos for m=1:1000 x(1)=exp(-2*m/1000)*sin(m*6*pi/1000); x(2)=exp(-2*m/1000)*cos(m*6*pi/1000); set(h,'xdata',x(1),'ydata',x(2)); drawnow end toc
Esta segunda alternativa es muy lenta, recuerda cada uno de los números que se han graficado y como la gráfica tiene toda la información es muy lenta al hacer espacio, es como crecer el arreglo mientras esta corriendo el programa. function []=espiral2() tic x=[0;0]; plot (x(1),x(2),'.') axis([-3 3 -2 2]) hold on for m=1:1000 x(1)=exp(-2*m/1000)*sin(m*6*pi/1000); x(2)=exp(-2*m/1000)*cos(m*6*pi/1000); plot (x(1),x(2),'.') drawnow end toc
7.5
Otras consideraciones
Si solo se trabaja con números reales, existen funciones para números reales que mejoran la eficiencia de los cálculos, cuando existe un uso masivo de ellas. log sqrt
pow
reallog
realpow
realsqrt
La ejecución de funciones es mas eficiente que la ejecución de scripts, estas ultimas son leídas y ejecutadas una línea a la vez, mientras que las funciones permanecen en memoria después de ser usadas la primera vez.
MATLAB para novatos
74
Los procesos compartidos en las computadoras pueden disminuir el desempeño, cuando los cálculos son muy intensivos, evite compartir los recursos de cómputo entre varias actividades. 7.6
Archivos
Ahora presentaremos tres temas; como guardar archivos mientras corre el programa, como leer archivos mientras se corre el programa y finalmente como extraer datos de un archivo .fig. Suponiendo que se tiene una matriz R con la información que se desea guardar, en este caso tres vectores de la misma longitud. Como esta instrucción se tomo de otro programa y solo tiene el fin de ilustrar el procedimiento no se espera que tenga sentido o que corra si se incluye en un script de matlab, tendría que modificarse esta matriz por una que estuviera bien definida. R=[real(t), real(y(:,4)+y(:,3)+y(:,1)+V1), real(V1./RS)]; m=m+1; fname=strcat('C:\Users\bigcheese\Documents\test4\' ,'data',num2str(m)); save (fname,'R','-ASCII')
En este caso el folder (’C:\Users\bigcheese\Documents\test4\' ) debe existir en la computadora y genera en el un nuevo archivo de nombre datam donde la m se sustituye por el número que tenga la variable m. En las siguientes lineas se forma el nombre de la direccion don de se encuentra el archivo y al leerlo se graba en la variable R, la siguiente linea pregunta por la longitud de las columnas de R. fname=strcat('C:\Users\bigcheese\Documents\test4\' ,'data',num2str(m)); R=load (fname,'-ASCII'); aa=max(R(:,1));
Finalmente vamos a extraer los valores de las variables de un archivo *.fig Archivos *.fig son binarios que contienen toda la información que reconstruye la gráfica, incluyendo los datos que la generan. Recordemos que el archivo puede tener mucha información oculta, tal que muchos datos no se muestran en el trazo, lo que no significa que esos datos se hayan perdido, por el contrario, puede ser muy complicado extraer los datos si la figura es muy compleja o tiene muchos datos ocultos. Primero suponemos que el archivo se encuentra en el directorio base de MATLAB, y se incluye en la variable a. toda la información esta ahora en a, de esta hay que buscarla y tomar los datos, en este caso x esta en la variable que se enlista a continuación, al igual que la variable y. Se puede explorar la variable a en el editor de variables, en general si la información es simple es muy probable que siga la secuencia que aquí se muestra. Los datos están en a, solo hay que encontrarlos. a = load('seno.fig', '-mat') >> x=a.hgS_070000.children(1,1).children(1,1).properties.XData;
MATLAB para novatos >> y=a.hgS_070000.children(1,1).children(1,1).properties.YData;
7.7
Ejercicios
1. Construir un arreglo de coordenadas en tres dimensiones de 5 unidades en cada dimensión sin emplear bucles. Crear una alternativa al programa clear all Lx=50; Ly=50; Lz=50; A=zeros([Lx*Ly*Lz,3]); for indicex=1:Lx for indicey=1:Ly for indicez=1:Lz A(indicez+(indicey-1)*Ly+(indicex-1)*Lz*Ly,:)=[indicex indicey indicez]; end end end
Posible Solución Lx=50;Ly=50;Lz=50; X=[repmat(1:Lx,Ly*Lz,1)]; Y=[repmat(1:Ly,Lz,Lx)]; Z=[repmat(1:Lz,1,Lx*Ly)]; INDX=[X(:),Y(:),Z(:)];
75
MATLAB para novatos
76
8 GUIs La idea de darle datos de entrada a un programa y esperar por la respuesta es la noción secuencial más común de programar. La utilidad del programa para otros existe cuando es amigable con el usuario, cuando indica qué es razonable proponer al programa, los botones, menús y valores predeterminados son muy útiles y cambian el estilo del programar a activado por eventos. GUI es el acrónimo de Graphics User Interface, Interfase gráfica para el usuario.
Usualmente el proceso es, pensar en construir el programa activado por eventos, pero construir las partes secuenciales del programa y probarlas antes de iniciar la construcción del GUI. Las correcciones posteriores son poco económicas, al grado de que es mejor iniciar de cero que corregir una estructura terminada.
8.1
Ejemplo de script para FFT Ejemplo de function para FFT Ejemplo de GUI para FFT Ejemplo de generar ejecutable de FFT Ejemplo de script para FFT
La siguiente secuencia de instrucciones define t, f, su transformada y sus gráficas. t=5:.01:15; y=5*cos(2*pi*10*t)+3*cos(2*pi*25*t); tmin=t(1); tmax=t(end); dt=t(2)-t(1); N= length(t); M=1024 Y=fft(y,M); r=abs(Y(2:M/2))*2/N; o=angle(Y(2:M/2))*2/N; f=(1:(M-2)/(M-4):M/2)/dt/M; figure (1) plot(t,y) title('oscilograma') grid on xlabel('tiempo') title('funcion') figure (2) plot(f,r) title('espectrograma') grid on xlabel('frecuencia') title('amplitud')
explicar el programa…
MATLAB para novatos 8.2
77
Ejemplo de function para FFT
Se tiene que definir la function ffcorrecta y llamarla después de tener los vectores (t, y), las gráficas se hacen fuera en otro script . % % % % % % % % % %
[r, o, f]=ffcorrecta(t, y, M) salidas de [amplitud, fase, frecuencia] entradas de (tiempo, funcion, M) el tiempo pude ser un vector como t=5:.01:15 la función puede ser y=5*cos(2*pi*10*t)+3*cos(2*pi*25*t) M es un entero igual o mayor que número de puntos discretos, de no introducirse, entonces M es el número de puntos discretos.
function [r, o, f]=ffcorrecta(t, y, M); if (nargin == 2) M=0; end tmin=t(1); tmax=t(end); dt=t(2)-t(1); N= length(t); if M
explicar las variantes del programa 8.3
Ejemplo de GUI para FFT
Existen dos maneras de crear una GUI, aquí presentaremos la mas fácil, emplear GUIDE como asistente. 1. escribir guide en command prompt 2. seleccionar la hoja en blanco, ver figura 7.1
MATLAB para novatos
78
Figura 8.1 Pantalla de decisión después de llamar guide. 3. Agregar los objetos que se deseen, ver figura 8.2 rectángulo rojo. Como ejemplo hagamos el ejercicio para FFT. Los rectángulos negros indican los textos dinámicos, que cambian en la ejecución del programa. El rectángulo verde indican los textos estáticos, que no se modifican durante la ejecución del programa. El rectángulo púrpura indica una barra deslizable, el rectángulo azul un botón y los rectángulos amarillo las regiones de las graficas.
Figura 8.2 pantalla de la figura donde se distribuyen los objetos que controlan la ejecución del programa. El rectángulo rojo indica el acceso directo a los objetos. En la figura 8.2 para cada objeto que se agregó, se selecciona con el botón derecho properties inspector para modificar tag (leyendas que aparecen en el programa) y string (leyenda que aparecen en la figura). Para nuestro ejemplo se han agregado y modificado las leyendas como se ve en la tabal 8.1.
Texto estático
Tabla 8.1 Leyendas en la figura tag sugerida tag definida string sugerida text1 text1t Static text text2 text2titulo Static text text3 text3y Static text text4 text4M Static text text5 text5instrucc Static text
string definida t= FFT y= valor de M Instrucciones…
MATLAB para novatos Texto dinámico
Barra deslizable Botón Gráfica
edit1 edit2 edit3 slider1 pushbutton1 axes1 axes2
79 edit1t Edit text edit2M Edit text edit3y Edit text slider1M button1START push button 1 axesoscilo axesespectro
5:.01:15 Edit text 5*cos(2*pi*10*t)… START
4.- Al terminar se salva con un nombre significativo y se generan dos archivos, *.fig y *.m. El archivo fig tiene toda la información del desplegado de los objetos, el archivo m tiene la programación de los eventos, véase la figura 8.3.
Figura 8.3 archivos *.fig y *.m generados al salvar el desplegado de los objetos en la figura. La ilustración de la derecha es después de presionar el botón START y haber terminado de programar el archivo *.m. La ilustración de la derecha muestra un código generado por guide al que hay que agregar lo que cada acción en los objetos realizara. Para concluir el programa solo es necesario agregar las funciones deseadas. Las modificaciones se agregan en las acciones de los objetos. En este programa solo es necesario decir lo que sucede al presionar el botón START, y la manera en que se interrelacionan el control deslizable y el texto dinámico que definen a la variable M. Adicionalmente es necesario poner atención en el intercambio de datos dentro de la GUI, de los textos dinámicos al botón START y del este botón a las graficas.
Modificación al pushbutton callback t=eval(char(get(handles.edit1t, 'String'))); y=eval(char(get(handles.edit3y, 'String'))); tmin=t(1); tmax=t(end); dt=t(2)-t(1); N= length(t); M=floor(get(handles.slider1M, 'Value')*2048); Y=fft(y,M); r=abs(Y(2:M/2))*2/N;
MATLAB para novatos
80
o=angle(Y(2:M/2))*2/N; f=(1:(M-2)/(M-4):M/2)/dt/M; plot(handles.axisoscilo,t,y) grid (handles.axisoscilo, 'on') xlabel(handles.axisoscilo, 'tiempo') ylabel(handles.axisoscilo, 'funcion') plot(handles.axisespectro,f,r) grid (handles.axisespectro, 'on') xlabel(handles.axisespectro, 'frecuencia') ylabel(handles.axisespectro, 'amplitud')
La línea 1 crea el eje t de los caracteres en edit1t. La línea 2 crea a la función de los caracteres en edit3y. La línea 7 define el valor de M del valor que tenga slider1M. Las ultimas 8 líneas definen donde de realizarán las acciones en las graficas.
Modificación al slider callback M=floor(get(hObject, 'Value')*2048); letrasM=num2str(M); set(handles.edit2M, 'string',letrasM)
Estas líneas toman el valor de slider y lo mandan a edit2M
Modificación al editM callback M=str2num(get(handles.edit2M, 'String')); if (M>=(get(handles.slider1M, 'Min')*2048)) & (M<=(get(handles.slider1M, 'Max')*2048)) set(handles.slider1M, 'Value',M/2048) else set(hObject,'String','error') end
Estas líneas leen el valor que se introduce como texto, lo compara con los valores posibles para slider , si esta contenido escribe el número y lo transfiere a slider , de no ser, escribe error. 8.4
Pasarlo a ejecutable
Con el MATLAB compiler se puede escribir mcc -m ffcorr.m se selecciona el compilador de los disponibles, si no se tiene alguno especial, MATLAB ofrece uno propio.
MATLAB para novatos
81
9 Conexión con C El propósito principal de esta sección son los MEX-files, que permiten manejar subrutinas en C o FORTRAN dentro de MATLAB, aunque también es posible prepara material de MATLAB para ser usado por otros programas e incluso preparar aplicaciones autocontenidas como se hizo al final de la sección anterior.
Como se comento anteriormente cuando se desea acelerar algún proceso, en particular bucles, es posible pedir a MATLAB que realice esa tarea con la ayuda de los archivos MEX. Estos son alternativas a los archivos m y están conformados además de la parte de computación por una plataforma de intercambio de datos, recuerde que para MATLAB todo es un arreglo mientras que en C es necesario especificar el contenido. Las diferencias con los archivos m son: estos son dependientes de la plataforma en que serán usados (mexw32 y mexw64 en Windows), el manejo de la memoria es responsabilidad del programador y no tiene documentación así que esta tiene que ser manejada de diferente manera. Como siempre, es necesario decir que MATLAB es altamente eficiente, y los archivos MEX no necesariamente son una alternativa de velocidad, pero si la subrutina que desea usar en MATLAB ya esta diseñada y es muy eficiente, se puede usar por este método como las subrutinas disponibles en www.netlib.com. Como ejemplo en el directorio de MATLAB esta un ejemplo del programa en C y en MATLAB, (C:\Program Files\MATLAB\R2006a\extern\examples\mex) yprime.c y yprime.m, la secuencia de instrucciones para generar el yprime.mexw32 a partir de yprime.c es la siguiente. Primero mostremos el archivo m y los resultados para y=[1 2 3 4] >> cd([matlabroot '\extern\examples\mex']) mex yprime.c >>
>> y=[1 2 3 4]; mu = 1/82.45; mus = 1-mu; r1 = norm([y(1)+mu, y(3)]); % Distance to the earth r2 = norm([y(1)-mus, y(3)]); % Distance to the moon yp(1) = y(2); yp(2) = 2*y(4) + y(1) - mus*(y(1)+mu)/r1^3 - mu*(y(1)-mus)/r2^3; yp(3) = y(4); yp(4) = -2*y(2) + y(3) - mus*y(3)/r1^3 - mu*y(3)/r2^3; yp =
MATLAB para novatos
82
2.0000 8.9685 4.0000 -1.0947 Ahora mostremos el archivo c #include
#include "mex.h" /* Input Arguments */ #defineT_IN prhs[0] #defineY_IN prhs[1] /* Output Arguments */ #defineYP_OUT
plhs[0]
#if !defined(MAX) #defineMAX(A, B) #endif
((A) > (B) ? (A) : (B))
#if !defined(MIN) #defineMIN(A, B) #endif
((A) < (B) ? (A) : (B))
static double mu = 1/82.45; static double mus = 1 - 1/82.45; static void yprime( double double double ) { double r1,r2;
yp[], *t, y[]
r1 = sqrt((y[0]+mu)*(y[0]+mu) + y[2]*y[2]); r2 = sqrt((y[0]-mus)*(y[0]-mus) + y[2]*y[2]); /* Print warning if dividing by zero. */ if (r1 == 0.0 || r2 == 0.0 ){ mexWarnMsgTxt("Division by zero!\n"); } yp[0] = y[1];
MATLAB para novatos yp[1] = 2*y[3]+y[0]-mus*(y[0]+mu)/(r1*r1*r1)-mu*(y[0]-mus)/(r2*r2*r2); yp[2] = y[3]; yp[3] = -2*y[1] + y[2] - mus*y[2]/(r1*r1*r1) - mu*y[2]/(r2*r2*r2); return; } void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] ) { double *yp; double *t,*y; unsigned int m,n; /* Check for proper number of arguments */ if (nrhs != 2) { mexErrMsgTxt("Two input arguments required."); } else if (nlhs > 1) { mexErrMsgTxt("Too many output arguments."); } /* Check the dimensions of Y. Y can be 4 X 1 or 1 X 4. */ m = mxGetM(Y_IN); n = mxGetN(Y_IN); if (!mxIsDouble(Y_IN) || mxIsComplex(Y_IN) || (MAX(m,n) != 4) || (MIN(m,n) != 1)) { mexErrMsgTxt("YPRIME requires that Y be a 4 x 1 vector."); } /* Create a matrix for the return argument */ YP_OUT = mxCreateDoubleMatrix(m, n, mxREAL); /* Assign pointers to the various parameters */ yp = mxGetPr(YP_OUT); t = mxGetPr(T_IN); y = mxGetPr(Y_IN); /* Do the actual computations in a subroutine */ yprime(yp,t,y); return; }
83
MATLAB para novatos
84
Después de configurar el compilador, en caso de disponer de alguno diferente al que provee MATLAB mex ( -setup ), y de definir el directorio de trabajo (cd([matlabroot '\extern\examples\mex']) ) se pide la generación del archivo mex. >>mex yprime.c
Que produce el archivo yprime.mexw32 que puede ser usado en MATLAB >> yprime(1,1:4) ans = 2.0000
8.9685
4.0000
-1.0947
>> Para los formatos definidos para crear un archivo c para convertirlo a archivo mex véase la ayuda de MATLAB en External Interfaces.
MATLAB para novatos
85
10 Simulink Aquí se presenta una visión superficial que puede permitir iniciarse en el modulo de Simulink, que básicamente es en modulo para la solución numérica de ecuaciones, hasta sistemas de ecuaciones diferenciales. Oficialmente permite la simulación dinámica de sistemas, que pueden ser lineales y nolineales.
La figura 10.1 identifica el acceso directo a Simulink, 10.2 muestra la biblioteca de funciones de Simulink que se despliegan inicialmente, agrupadas por categorías.
Figura 10.1 Acceso directo a Simulink. 10.1 Simulink para resolver ecuaciones algebraicas.
Con Simulink se puede pedir a MATLAB que ejecute las mismas instrucciones, ahora presentadas de manera grafica. 1 3 b b 2 4ac Por ejemplo y la raíz cuadrática mayor con a=b=c=1 da y i . 2a 2 2 Se puede escribir en matlab >> a=1;b=1;c=1; >> y=(-b+sqrt(b^2-4*a*c))/(2*a) y =
MATLAB para novatos
86
-0.5000 + 0.8660i
Figura 10.2 Representación en Simulink de la raíz cuadrática mayor. A la derecha esta la ventana de edición de dos objetos; la función matemática (raíz cuadrada) que puede seleccionar otras funciones matemáticas y la ganancia que multiplica el valor que entra por el factor de ganancia. En Simulink, vea la figura 10.2, se ve de esta manera, donde los objetos se toman de la biblioteca de Simulink, ver la figura 10.3. Los objetos son claros en cuanto a que tienen entradas y salidas, las entradas son los argumentos y las salidas son los resultados. Para este tipo de ejemplos los grupos de la biblioteca mas empleados son las funciones matemáticas (Math Operations), las señales empleadas como fuentes (Sources) y los exhibidores (Sinks). Los objetos se conectan entre sus puertos y para ejecutar la secuencia se presiona el botón encerrado en el círculo verde en la figura 10.2.
Figura 10.3 Biblioteca de Simulink, la columna de la izquierda muestra las categorías.
MATLAB para novatos
87
10.2 Simulink para resolver ecuaciones diferenciales
Como ejemplo de uso del Simulink como herramienta para la solución de ecuaciones diferenciales se repite el ejemplo de la sección 5 de la siguiente ecuación diferencia.
2
d 2 f 2
dt
1
df dt
5 f 3 .05t 3 5 cos(1.5t )
con las condiciones iniciales f (0)=1.3 y
df dt 0
8.
La figura 10.4 muestra este ejemplo. Es interesante notar el ambiente grafico de trabajo en Simulink, donde las jerarquías se definen por medio de flechas y los operadores por medio de cajas que realizan las funciones deseadas. El operador derivada se muestra con la etiqueta “integrator”, la manera de trabajar en Simulink es despejar la derivada de orden mayor que en la figura 9.4 se encuentra a la derecha del cuadro rojo e igualarla a g(t) y a las derivadas de orden menor.
Figura 10.4 Editor de Simulink donde se muestra el ejemplo de una ecuación diferencial simple. Los osciloscopios son receptores de datos, en esta simulación para 10 segundos, círculo verde en la figura 10.4 y el botón para dar inicio a la simulación se encuentra señalado en la figura 10.4 por medio del círculo azul.
MATLAB para novatos
88
Figura 10.5 Simulación de los primeros 10 segundos, de izquierda a derecha se muestra la función f, su primera derivada y su segunda derivada. La figura 10.5 muestra los resultados de la simulación para 10 segundos, de la función y sus dos derivadas, note en la figura las condiciones iniciales. 10.3 Simulink para sistemas de ecuaciones
Veamos ahora un ejemplo para un sistema de ecuaciones como el siguiente, I fuente1 I fuente 2 V 1 L1 V 2 L2
V 1 R1
V 2 R2
dI L1
dI L 2
I L1
V 1 V 2
I L 2
R3 V 2
C 1
V 1
R3
dV 1
C 2
dt dV 2 dt
dt dt
10.6 Programa para el sistema de ecuaciones, los módulos se describen en la figura 10.7
MATLAB para novatos
89
10.7 Modulos empleados en la figura 10.6. En las figuras 10.6 y 10.7 se presenta la información que describe a las ecuaciones, en este caso se hizo uso extensivo de subsystem. Este modulo esta en la biblioteca de Simulink y al agregarlo al proyecto solo incluye un conector de entrada y un conector de salida, se pueden agregar más y construir la estructura tan elaborada como sea necesario. En estos ejemplos simples no se ha hecho hincapié en configuración de parámetros de la simulación, la figura 10.8 muestra la ventana que se despliega al seleccionar la pestaña de simulación en el editor de Simulink y seleccionar parámetros de configuración.
Figura 10.8 Configuración de los parámetros de la simulación. El algoritmo de solución ( solver ) se puede cambiar ode45 es el mas recomendado para iniciar, se debe consulta la ayuda para decidir sobre las opciones de simulación. Las mas solicitadas son la elección del tiempo de simulación, el delta de la variable independiente y las tolerancias del error, vease la figura 10.8.
MATLAB para novatos
90
En esta introducción no se ha hecho hincapié en la transferencia de datos entre MATLAB y Simulink, pero es muy versátil al igual que guardar los resultados.
MATLAB para novatos
11
91
SimMechanics
SimMechanics es un ambiente para modelar sistemas dinámicos que emplea diagramas a bloques. Aquí se diseñan y simulan mecanismos rígidos y su movimiento de acuerdo a las leyes dinámicas de Newton. Con SimMechanics se puede modelar y similar sistemas mecánicos con un conjunto de herramientas que permiten especificar cuerpos, su posible movimiento y sus restricciones dinámicas, el sistema mecánico se representa por medio de conexiones entre bloques. Aquí es posible exhibir y animar representaciones simplificadas de la simulación. SimMechanics trabaja en el ambiente de Simulink y es completamente compatible con el. La diferencia es que Simulink representa operaciones matemáticas mientras que SimMechanics representa componentes físicos y relaciones entre ellos. Si se quisiera representar un péndulo simple, idealizado. l
g
m
b=coeficiente que se opone a la velocidad Figura 11.1 Esquema del péndulo simplificado. Con las consideraciones de que el punto de apoyo de la cuerda no se desplaza, que el alambre que sostiene la masa no tiene peso, ni cambia de longitud, ni frena al péndulo por fricción, que la masa es puntual y que la única fuente de fricción es la rotula con el apoyo. Estas consideraciones en general se omiten en SimMechanics y se puede hacer una simulación menos idealizada. Es importante mencionar algunos preliminares; las coordenada se pueden fijar en función de un punto universal de referencia, o referidas al centro de masa o referidas al punto de unión con el elemento adjunto. La geometría se define por el momento de inercia, una esfera homogénea de radio r tiene un momento de inercia,
1 0 0 I r 0 1 0 r eye(3) 0 0 1 un objeto puntual tiene una matriz cero.
MATLAB para novatos
92
Para iniciar, crear un nuevo modelo, se parte de la biblioteca de Simulink y se transfieren los iconos del modelo que se desean. En la figura 11.2 se muestran los iconos disponibles para SimMechanics, en particular en la sección de cuerpos, los primeros tres serán comunes en todos los proyectos, Body representa cualquier objeto físico representado por masa y geometría; Ground señala el punto coordenado de referencia y Machine Environment se emplea para definir la aceleración de la gravedad.
Figura 11.2 Biblioteca de SimMechanics dentro de la Simulink.
Figura 11.3 Proyecto de un péndulo simplificado.
