Guia Rapida Para Matlab y Simulink

INTRODUCCIÓN RÁPIDA A MATLAB Y SIMULINK PARA CIENCIA E INGENIERÍA

MANUEL GIL RODRÍGUEZ Científico Titular del Consejo Superior de Investigaciones Científicas (CSIC)


© Manuel Gil Rodríguez, 2003

Reservados todos los derechos. «No está permitida la reproducción total o parcial de este libro, ni su tratamiento informático, ni la transmisión de ninguna forma o por cualquier medio, ya sea electrónico, mecánico, por fotocopia, por registro u otros métodos, sin el permiso previo y por escrito de los titulares del Copyright.»

Ediciones Díaz de Santos, S. A. Doña Juana I de Castilla, 22. 28027 MADRID

E-mail: [email protected] Internet://http:www.diazdesantos.es/ediciones

ISBN: 84-7978-596-9 Depósito legal: M. 43.971-2003

Diseño de cubierta: Angel Calvete Fotocomposición e impresión: Fernández Ciudad, S. L. Encuadernación: Rústica-Hilo, S. L. Impreso en España

Índice

Presentación ...........................................................................................

XI

1. PRIMEROS PASOS EN MATLAB ................................................ 1.1. Introducción ........................................................................... 1.2. Comenzando .......................................................................... 1.3. Espacio de trabajo .................................................................. 1.4. Variables ................................................................................ 1.5. Formato de números .............................................................. 1.6. Programas ............................................................................... 1.7. Funciones ............................................................................... 1.7.1. Reglas de construcción de funciones .......................... 1.7.2. Funciones en línea ...................................................... 1.7.3. Ejemplo de función recursiva ..................................... 1.8. Números complejos ............................................................... 1.9. Manejo de vectores y matrices ............................................... 1.10. Polinomios .............................................................................. 1.10.1. Multiplicación y división de polinomios .................. 1.10.2. Desarrollo en fracciones simples .............................. 1.10.3. Derivadas de polinomios .......................................... 1.10.4. Integración de polinomios ........................................ 1.10.5. Interpolación polinomial............................................ 1.11. matlabpath .............................................................................. 1.12. lookfor..................................................................................... 1.13. LATEX ................................................................................... 1.14. Funciones del tiempo ............................................................. 1.15. Intercambio de datos ..............................................................

1 1 1 2 3 5 7 8 9 11 11 12 13 16 16 17 17 18 18 19 19 19 20 21

VII

VIII

ÍNDICE

2. CÁLCULO SIMBÓLICO ............................................................... 2.1. Introducción ........................................................................... 2.2. Objetos y expresiones simbólicas .......................................... 2.3. Ejemplos de cálculo simbólico .............................................. 2.3.1. Derivadas e integrales ................................................. 2.3.2. Sistemas de ecuaciones ............................................... 2.3.3. Ecuaciones diferenciales.............................................. 2.4. Transformación de Laplace e inversa .................................... 2.5. Límites ................................................................................... 2.6. Series de Taylor y Mac Laurin .............................................. 2.7. Invocando a Maple V.............................................................. 2.7.1. Transformación de Laplace e inversa con Maple V .... 2.7.2. Resolución de ecuaciones diferenciales con Maple V. 2.7.3. Resolución de ecuaciones diferenciales con la transformación de Laplace ..................................................

25 25 25 27 27 28 29 30 31 31 31 31 32

3. SENTENCIAS DE CONTROL DE FLUJO ................................... 3.1. input ....................................................................................... 3.2. if - else - end .......................................................................... 3.3. while - end ............................................................................. 3.4. for - end .................................................................................. 3.5. continue .................................................................................. 3.6. break ....................................................................................... 3.7. switch end ...............................................................................

35 35 36 38 38 40 40 40

4. GRÁFICOS EN MATLAB ............................................................. 4.1. Tipos de gráficos .................................................................... 4.2. Utilidades de gráficos ............................................................ 4.3. TEXtos en gráficos ................................................................ 4.4. LaPrint ................................................................................... 4.5. Estilos de líneas, marcas y colores ......................................... 4.6. area .........................................................................................

41 41 44 45 45 47 48

5. APLICACIONES DE CÁLCULO NUMÉRICO............................. 5.1. Integración numérica ............................................................. 5.2. Mínimos, ceros y optimización............................................... 5.2.1. Mínimos y ceros de funciones .................................... 5.2.2. Resolución de ecuaciones no lineales ......................... 5.2.3. Minimización y ajuste de datos .................................. 5.3. Integración numérica de ecuaciones diferenciales ................. 5.3.1. Método de Runge-Kuta ............................................... 5.4. Modelo dinámico de un tanque .............................................. 5.5. Determinación de retrasos y derivadas .................................. 5.6. Ajuste de datos experimentales a una recta ...........................

51 51 51 51 52 53 57 59 62 67 72

33

ÍNDICE

IX

5.6.1. Ajuste de funciones no líneales por linealización ....... 5.7. Anális Espectral ..................................................................... 5.8. Evitando la división por cero y rebose ...................................

74 78 80

6. SIMULINK ...................................................................................... 6.1. Introducción a Simulink ......................................................... 6.2. Construcción de un modelo muy sencillo .............................. 6.3. Solución Simulink de una ecuación diferencial ..................... 6.4. Simulación dinámica de un ecualizador ................................

81 81 85 86 88

BIBLIOGRAFÍA .....................................................................................

93

ÍNDICE ALFABÉTICO ..........................................................................

95

Presentación

MatLab (MATrix LABoratory) comenzó siendo un programa interactivo de análisis integrado, especializado en cálculos matriciales. En muy poco tiempo tuvo una gran difusión, facilitado por su potencia de cálculo y facilidad de uso. La primera versión de MatLab, de éxito generalizado, se lanzó el 3 de febrero de 1989 para MS-DOS, Mac y Workstations, era la versión 3.5 y venía en tres discos de 5 1/4”, estaban disponibles las toolboxes de Signal Processing, Control Systems, System Identification, State-Space Identification, RobustControl, Spline, Chemometrics y Optim. MatLab incluye funciones elementales de cálculo, de Bessel, de complejos, transformadas e inversas de Laplace y Fourier, filtros para procesado de señales, Max, Min, Sum, Product, Cumulative, Cumulative Product, Mean, Median, Sort... Las funciones matriciales incluyen el determinante, el inverso, valores y vectores propios, descomposiciones, factorizaciones... Desde las primeras versiones, MatLab disponía de potentes herramientas de representaciones gráficas, ya en 2-d como 3-d, que en el transcurso de las sucesivas versiones se vieron muy potenciadas, permitiendo exportar los cuadros gráficos a otras aplicaciones. Esta obra abarca lo esencial de MatLab y Simulink, que ya no cambian con las versiones que van saliendo. Está especialmente indicada para alumnos de Ciencias, Ingeniería, Postgrado o profesionales que deseen familiarizarse, en un tiempo mínimo con las principales herramientas de cálculo suministradas por MatLab a través de ejercicios prácticos a lo largo de este manual. Manuel GIL RODRÍGUEZ

XI

1 Primeros pasos en MatLab

1.1. INTRODUCCIÓN MatLab, desde las primeras versiones dispuso de help y demo, para iniciación rápida. La información suministrada a través de los menús de estas ayudas, correspondientes a las últimas versiones, crecieron de forma exponencial, siendo de utilidad práctica disponer de un libro resumen de MatLab, en donde se encuentren los comandos de uso más frecuente, a la vez que se muestren sus aplicaciones prácticas en ejercicios, desde lo más sencillo, hasta otros de mayor complejidad. Este Capítulo es adecuado para principiantes absolutos, y de afianzamiento a los ya iniciados. 1.2. COMENZANDO Al arrancar MatLab, presenta una pantalla dividida en varias ventanas, configurables desde Desktop Layout del menu de View; en una de las ventanas estará el cursor parpadeando a la derecha de «>>», es la ventana de comandos desde donde se ejecutan los mismos, las otras son informativas: >> 3 + 4 ans = 7 >> 3*5 ans = 15 >> 15/3 ans = 5 1

2


>> 15\3 ans = 0.2000 >> 2ˆ3 ans = 8 >> sin(2*pi*30/360) ans = 0.5000

Figura 1.1. Ventanas de MatLab.

1.3. ESPACIO DE TRABAJO Los datos y variables usadas residen en el espacio de trabajo, workspace, accesible desde la ventana de comandos. Para inspeccionar el contenido de este espacio se utilizan los comandos who y whos. Los archivos directamente accesibles desde el espacio de trabajo, se muestran mediante what. En el siguiente ejemplo se muestran sus características: >> t = linspace(1,10,4) % Crea un vector de 4 elementos desde % 1 a 10. t = 1 4 7 10

PRIMEROS PASOS EN MATLAB

>> t = t(:) t = 1 4 7 10

3

% Crear el vector t en columna.

>> A = 2*t; B = 2; >> who Your variables are: t

A

B

>> whos Name t A B

Size 1 × 4 1 × 4 1 × l

Bytes 32 32 8

Class double array double array double array

Grand total is 9 elements using 72 bytes >> what M-files in the current directory D:\MatLab\work AjusNL Datos EcudifP AjusN1Fun Ecudif Fun MAT-files in the current directory D:\MatLab\work Datos MDL-files in the current directory D:\MatLab\work Bcont Bfuntab BnoLineal Bseñsis Bfuentes Bmat Bsalidas Ecu2

En el espacio de trabajo se crearon 3 variables, 2 de 4 elemenos, y una de 1 elemento, de modo que son 9 elementos a 8 bytes por elemento, lo que hace un total de 72 bytes. A partir de la versión 6, release 12, MatLab incorporó workspace, potenciando la capacidad de whos. 1.4. VARIABLES En MatLab no es necesario hacer declaraciones previas acerca de las variables. El contenido de las variables de caracteres ha de ir delimitado por el signo «'». >> numero_de_visitantes = 25 numero_de_visitantes = 25

4


De ese modo se crean variables numéricas, numero_de_visitantes, que almacenan su valores en una matriz, en este caso la matriz es de 1 × 1, y su valor es 25. » size(numero_de_visitantes) % Dimensión de variables. ans = 1

1

» Nombre='Pepe'; » size (Nombre)

% Variable de caracteres.

ans = 1

4

Los nombres de las variables deben seguir estas reglas: 1. Se forman con las letras del abecedario, los dígitos 0 a 9 y el signo «_», distinguiéndose mayúsculas de minúsculas. 2. Los nombres de las variables han de comenzar por una letra y no deben contener espacios en blanco. 3. Los nombres de las variables no pueden coincidir con los nombres de las keywords, nombres reservados. La lista de los nombres reservados se obtiene por medio de iskeyword: » iskeyword ans = 'break' 'case' 'catch' 'continue' 'else' 'elseif' 'end' 'for' 'function' 'global' 'if' 'otherwise' 'persistent' 'return' 'swirch' 'try' 'while'


5

>> if = 5 ??? if = 5 Error: Expected a variable, function, or constant, found”=”.

Los nombres de las variables pueden ser tan extensos como se quiera, pero MatLab sólo reconoce los 31 primeros caracteres. Las variables se eliminan del espacio de trabajo con el comando clear: clear clear variables clear global clear functions clear all clear pipo*

Elimina las variables del espacio de trabajo. Es equivalente al comando anterior. Elimina las variables globales. Elimina todas las funciones compiladas. Elimina todas las variables, globales y funciones. Elimina las variables que empiezan por pipo.

MatLab suministra amplia información adicional mediante help clear. 1.5. FORMATO DE NÚMEROS MatLab presenta los resultados numéricos en varios formatos, según se expresa a continuación: >> help format FORMAT Set output format. All computations in MATLAB are done in double precision. FORMAT may be used to switch between different output display formats as follows: FORMAT Default. Same as SHORT. FORMAT SHORT Scaled fixed point format with 5 digits. FORMAT LONG Scaled fixed point format with 15 digits. FORMAT SHORT E Floating point format with 5 digits. FORMAT LONG E Floating point format with 15 digits. FORMAT SHORT G Best of fixed or floating point format with 5 digits. FORMAT LONG G Best of fixed or floating point format with 15 digits. FORMAT HEX Hexadecimal format. FORMAT + The symbols +, – and blank are printed for positive, negative and zero elements. Imaginary parts are ignored. FORMAT BANK Fixed format for dollars and cents. FORMAT RAT Approximation by ratio of small integers. Spacing: FORMAT COMPACT Suppress extra line-feeds. FORMAT LOOSE Puts the extra line-feeds back in.

6


Al mostrar resultados numéricos, MatLab sigue estas dos reglas: 1. MatLab intenta mostrar números enteros. Si el entero es muy grande, se presenta en formato exponencial, con 5 cifras significativas. 2. Los números con decimales se muestran con 4 o 5 cifras significativas. Los números en valor absoluto menores de 0,01 y mayores de 1.000, se muestran en formato exponencial. A continuación se muestran ejemplos demostrativos de formatos numéricos: >> sqrt(2) ans = 1.4142 >> format long >> sqrt(2) ans = 1.41421356237310 >> format Long e >> sqrt (2) ans = 1.41421356237310e+00 >> format short >> A= [10000 0.0001] ans = 1. 0e+04 * 1.0000 0.0000 >> format short g >> A A = 10000 >> format rat >> A

0.0001

A = 10000 1/10000


7

1.6. PROGRAMAS MatLab acepta comandos directos, para inmediatamente producir el resultado o ejecutar una serie de comandos almacenados en un archivo, con la extensión «.m». Un archivo.m, consiste en una secuencia de sentencias MatLab, posiblemente incluyendo referencias a otros archivo.m, o recursivamente a sí mismo. A estos archivos los llamamos programas MatLab, en inglés scripts. Las variables de los programas se mantienen en el espacio de trabajo, pudiendo ser invocadas en cualquier momento para ver su contenido. En una sentencia, lo que sigue a % no se ejecuta, se considera un comentario. Si se desea construir una tabla con inversos, cuadrados y raíces cuadradas de 1 a 10, se edita un archivo, Numeros.m, con cualquier editor, tal como el bloc de notas del sistema operativo, o con el editor propio de MatLab, según: % ----------------Numeros.m ----------------------------------x=1:10; % Crea un vector de 1 a 10 de 1 en 1. Vector en línea. x=x'; % Transposición. Vector en columna. x=[x,1../x,x.ˆ2,sqrt(x)]; % Matriz de 4 columnas. % ------------------------------------------------------------

El programa se invoca ejecutando «Numeros». Como en el programa todas las sentencias se finalizaron con «;», no se muestra ningún valor numérico. Al ejecutar «x», se obtendrá la tabla desada: >> Numeros >> x x = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000

1.0000 1.0000 0.5000 4.0000 0.3333 9.0000 0.2500 16.0000 0.2000 25.0000 0.1667 36.0000 0.1429 49.0000 0.1250 64.0000 0.1111 81.0000 0.1000 100.0000

1.0000 1.4142 1.7321 2.0000 2.2361 2.4495 2.6458 2.8284 3.0000 3.1623

Pulsando la tecla ↑, se consiguen las líneas de los comandos previamente ejecutadas.

8


1.7. FUNCIONES El otro tipo de archivos utilizado por MatLab son las funciones, cuya primera característica es que sus variables son locales en su entorno y no definidas en el espacio de trabajo, ni en otras funciones. Buena parte de la potencia de MatLab se basa en su extenso conjunto de funciones, las básicas y las distribuidas de forma separada para aplicaciones específicas, MatLab toolboxes, y otras que desarrollan los usuarios. Las funciones toman unas variables de entrada para calcular unos datos de salida, sea: Fun( x ) =

1 1 + −5 2 ( x − 1) + 0,1 ( x − 3)2 + 0, 2

almacenado en el archivo «Fun.m», cuyo contenido es: function y=Fun(x) % -----------------Fun.m -------------% Esto es un ejemplo % de una función. y=1../((x–1).ˆ2+0.1)+1../((x–3).ˆ2+0.2)–5; %-------------------------------------

Para evaluar Fun gráficamente, se lanza con las siguientes instrucciones: >> x=–2:0.01:6; % Vector de –2 a 6, a incrementos de 0,01. >> y=Fun(x); % Guardando el vector Fun(x) en y. >> plot(x,y),grid % Representación con rejilla.

En Fun(x), x es el argumento o entrada de Fun, para dar unos resultados de salida que se almacenan en y, que se muestran gráficamente. Hay funciones del sistema o construidas por un usuario, que toman diferente número de argumentos de entrada que de salida. Así la función max, toma un vector de argumentos y puede suministrar una o dos salidas, según se use: >> A = [1 2 1 5 2 3]; >> max(A) % Suministrará el valor máximo de A. ans = 5 >> [X, i] = max(A) % X, valor máximo, i posición del máximo. X = 5 i = 4


9

4

2

0

–2

–4

–6 –2

–1

0

1

2

3

4

5

6

Figura 1.2. Representación gráfica de Fun.

1.7.1. Reglas de construcción de funciones 1. El nombre de la función y del archivo deben ser idénticos. 2. Los nombres de las funciones se rigen por las normas de los nombres de las variables. 3. La primera línea ejecutable de una función debe ser la línea de declaración de función. 4. Las variables del interior de las funciones son variables locales. 5. El conjunto de líneas consecutivas de comentarios que siguen a function, componen el texto de ayuda de esa función, obtenible mediante help y lookfor. >> help Fun ----------------Fun.m ------------Esto es un ejemplo de una función.

6. Una función termina al encontrar un retun o al llegar a la última línea de la función. 7. Si una función llama a un programa, éste es evaluado en el espacio de trabajo de la función y no en el workspace de MatLab. 8. Cada función tiene su espacio de trabajo separado del de MatLab, de modo que la conexión entre estos ambientes se realiza a través de las variables de entrada y salida de la función.

10


9. Para compartir variables del interior de las funciones y del espacio de trabajo, se declaran variables globales donde se necesiten mediante la instrucción global. >> global X a b

% Declaración de variables globales

Para facilitar el manejo de funciones, MatLab incorporó recientemente @, y feval, para mejorar eval, cuya utilidad se expone en el siguiente ejemplo: >> F = @Fun

% Creación directa de F.

F = @Fun >> feval(F,2) ans = -3.2576

Se consigue el mismo resultado con: >> eval('Fun(2)') ans = -3.2576 >> Fun(2) ans = -3.2576 >> F(2)=@cos

% Creación directa de F(2).

F = @Fun

@cos

La eficiencia de feval es considerablemente superior a eval, ya que el primero evalúa directamente lo que se invoca, mientras que eval llama al interpretador completo de MatLab. La diferencia en tiempo de ejecución de ambas funciones se pone de manifiesto con: >> tic, for i = 1:100000, a = eval('Fun(i)'); end, toc elapsed_time = 14.3210 >> tic, for i = 1:100000, a = feval('Fun',i); end, toc elapsed_time = 4.0960


11

1.7.2. Funciones en línea Un segundo modo de definir funciones, sin editar archivos, se logra con inline: >> syms x y >> f = inline('x.ˆ2 + y.ˆ2') f = Inline function: f(x) = x.ˆ2 + y.ˆ2 >> f(3,4) ans = 25 >> feval(f,3,4) ans = 25

1.7.3. Ejemplo de función recursiva En muchas aplicaciones se presenta la recursividad, función que en su interior se llama a sí misma. El ejemplo más secillo de recursividad es el cálculo del factorial de un número: N! = N × (N – 1)! definiendo 1! como 1. La función Factorial.m, toma un número de entrada, y suministra como salida su factorial: function f = Factorial(N) % Esta función calcula el factorial de la parte % entera de un número. n = fix(N); % n toma la parte entera de N. if n > 1 f = n*Factorial(n – 1); else f = 1; end

Esta sencilla función toma un número y calcula el factorial de su parte entera. >> Factorial(3.3) ans = 6

12


1.8. NÚMEROS COMPLEJOS MatLab admite operaciones con números complejos, permitiendo usar indistintamente la i y la j, según se muestra en lo siguiente: >> a = sqrt(–1) a = 0 + 1. 0000i >> conj(a) ans = 0 – 1.0000i >> sqrt(a) ans = 0.7071

+0.7071i

>> exp(2i) ans = -0.4161 + 0.9093i >> A = (3 + 4i)*(2 – j) A = 10. 0000 + 5.0000i >> r = real(A) r = 10 >> I = imag(A) I = 5 >> r = abs(A) r = 11.1803 >> Angulo = angle(A)


13

Angulo = 0.4636 >> Angulo = atan2(imag(A),real(A)) 0.4636 >> Aa = r*exp(Angulo*i) Aa = 10.0000 + 5.0000i

1.9. MANEJO DE VECTORES Y MATRICES La forma más sencilla de crear un vector es mediante el uso de [], vector en línea, o con []', vector en columna. Los elementos se separan por espacios o comas, el «;» se reserva para anexar en columna: >> t = [3 5 7, 8, 9] t = 3 5 7 8 9

También se generan vectores mediante las instrucciones linspace y logspace, ambos con dos o tres argumentos, y con «:», ya mencionado: >> x = logspace(0,2,5) % Vector de 5 componentes de 10ˆ0 a 10ˆ2. x = 1.000

3.1623

10.0000

31.6228

100.000

Con el siguiente ejemplo se muestra la creación y manejo de matrices: >> x = 0:4 x = 0

1

2

3

4

>> y = x.ˆ2

% El punto antes del exponente % hace que la exponenciación % sea elemento a elemento.

y = 0

1

4

>> a = [x;y]

9

16 % Crear una matriz anexando % vectores.

14


a = 0 0

1 1

2 4

3 9

4 16

>> A = a'

% Crear la matriz A, transpues% ta de a.

A = 0 1 2 3 4

0 1 4 9 16

» B = [A; 5 25]

% Añadir una línea a una ma% triz.

B = 0 1 2 3 4 5

0 1 4 9 16 25

>> C = reshape(B,3,4)

% Reconfigurar la matriz B % con 3 líneas y 4 columnas.

C = 0 1 2

3 4 25

0 1 4

9 16 50

>> C(2,:)=[]

% % % %

C = 0 2

3 25

0 4

9 50

>> C(:,3) = [] C = 0 2

3 25

9 50

>> A = [1 2; 3 5] A = 1 3

2 5

Eliminar la 2a línea. Los dos puntos indican para todos los valores de esa dimensión.

% Eliminar la 3a columna, % 0 4.


>> Aˆ2

15

% Diferencia entre ˆ y .ˆ

ans = 7 18

12 31

>> A.ˆ2 ans = 1 9

4 25

>> 1../A ans = 1.0000 0.3333

0.5000 0.2000

>> det(A) ans = -1 >> inv(A) ans = -5.0000 3.0000

2.0000 -1.0000

Las exponenciaciones, elemento a elemento y matricial, también se realizan con los comandos power y mpower, según: >> x =[2 3 4]; >> y = power(x,2)

% Equivalente a y = x.ˆ2

y = 4

9

16

>> x = [2 3;1 4]; >> y = mpower(x,2) y 7 6

18 19

% Equivalente a y = xˆ2

16


1. 10. POLINOMIOS Sea el polinomio p = x2 – 5x + 6, con vector de coeficientes C: >> C = [1 -5 6];

Las raíces de este polinomio se obtienen, mediante: >> r = roots(C) r = 3 2

poly sobre las raíces devuelve los coeficientes del polinomio: >> poly(r) ans = 1

-5

6

polyval evalúa el polinomio sobre un valor: >> polyval(C,5) ans = 6

1.10.1. Multiplicación y división de polinomios La multiplicación de dos polinomios, Pol1 = x2 + 2x + 3 por Pol2 = x – 1, se efectúa por medio de conv: >> Pol1 = [1 2 3]; >> Pol2 = [1 -1]; >> PolProd = conv(Pol1,Pol2) PolProd = 1

1

1

-3

Cuyo resultado equivale a x3 + x2 + x – 3. La división de polinomios se realiza mediante deconv: >> PolDiv = deconv(PolProd,Pol2) PolDiv = 1

2

3


17

1.10.2. Desarrollo en fracciones simples Mediante residue, aplicado a dos polinomios, se obtienen fracciones simples, cuya suma es equivalente al cociente de los polinomios: >> P1 = [5 -20 9]; >> P2 = [1 -5 4]; >> [r,s,t] = residue(Pl,P2) r = 3 2 s = 4 1 t = 5

Representando r, un vector columna con los numeradores de las fracciones, s las raíces de cada denominador, y t los coeficientes del término independiente: 5 x2 − 20 x + 9 3 2 = + +5 x − 4 x −1 x2 − 5 x + 4

Si se conocen los numeradores, las raíces de los denominadores y los términos independientes, se pueden generar los polinomios: >> [p1,p2] = residue(r,s,t) p1 = 5

-20

9

1

-5

4

p2 =

Siendo p1 y p2, los polinomios obtenidos, con los numeradores, raíces y términos independientes. 1.10.3. Derivadas de polinomios Las derivadas de los polinomios se obtienen con polyder: >> derp1 = polyder(p1) derp1 = 10 - 20

es decir, 10x – 20

18


1.10.4. Integración de polinomios La integral de un polinomio se obtiene mediante polyint(p), o con polyint(p,C); en el primer caso se supone que la constante de integración es 0, y en el segundo C: >> polyint(p2) ans = 0.3333

-2.5000

4.0000

0

4.0000

2.0000

>> polyint(p2,2) ans = 0.3333

-2.5000

Equivalente a: 1 3 5 3 x − x + 4x + 2 3 2

1.10.5. Interpolación polinomial Matlab permite varios modos de interpolación, relacionados y descritos en help interp1; a continuación se muestran unos ejemplos: >> >> do >>

x = [2 4 6]; y = power(x,2); de y5 = interp1(x,y,5)

% Por defecto se utiliza el méto% interpolación lineal.

y5 = 26 >> y5c = interp1(x,y,5,'cubic') y5c = 24.8750 >> y5s = interp1(x,y,5,'spline') y5s = 25

% Interpolación cúbica.


19

1.11. matlabpath El path de MatLab, llamado matlabpath, establece el camino para buscar variables, programas y funciones de MatLab que sean llamados directamente desde la ventana de comandos, o durante la ejecución de programas. Al invocar un comando, MatLab lo busca, y ejecuta el que primero coincida con el nombre invocado, según el siguiente orden: 1. 2. 3. 4.

Si es una variable del workspace. Si es una función incorporada. Si es un archivo.m presente en el directorio actual. Lo busca siguiendo el orden establecido en matlabpath.

El comando which, aplicado sobre una función muestra su ruta: >> which Fun d:/MatLab6p5/work/Fun.m

1. 12. lookfor lookfor aplicado a una variable de caracteres, busca esos caracteres en la primera línea de comentarios de los archivo.m encontrados en el matlabpath. >> lookfor Fun.m Fun.m: %---------------- Fun.m --------------------

1.13. LATEX La función latex(A) devuelve la representación LATEX de una expresión simbólica: >> syms x >> A = taylor(exp(–x)) A = 1 - x + 1/2*xˆ2 - 1/6*xˆ3 + 1/24*xˆ4 - 1/120*xˆ5 >> pretty(A) 2 3 4 5 1 - x + 1/2 x - 1/6 x + 1/24 x - 1/120 x >> latex(A) ans = 1 - x + 1/2\,{x}ˆ{2}-1/6\,{x}ˆ{3}+1/24\,{x}ˆ{4}{\frac{1}{120}}\, {X}^{5}

20


1.14. FUNCIONES DEL TIEMPO MatLab suministra datos relacionados con el tiempo, como una variable de caracteres, 12-Oct-1498, como una variable numérica, 739476 o como un vector, 1789 07 14 0 0 0. La función clock responde con: >> Tiempo = clock Tiempo = 1.0*e + 03 * 2.0030 0.0010

0.0290

0.0160

0.0400

0.0502

Esos datos responden a Tiempo[año mes día hora minutos segundos]. now devuelve la fecha y hora en un número, datestr convierte una fecha numérica en string, datevec individualiza los componentes de fecha y hora, datenum convierte una fecha en un número: >> now ans = 7.3161e + 05 >> datestr(7.3161e + 05) 29-Jan-2003 >> datevec(7.3161e + 05) 2003

1

29

0

0

0

>> datenum(date) 7.3161e + 05

La función date devuelve una variable de caracteres en el formato ddmmm-yyyy: >> date ans = 29-Jan-2003

Numerosas transformaciones de fechas pueden encontrarse mediante help datestr.


21

Para temporizaciones se usan tic, para comienzo, y toc para finalizar la temporización y mostrar el resultado: >> T = clock; tic, for i = 1:1000000, a = sqrt(i); end, ... toc, Tt = etime(clock,T) elapsed_time = 3.6174 Tt = 3.6642

Los tres puntos seguidos indican continuación en la siguiente línea. La ejecución de las raíces cuadradas de 1 a 1 millón, en un Pentium III a 866 Mhz, con Linux Suse 8.2 y Student MatLab 12, tarda 3.62 segundos. La función etime devuelve el tiempo, en segundos, transcurrido entre dos valores del tiempo: >> etime(Tiempo,clock) ans = -1.2835e + 03

1.15. INTERCAMBIO DE DATOS A partir de la versión 6, taınbién numerada como 12, MatLab presenta una ventana de Historia de Comandos, en donde se listan las órdenes efectuadas, de modo que en una sesión nueva se pueden buscar comandos ejecutados en sesiones precedentes. Para guardar el espacio de trabajo en un archivo, para posteriormente recuperarlo, se utilizan los comandos load y save, según: >> clear >> X=rand(2,3) X = 0.3046 0.1897

0.1934 0.6822

>> Y=round(X) Y = 0 0

0 1

0 1

0.3028 0.5417

22

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


save mi_archivo clear who load mi_archivo who

Your variables are: X Y >> X X = 0.3046 0.1897

0.1934 0.6822

0.3028 0.5417

>> Y Y = 0 0

0 1

0 1

round se relaciona con ceil y floor, cuyas funciones se intuyen por su significado: >> A=[1.5 2.49; 0.2 9.99] A = 1.5 0.2

2.49 9.99

>> ceil(A) ans = 2 3 1 10 >> floor(A) ans = 1 0

2 9

save y load permiten salvar o cargar variables determinadas, separadas por espacios, a la vez que admiten el formato ASCII, añadiendo al final de estos comandos -ascii, como se especifica en el siguiente ejemplo: >> save exp_mayo.dat X -ascii


23

En el archivo exp_mayo.dat, se guarda la variable X en formato ASCII. Para guardar la sesión de trabajo, se utiliza el comando diary, que almacena una copia de todas las entradas realizadas desde el teclaclo en un archivo, en el directorio actual, en formato ASCII:

2 Cálculo simbólico

2.1. INTRODUCCIÓN MatLab se caracterizó desde un principio por ser muy potente en cálculo numérico, mientras que el cálculo simbólico fue incorporado como una toolbox, cuando MathWorks, empresa que comercializa MatLab, se extendió internacionalmente. La Symbolic Math Toolbox es una colección de herramientas para MatLab, que se utilizan para manejar y resolver expresiones simbólicas. Las herramientas simbólicas disponibles más usadas son; combinar, simplificar, factorizar, derivar, integrar, límites, resolución de sistemas de ecuaciones algebraicas o diferenciales, transformaciones integrales, la mayoría de las operaciones del álgebra lineal... Estas herramientas de cálculo simbólico son parte del programa Maple V, comercializado por Waterloo Maple Software Inc. 2.2. OBJETOS Y EXPRESIONES SIMBÓLICAS En Matlab hay dos tipos de objetos, numéricos y literales, strings. La Symbolic Math Toolbox usa objetos simbólicos para representar variables y operadores, por ejemplo: >> x=sym('x') x = x

25

26


Un objeto numérico puede convertirse en simbólico según: >> M=magic(2)

% M, variable numérica.

M = 1 4

3 2

>> N=sym(M)

% N, variable simbólica.

N = [1, 3] [4, 2] >> syms a b c d >> Mat=[a,b;c,d]

% Definición de simbólicos.

Mat = [a, b] [c, d] >> det(Mat)

% Cálculo del determinante.

ans = a*d-b*c >> M=(a-b)/(c+d) M = (a-b)/(c+d) >> pretty(M) a - b c + d

Para la simplificación y transformaciones de expresiones se utilizan los operadores collect, expand, horner, factor, simple y simplify, algunos de los cuales se aplican en lo que sigue. 1/2

2 15

1/2

/5 sin(15

/2 t) exp(–3/2 t)

80 70 60 50

40 30 20 10 0 –10 –3

–2

–1

0 t

1

2

3

Figura 2.1. Representación gráfica de y.

CÁLCULO SIMBÓLICO

27

2.3. EJEMPLOS DE CÁLCULO SIMBÓLICO 2.3.1. Derivadas e integrales Se crea la función y, dependiente de t, según: >> t=sym('t'); >> y=sym('2*15^(1/2)/5*sin(15^(1/2)/2*t)*e^(-3/2*t)'); >> pretty (y)

2/5 15

1/2 1/2 (- 3/2 t) sin(1/2 15 t) e

La representación gráfica de y, para valores de -π a π, se consigue con la instrucción: >> ezplot(y,[-pi pi])

La derivada y’, almacenada en dy, simplificada y factorizada se obtiene según: >> y Y= 2*15^(1/2)/5*sin(15^(1/2)/2*t)*exp(-3/2*t) >> >> >> >>

dy=diff(y); dy=simplify(dy); dy=factor(dy); pretty(dy)

1/2 1/2 1/2 - 3/5 exp(-3/2 t) (-5 cos(1/2 15 t)+15 sin(1/2 15 t))

La representación gráfica de dy se obtiene, como en el caso anterior, mediante ezplot. La integral de la derivada es la propia función, por lo tanto, integrando dy, ha de obtenerse una expresión idéntica a la de y: >> ezplot(dy, [-pi pi]) >> Intdy=int(dy,t); >> Intdy=simplify(Intdy) Intdy = 2/5*15^(1/2)*sin(1/2*15^(1/2)*t)*exp(-3/2*t)

28


1/2

3/5 exp(3/2 t) (- 5 cos(1/2 15

1/2

sin(1/2 15

1

2

t)+15

1/2

t))

60 40 20 0 –20 –40 –60 –80 –100 –3

–2

–1

0 t

3

Figura 2.2. Representación gráfica de dy.

2.3.2. Sistemas de ecuaciones Los sistemas de ecuaciones se resuelven mediante la instrucción solve, tomando como argumentos el primer miembro de las ecuaciones igualadas a 0, según: x+ y= 5 x – y = –1 ez·x = 7,389 >> syms x y z >> [x,y,z]=solve(x+y-5,x-y+1,exp(x*z)-7.389) x = 2 y = 3 z = 1/2*log(7389/1000) >> z=double(z) z = 1.0000

Si las ecuaciones a resolver no tienen solución analítica, solve devuelve una expresión numérica, que se ejecuta mediante double, transformación a doble precisión.

CÁLCULO SIMBÓLICO

29

2.3.3. Ecuaciones diferenciales La función dsolve calcula las soluciones simbólicas de ecuaciones diferenciales ordinarias. Los argumentos de dsolve deben ser expresiones de caracteres, strings, conteniendo el signo «=». Para indicar la derivada primera se utiliza el signo D. para la derivada segunda se utiliza D2, y así sucesivamente. La sintaxis de esta operación se expresa como: r=dsolve('Ecu1,Ecu2,...','Cond1,Cond2,...','x')

siendo x la variable independiente; si no se expresa se utiliza t por defecto. Obteniendo y'', por derivación de y', se puede componer con y e y', la siguiente ecuación diferencial: d2 y dy +3 + 6y = 0 2 dt dt

La solución a esta ecuación diferencial se obtiene según: >> yc=dsolve('D2y+3*Dy+6*y=0,Dy(0)=3,y(0)=0') yc = 2/5*15^(1/2)*exp(-3/2*t)*sin(1/2*15^(1/2)*t)

La especificación de las constantes iniciales es opcional, a continuación se muestran ejemplos con y sin estas constantes: >> r=dsolve('D2y+3*Dy+2*y=0') r = C1*exp(-t)+C2*exp(-2*t) >> s=dsolve('D2y+3*Dy+2*y=0','Dy(0)=l,y(0)=0') s = exp(-t)-exp(-2*t) >> S=dsolve('D2y+3*Dy+2*y=cos(t)','Dy(0)=l,y(0)=0') S = 1/l0*cos(t)+3/l0*sin(t)+1/2*exp(t)-3/5*exp(-2*t)

30


2.4. TRANSFORMACIÓN DE LAPLACE E INVERSA La transformación de Laplace calcula la integral: ∞

(s) = ∫0 f (t) ⋅ e – s⋅t dt

permite transformar f(t) en el dominio del tiempo, a (s), en el dominio de la variable compleja >> syms a w t s >> F=cos(w*t)*exp(–a*t) F = cos(w*t)*exp(–a*t) >> L=laplace(F,t,s) L = (s+a)/((s+a)^2+w^2) >> pretty(L) (s + a) ----------2 2 (s + a) + w >> L1=laplace(exp(-t),t,s) L1 = 1/(l+s) >> pretty(L1) 1 ---1 + s

La transformada inversa se ejecuta sobre expresiones de la variable compleja s, para volver al dominio del tiempo. >> I=ilaplace(L,s,t) I = cos(w*t)*exp(-a*t) >> I1=ilaplace(Ll,s,t) I1 = exp(-t)

CÁLCULO SIMBÓLICO

31

2.5. LÍMITES >> limit(1/x,x,0,'left') ans = -inf

El límite de 1/x, cuando la variable x tiende a 0 por la izquierda, es –∞. >> Lim=limit((x^3-1)/(x^2-l),x,1) Lim = 3/2

2.6. SERIES DE TAYLOR Y MAC LAURIN Para obtener desarrollos de series de Taylor y Mac Laurin, se utiliza indistintamente taylor, con 3 o 4 argumentos: >> syms x >> f=exp(-x); >> pretty(taylor(f,x,3,8)) 2 exp(-8) - exp(-8) (x - 8) + 1/2 exp(-8) (x - 8) >> pretty(taylor(f,x,3)) 2 1 - x + 1/2 x

2.7. INVOCANDO A MAPLE V Desde MatLab se accede a Maple V con la función maple, que toma como argumento la expresión con sintaxis de Maple V, mediante las instrucciones genéricas: r=maple('Sentencia Maple') r=maple('Función',Argl,Arg2,...)

2.7.1. Transformación de Laplace e inversa con Maple V Para comenzar a utilizar las funciones de MapleV, relacionadas con la transformación de Laplace, es necesario cargar previamente el paquete de transformaciones integrales with(inttrans): >> maple('with(inttrans)') % Se carga el paquete de ... % transformaciones integrales.

32


ans = [addtable, fourier, fouriercos, fouriersin, hankel, hilbert,... invfourier, invhilbert, invlaplace, invmellin, ... laplace, mellin, savetable] >> maple('laplace(y(t)=t^2+sin(t),t,s)') ans = laplace(y(t),t,s) = 2/s^3+1/(s^2+1) >> M1=maple('invlaplace((s+1)/(s*(s^2+s+1)),s,t)') M1 = 1+1/3*exp(-1/2*t)*3^(1/2)*sin(1/2*3^(1/2)*t)-exp(-1/2*t)*... cos(1/2*3^(1/2)*t)

2.7.2. Resolución de ecuaciones diferenciales con Maple V Para resolver la ecuación diferencial: ÿ + 3 y· + 6y = 0 con condiciones iniciales: y· (0) = 3 y (0) = 0 se procede con las siguientes instrucciones: >> maple('eq:=diff(y(t),t$2)+3*diff(y(t),t)+6*y(t)=0') ans = eq := diff(y(t),'$'(t,2))+3*diff(y(t),t)+6*y(t) = 0 >> maple('ini:=y(0)=0,D(y)(0)=3') ans = ini := y(0) = 0, D(y)(0) = 3 >> maple('Sol:=dsolve({eq,ini},{y(t)})') ans = Sol := y(t) = 2/5*15^(1/2)*exp(-3/2*t)*sin(1/2*15^(1/2)*t)

CÁLCULO SIMBÓLICO

2.7.3. Resolución de ecuaciones diferenciales con la transformación de Laplace Para resolver la ecuación diferencial: ÿ – 2 y· – 3y = 0 con condiciones iniciales: y· (0) = 1 y (0) = 0 se procede con las siguientes instrucciones: >> maple('eq:=diff(y(t),t$2)-2*diff(y(t),t)-3*y(t)=0') ans = eq := diff(y(t),'$'(t,2))-2*diff(y(t),t)-3*y(t) = 0 >> maple('ini:=y(0)=0,D(y)(0)=1') ans = ini := y(0) = 0, D(y)(0) = 1 >> maple('La:=laplace(eq,t,s)') ans = La := s*(s*laplace(y(t),t,s)-y(0))-D(y)(0)-2*s* ... laplace(y(t),t,s)+2*y(0)-3*laplace(y(t),t,s) = 0 >> Maple('Sol:=subs(ini,{La})') ans = Sol := {s^2*laplace(y(t),t,s)-1-2*s*1aplace(y(t),t,s) ... -3*laplace(y(t),t,s) = 0} >> maple('Sol1:=solve(Sol,{laplace(y(t),t,s)})') ans = Sol1 := {laplace(y(t),t,s) = 1/(s^2-2*s-3)} >> maple('Solf:=invlaplace(Sol1,s,t)') Solf =

33

34


{y(t) = 1/16*16^(1/2)*(exp((1+1/2*16^(1/2))*t)- ... exp((1-1/2*16^(1/2))*t))} >> maple('simplify(Solf)') ans = {y(t) = -1/4*(-1+exp(-4*t))*exp(3*t)}

1 y (t ) = – (e –4 t – 1) e 3t 4

3 Sentencias de control de flujo

3.1. input La forma de input, se indica en los ejemplos siguientes, según se trate de variables literales o numéricas: >> R=input('> Cuál es tu nombre ? ','s') > Cual es tu nombre ? Pepe >> ['Hola',R] ans = Hola Pepe >> ['Hola';R] Hola Pepe >> P=input('> Dime el radio de la rueda ? ') > Dime el radio de la rueda ? 5 >> P+1 ans = 6

35

36


3.2. if - else - end El salto condicional de flujo más sencillo, se construye de la siguiente manera: if condición comandos end

Los comandos entre if y end, se ejecutarán, si la condición es verdad, uno. Si la condición es falsa, cero, se puentearán los citados comandos. Operadores relacionales eq ne lt gt le ge

== ~= < > <= >=

Igual No igual Menor que Mayor que Menor que o igual Mayor que o igual

En la primera de las siguientes instrucciones se responde con «hola», ya que es verdad que 'a' es igual a 'a'. En la segunda esto se obvia, dado que la condición es falsa. >> if 'a'=='a', 'hola', end ans = hola >> if 'a'=='b', 'hola', end

Para evaluar dos alternativas, la construcción if - else - end toma la siguiente construcción: if condición Ejecución de comandos de esta zona, si la condición es verdad else Ejecución de comandos de esta zona, si la condición es falsa end

Si hubiese tres alternativas, la construcción sería:

SENTENCIAS DE CONTROL DE FLUJO

37

if condición 1 Ejecución de instrucciones si la condición 1 es verdad elseif condición 2 Ejecución de instrucciones si la condición 2 es verdad elseif condición 3 Ejecución de instrucciones si la condición 3 es verdad else Ejecución de instrucciones si ninguna condición es verdad end

En un comercio se vende un vino a un precio condicionado por la cantidad requerida. Hasta 5 botellas el precio unitario es de 6 €, desde 6 a 12 botellas el precio es de 5,5 €, y a partir de 13, a 5 € la botella. Elaborar un programa, que pregunte cuántas botellas se desean e indique el precio unitario y el total del gasto. %------------ Vino.m -------------C=input('Cuantas Botellas ?'); if C<=5 Pu=6; Pt=Pu*C; elseif C<=12 Pu=5.5; Pt=Pu*C; else Pu=5; Pt=Pu*C; end Pu Pt % ----------------------------------

A continuación se muestra la ejecución de este programa: >> Vino Cuantas Botellas ? 8 Pu = 5.5000 Pt = 44

38


3.3. while - end La sentencia while - end funciona según: while

Condición

Comandos

end

Se ejecutarán los comandos de este bucle, mientras la Condición sea verdad. Se expone a continuación la generación de una tabla que suministre los inversos, cuadrados y raíces cuadradas del 1 al 3: % ------------------ Tabla.m -----------------I=0; while I<4 I=I+1; J(I)=I; A(I)=1/I; B(I)=power(I,2); C(I)=sqrt(I); end D=[I;J;A;B;C]; D=reskape(D,3,4); % -------------------------------------------->> Tabla >> D D = 1.0000 2.0000 3.0000

1.0000 0.5000 0.3333

1.0000 4.0000 9.0000

1.0000 1.4142 1.7321

3.4. for - end El bucle o lazo de control for, permite realizar un conjunto de instrucciones iguales, variando uno o varios subíndices, su forma genérica es: for I = vector Comandos end

A continuación se muestra la generación de la tabla anterior utilizando for.

SENTENCIAS DE CONTROL DE FLUJO

39

% --------------- Tabla1. m ----------for I=1:3 J(I)=I; A(I)=1/I; B(I)=power(I,2); C(I)=sqrt(I); end D=[I;J;A;B;C]; D=reshape (D,3,4); % ------------------------------------>> Tabla1 >> D D = 1.0000 2.0000 3.0000

1.0000 0.5000 0.3333

1.0000 4.0000 9.0000

1.0000 1.4142 1.7321

El funcionamiento de for, en el programa anterior, comienza tomando I el valor de 1, recorriendo los comandos del bucle de arriba hacia abajo, creándose J(1) con el valor 1, A(1) con el valor del inverso de 1... Al llegar al final, I se incrementa en 1, tomando el valor 2, y el flujo de cálculo vuelve a cabecera, donde se comprueba si I rebasa o no el valor límite, en este caso 3, de modo que continuará el ciclo, se creará J(2) con el valor 2, A(2) con el valor de 1/2 y así sucesivamente. Cuando I alcance el valor de 4 ya no se ejecutarán las sentencias del interior del for, se continuará con lo que venga después de end. A continuación se muestra una estructura de un for dentro de otro, anidamiento, para ejercitar el funcionamiento de los subíndices. Para cada variación de un valor del índice del for externo, se ejecutan todas las variaciones del interno. El funcionamiento del siguiente ejemplo comienza con I=1 y J, tomando los valores de 1, 2, 3 y 4, entonces I pasa al valor 2, y J vuelve a tomar los valores de 1, 2, 3 y 4, y así hasta completar la última vuelta con I=3. % ------------- Tabla2.m ---------------------A=[]; % Crear una matriz vacia. x=[1 -1 2 0.5]; for I=1:3 for J=1:4 A(I,J)=I^x(J); end end % ---------------------------------------------

En la tabla formada, la I varía de 1 a 3 al descender por cada columna, y la J varía en horizontal, del principio al final de cada línea, tomando siempre los

40


valores de 1 para la primera columna, 2 para la segunda, 3 para la tercera y 4 para la cuarta y así en cada línea. En la primera línea del programa Tabla2 se crea la matriz A sin contenido, para seguridad, siendo equivalente a clear A. Se ejecuta Tabla2 según: >> Tabla2 >> A A = 1.0000 2.0000 3.0000

1.0000 0.5000 0.3333

1.0000 4.0000 9.0000

1.0000 1.4142 1.7321

3.5. continue Si aparece un continue en un lazo for end o while end, el cálculo pasa a la siguiente iteración del end de ese bucle. 3.6. break Si aparece un break en un lazo for end o while end, el cálculo se puentea a la siguiente instrucción del end de ese bucle, es decir se finaliza ese bucle. 3.7. switch end switch ejecuta un grupo determinado de sentencias basado en el valor de una variable o expresión: switch Expresión % Escalar o de caracteres. case Valor 1 Sentencias % Se ejecutan si Expresión igual Valor 1. case Valor 2 Sentencias . . . otherwise Sentencias % Se ejecutan para Valor no contemplado. end

4 Gráficos en MatLab

4.1. TIPOS DE GRÁFICOS Desde las primeras versiones, MatLab traía suficientes utilidades gráficas, que en las versiones posteriores fueron incorporando cantidades ingentes de nuevas facilidades. Se reseñan los principales tipos de construcción de gráficos: — fplot. Para representación de funciones: fplot('Fun',[-pi pi]). — plot. Representación de x frente a y: plot(x,y). — plotyy. Representación en los ejes opuestos de ordenadas: plotyy (x1,y1,x2,y2). — plotmatrix. Matriz de representaciones: plotmatrix(x,y). — bar. Representación con barras: bar(x,y,ancho,tipo). — stairs. Representación en escalones: stairs(x,y). — errorbar. Representación acompañada de un parámetro de desviación: errorbar(x,y,e). — stem. Representación discreta: stem(x,y). — pie. Representación en tarta: pie(x). — plot3. Representación en 3-d: plot3(x,y,z). — semilogy. Representación semilogarítmica en el eje y: semilogy (x,y). — semilogx. Representación en el eje x: semilogx(x,y). — loglog. Representación logarítmica en los dos ejes: loglog(x,y). En las Figuras 4.1 y 4.2, se representan los principales tipos de gráficos mencionados. 41

42


6

5

4 >> subplot(221) >> fplot(’Fun’,[pi pi]) 2

>> subplot(222) >> stairs(x,Fun(x),’*’)

0

0

–2 –4 –6 –2

0

2

5

–5

–2

0

2

5

1

0

0

>> subplot(223) >> plot(x,Fun(x),’:.’) 0

–5

–5

–2

0

2

–4

–2

0

–1

2

4

Figura 4.1. Diferentes tipos de gráficos. 5

6

>> subplot(222) 4 >> errorbar(x,Fun(x),e)

>> subplot(221) >> bar(x,Fun(x),x,.2)

2 0

0 –2 –4

–5

–6 –2

5

0

2

–2

0

>> subplot(223) >>stem(x,Fun(x))

120 150

2

90 0,5 0,4 60 0,3 0,2

180

0

0

330

210 240 –5 –2

0

30

2

Figura 4.2. Matriz de gráficos.

270

300

GRÁFICOS EN MATLAB

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

43

x=-pi:.25:pi; subplot(221), fplot('fun',[-pi pi]) axis([-pi pi -6 6]) subplot(222), stairs(x,fun(x),’*-’) axis([-pi pi -5 5]) subplot(223), plot(x,fun(x),':.') axis([-pi pi -5 5]) subplot(224), plotyy(x,fun(x),x,x.^2) subplot(221), bar(x,fun(x),0.2) axis([-pi pi -5 5]) e=rand(length(x),1) subplot(222), errorbar(x,fun(x),e) axis([-pi pi -6 6]) subplot(223), stem(x,fun(x)) axis([-pi pi -5 5]) subplot(224), polar(x,sin(2*x).*cos(2*x))

La representación plotmatrix, se muestra en la Figura 4.3, obtenida con las instrucciones: x=[-3:.1:4]'; y=Fun(x); y1=sin(x); y2=tan(x); A=[y,y1,y2]; plotmatrix(x,A,'h-') Título del gráfico 6 Función

4 2 0 –2 –4 1 Seno de X

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

Seno de X 0,5 0

–0,5 –1 20

Tangente de X 0 –20 –3

–2

–1

0 1 Eje de las X

Figura 4.3. Matriz de gráficos.

2

3

4

44


35 30 25 20 15 10 5 0 1 1

0,5 0,5

0

0

0,5

0,5 1

1

Figura 4.4. Gráfico en 3-d.

La representación en 3-d, se efectúa según el ejemplo: >> t = 0:pi/l00:l0*pi; >> plot3(cos(t).*exp(-0.05*t),sin(t).*exp(-0.05*t),t);

4.2. UTILIDADES DE GRÁFICOS Las principales utilidades de los gráficos son: • line([xmin xmax],[ymin ymax],'Color','k','LineWidth',2) • xlabel('Leyenda del eje x'),ylabel('Leyenda del eje y'). • text(x,y,'Leyenda'). • gtext('Leyenda'), posicionada con el ratón. • title('Título del gráfico'). • axis([xmin xmax ymin ymax]). • hold on/off, para que el siguiente gráfico se presente sobre el anterior. • subplot(abc), para representar una matriz de gráficos. • axes('Posición',[x y],'XColor','r','YColor','b') • clf, borrar figuras. • legend, para insertar leyendas, legend(leyenda1,leyenda2..., posición) • datetick(eje,formato de fecha), eje puede ser 'x', 'y' o 'z', el forma-

to de fecha es un número del 0 al 28, 29 posibilidades de expresar las fechas, véase help datetick.

GRÁFICOS EN MATLAB

6%

45

a=[1 3 5 7] pie(a,a==a(2))

19%

44%

31%

Figura 4.5. Gráfico en tarta.

4.3. TEXTOS EN GRÁFICOS Para poder mostrar ecuaciones en cuadros gráficos, ya sea en las leyendas de los ejes, o en otras partes del cuadro, MatLab incorpora un subconjunto de comandos LATEX, y una colección de símbolos de uso frecuente en las ecuaciones, entre otros se encuentra el alfabeto griego, minúsculas y mayúsculas. El conjunto de caracteres puede verse en la documentación en línea de MatLab, en la sección de Handle Graphics. La sintaxis de LATEX utiliza los argumentos entre llaves, «{ }». Los subíndices deben ir precedidos de «_», los superíndices van precedidos de «^». Para conseguir el símbolo de las letras griegas, ha de escribirse su nombre en inglés precedido de «\», con la primera letra en minúscula o mayúscula. En el siguiente ejemplo se muestra cómo poner varias líneas de texto en las leyendas de los ejes, cómo incorporar una ecuación en un eje y en el interior de cuadros gráficos. >> >> >> >>

x=-10:0.01:10; y=sin(x).*exp(-0.2*x); plot(x, y) xlabel({'\fontsize{12}Eje de las X','Distancia \mum'}) ylabel('\fontsize{16}y=sin(x)\cdote^{-\phi \cdot x}') text(- 10.5,'\fontsize{20}y= \mu_{\barpi} \cdot ... sin(x) \cdot e^{-\phi \cdot x}')

46


5 4 3

y=sin(x) ⋅ e - Φ ⋅ x

2

y= μ ϖ ⋅ sin(x) ⋅ e -Φ ⋅ x

1 0 –1 –2 –3 –4 –5 –10

–8

–6

–4

–2

0

2

4

6

8

10

Eje de las X Distancia μ m

Figura 4.6. Ecuaciones en el interior de cuadros gráficos.

4.4. LAPRINT Al incluir cuadros gráficos procedentes de MatLab, en textos LATEX, escalados ya con resizebox o scalebox, suele haber recortes indeseados y desproporción en los textos de los cuadros gráficos en los documentos finales.

6 4 2 0 –2 –4 –6 1

2

3

4

5

6

7

8 1

2

3

4

5

6

Figura 4.7. Gráfico de barras, bar3(peaks(8)).

7

8

GRÁFICOS EN MATLAB

47

LaPrint es un archivo, laprint.m, cuya versión más reciente puede obtenerse en www.uni-kassel.dr/~linne, que reemplaza todas las anotaciones de una figura de MatLab por marcas, salvando la figura como un archivo.eps, a la vez que crea un archivo.tex, para reproducir la figura original, usando epsfig y psfrag, de modo que la figura incorporada en el documento final sea idéntica, incluidas sus fuentes de texto, a la figura original en MatLab. Para utilizar laprint.m, es necesario copiarlo a cualquier carpeta del matlabpath. Una vez que se ha construido el cuadro gráfico en MatLab, se ejecuta laprint.m en la ventana de MatLab. El archivo.tex producido, se inserta en el documento LATEX, en donde convenga, y el archivo.eps en la carpeta correspondiente. 4.5. ESTILOS DE LÍNEAS, MARCAS Y COLORES Los estilos de líneas, marcas y colores se presentan en la Tabla 4.1. Símbolo b g r c m y k w

Color Azul Verde Rojo Cian Magenta Amarillo Negro Blanco

Símbolo

Marca

Símbolo

· o x + * s d ∧ M < > p h

· o x + * ⵧ 䉫䉭䉮䉰䉯夽夹

– : –· ––

Estilo de línea Continua Punteada Trazo-Punto Trazo-Trazo

Cuadro 4.1. Características de las líneas gráficas.

En la Figura 4.8 se muestra un ejemplo de modificación de atributos de líneas como son: color de línea, grueso... % ----------------Lineas.m------------------------------------t=linspace(-4*pi,4*pi,1000); y=sin(t).*exp(-.25*abs(t)); z=cos(t).*exp(-.25*abs(t)); plot(t,y,'LineWidth',2,'Color','k'), grid axis([-4*pi,4*pi,-.7,1]) text(2,sin(2).*exp(-.25*abs(2)),... '\fontsize{16}\leftarrow sin(t)\cdote^{-0.25\cdotabs(t)} ')

48


xlabel('\fontsize{l5} t') ylabel('\fontsize{l5} sin(t)\cdote^{-0.25\cdotabs(t)}') hold on a=plot(t,z); text(0.25,cos(0.25).*exp(-.25*abs(0.25)),... '\fontsize{14}\leftarrow cos(t)\cdote^{-0.25\cdotabs(t)} ') set(a,'LineWidth',0.5,'Color','k') t1=-11; y1=sin(t1).*exp(-.25*abs(t1)); t2=0; y2=cos(t2).*exp(-.25*abs(t2)); line([t1 t2],[y1 y2],'LineWidth',8,'Color',[1 0 1])

1

← cos(t)⋅e0,25 ⋅abs(t) 0,8 0,6

← sin(t)⋅e 0,25⋅abs(t)

sin(t)⋅e 0,25⋅ abs(t)

0,4 0,2 0 –0,2 –0,4 –0,6 –10

–5

0 t

5

10

Figura 4.8. Gruesos de líneas, color... hold of % -----------------------------------------------------------------

4.6. area area es una función análoga a plot, rellena el espacio comprendido entre 0 y una línea gráfica. En la Figura 4.9, se rellena la superficie limitada entre 0 y Fun(x), gráfico superior, y entre –5 y Fun(x), gráfico inferior:

GRÁFICOS EN MATLAB

49

6 4

subplot(211) area(x,Fun(x))

2 0 –2 –4 –6 –2

–1

0

1

2

3

4

5

0

1

2

3

4

5

6 4 2

subplot(212) area(x,Fun(x),5)

0 –2 –4 –6 –2

–1

Figura 4.9. Relleno de áreas en gráficos.

33%

a=[3 5 7 9 12]; pie3(a,a==a(4))

8%

14%

25%

19%

Valor destacado

Figura 4.10. Gráfico en tarta en 3-d.

5 Aplicaciones de cálculo numérico

5.1. INTEGRACIÓN NUMÉRICA Matlab dispone de dos funciones para integración, quad, basado en la regla de Simpson, y quadl, evaluación mediante la técnica de cuadratura de Lobatto. En el siguiente ejemplo se muestra el cálculo de la integral de la función Fun, previamente definida, desde el límite inferior 0, hasta 5. >> A=quad(F,0,5)

% Integración por la regla de Simpson. % Siendo F=@Fun, definida en Cap. 1.

A = -10.0814

5.2. MÍNIMOS, CEROS Y OPTIMIZACIÓN Las principales funciones de MatLab para optimización son: fminbnd, que encuentra el mínimo de una función de una variable; fminsearch, correspondiente a fmins de versiones anteriores, calcula mínimos de funciones multivariables; lsqcurvefiting para ajustes de datos a diferentes tipos de funciones; y fzero, que encuentra ceros de funciones de una variable. 5.2.1. Mínimos y ceros de funciones Para encontrar un mínimo de la función, Fun, entre 0 y 3, se procede según: 51

52


>> x=fminbnd(F,0,3) x = 2.0351

La siguiente instrucción busca un cero en las proximidades de 4: >> fzero('Fun',4) ans = 3.0965 >> fzero(F,4) ans = 3.0965

5.2.2. Resolución de ecuaciones no lineales La resolución de ecuaciones no lineales se efectúa con la instrucción fsolve, cuya aplicación se realiza según el ejemplo que sigue. Sea el sistema: Sea el sistema

⎫ ⎪ ⎬ ⎪ = 7, 389⎭

x+y=5 x − y = –1 e z– x

La solución se consigue con el programa EcusP.m % --------------- EcusP.m ---------------------------clear Opciones=optimset ('MaxFunEvans',5000,'GradConstr','on'... 'TolCon',1e–5,'TolFun',1e–6,'TolX'1e–5); x=fsolve('Ecus',[40 50 –5]',Opciones) % ------------------------------------------------------------

En la fución Ecus se definen las ecuaciones a resolver, según el listado: % --------------- Ecus.m ---------------------------function q=Ecus(p) x=p(1); y=p(2); z=p(3); q(1)=x+y–5; q(2)=x–y+1; q(3)=exp(x.*z)–7.389; q=[q(1);q(2);q(3); % ------------------------------------------------------------

APLICACIONES DE CÁLCULO NUMÉRICO

53

La ejecución de EcusP produce el siguiente resultado: > EcusP Optimización terminated successfully: Relative function value changing by less than OPTIONS.TolFun x = 2.000 3.000 1.000 >>

5.2.3. Minimización y ajuste de datos El objetivo de fminsearch es encontrar un mínimo de una función multivariable sin restricciones. Se invoca de la forma: pk=fminsearch('pHAjusMin',[3 l0],Opciones)

En este caso, se llama a la función a minimizar, pHAjusMin, suministrando a continuación los valores iniciales de los parámetros a minimizar, siendo opcional añadir una variable de opciones, seleccionables con optimset. La aplicación y funcionamiento de fminsearch se comprende fácilmente con el siguiente ejemplo de ajuste de datos experimentales a una función, ecuación (5.1). En el proceso de lodos activos, la tasa de crecimiento de los lodos es función del pH, datos experimentales. Las reacciones de biodegradación, actividad microbiana, transcurren a un pH óptimo, próximo a la neutralidad y disminuyen hasta anularse, al separarse a zonas agresivas, como son las zonas de pHs ácidos o alcalinos. El efecto del pH en la velocidad específica de generación de lodos se representa adecuadamente por una función tipo campana: μ ( pH ) =

1 1 + 10

pk1 − pH

+ 10 pH − pk2

(5.1)

Las constantes pk1 y pk2, se calculan por ajuste de datos, representando los valores en los que μ(pH), tiene el valor de 0,5.

54


1 0,9 pk1 = 5,2036 0,8 pk2 = 10,1128 0,7

μ(pH)

0,6 0,5 0,4 0,3 0,2 0,1 0

0

2

4

6

8

10

12

14

pH

Figura 5.1. Ajuste de la tasa de crecimiento frente al pH.

El listado del ajuste mediante la ecuación (5.1), se presenta a continuación: % --------------- pHAjus.m ---------------------------clear,clf global pH mu pH=[O 1 2 3 4 5 6 7 8 9 10 11 12 13 14]'; mu=[O .01 .02 .04 .12 .4 .8 1 .96 .9 .6 .05 .01 ... .01 0]'; pk=fminsearch('pHAjusMin',[3 10]); plot(pH,mu,'*') hold on pHc=0:.1:14; muc=1../(1+1O..ˆ(pk(1)-pHc)+10..ˆ(pHc-pk(2))); plot(pHc,muc,'-') mu1=1../(1+10..ˆ(pk(1)-pk(1))+10..ˆ(pk(1)-pk(2))); mu2=1../(1+10..ˆ(pk(1)-pk(2))+10..ˆ(pk(2)-pk(2))); line([pk(1) pk(1)],[0 mu1],'Color','k') line([pk(2) pk(2)],[0 mu2],'Color','k') text(1,0.85,'pk1 = '), text(2.5,.85,num2str(pk(1))) text(1,0.75,'pk2 = '), text(2.5,.75,num2str(pk(2))) xlabel('pH'), ylabel('\mu(pH)'), grid, hold off %--------------------------------------------------------

El listado de la función a minimizar es: function q=pHAjusMin(p) %----------------- pHAjusMin.m -------------------% Función llamada por pHAjus.m %--------------------------------------------------


55

global pH mu pk1=p(1); pk2=p(2); muc=1../(1+10..ˆ(pk1-pH)+10..ˆ(pH-pk2)); q=sum((mu-muc).ˆ2); %--------------------------------------------------

En la función pHAjusMin.m, se define q, parámetro a minimizar; en este caso minimizar la suma de los cuadrados de las desviaciones de los μ experimentales a los calculados. La función lsqcurvefit es más específica para el ajuste de datos; se ejecuta según: X=lsqcurvefit('FUN',XO,X,Y,LI,LS,Opciones)

Se invoca a la función FUN, en la que se define la ecuación de ajuste, XO representa los valores iniciales de partida, X e Y son los datos experimentales, LI representa el límite inferior de los valores de los parámetros a calcular, LS es el límite superior, y Opciones es un conjunto de valores en los que se definen parámetros del cálculo, siendo opcional su especificación. A continuación se muestra un ejemplo de la utilización de lsqcurvefit: >> >> >> >> >> >>

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

x=[-pi*2:.5:pi*2]'; A=2; B=.2; y=A*sin(x).*exp(-B*x); r=rand(length(x),1)-0.5; Y=y+r; Opciones=optimset('Display','iter','Diagnostics','on',... 'TolX',1e-29,'TolFun',1e-29,'LargeScale','on',... 'MaxFunEvals',100) est=lsqcurvefit('AjusNlFun',[0 0],x,Y,[-5 -5],[10 5],... Opciones) Ac=est(1); Bc=est(2); xc=linspace(x(1),x(length(x))); yc=Ac*sin(xc).*exp(-Bc*xc); plot(x,y,'o'), hold on, plot(xc,yc,'-'), hold off text(0,5,'A='), text(1,5,num2str(A)) text(0,4,'B='), text(1,4,num2str(B)) text(4,5,'Ac='), text(5,5,num2str(Ac)) text(4,4,'Bc='), text(5,4,num2str(Bc)) xlabel('X'), ylabel('A.sen X ·eˆ{-B·X}')

La función llamada se muestra en el siguiente archivo: function F=AjusNlFun(p,x) %---------------AjusNlFun.m--------------A=p(1); B=p(2); F=A*sin(x).*exp(-B*x); %-----------------------------------------

56


En la Figura 5.2 se muestra el resultado de este ajuste, en la que se muestran los valores de partida y los encontrados, teniendo en cuenta que a cada Y, se le sumó aleatoriamente una cantidad comprendida en el intervalo ±0.5. 6 5

A= 2

Ac=2,096

4

B=0,2

Bc=0,19176

A·sen X·e –B·X

3 2 1 0 –1 –2 –3 –8

–6

–4

–2

0 X

2

4

6

8

Figura 5.2. Ajuste de datos con Isqcurvefit.

En el ajuste de datos experimentales a funciones es frecuente probar con polinomios, para ello MatLab desarrolló la función polyfit, cuya utilización se muestra en el ejemplo siguiente. 30 25 Ac=0,95551 Bc= 4,6119 Cc= 5,6322

x 2+5⋅x+6

20 15 10 5 0 –5 –8

–7

–6

–5

–4

–3 x

–2

–1

Figura 5.3. Ajuste de datos con polyfit.

0

1

2


57

Se definen 10 valores de x, linealmente espaciados, a los que se le aplica una función polinómica para obtener unas y, a las que se les suma un ruido aleatorio de ± 2.5, con lo que se obtienen Y valores. Con la función polyfit sobre x e Y, se pretende recalcular los coeficientes del polinomio original, según las siguientes instrucciones: >> >> >> >> >> >> >> >> >>

x=linspace(-8,2,10)'; A=1; B=5; C=6; y=A*x.ˆ2+B*x+C; Y=y+5*rand(length(x),1)-2.5; Est=polyfit(x,Y,2); xx=x(l):0.1:x(length(x)); plot(x,y,'*',xx,Est(1)*xx.ˆ2+Est(2)*xx+Est(3)) text(-5,20,'Ac='), text(-4.5,20,num2str(Est(1))) text(-5,18,'Bc='), text(-4.5,18,num2str(Est(2))) text(-5,16,'Cc='), text(-4.5,16,num2str(Est(3)))

5.3. INTEGRACIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES Las ecuaciones diferenciales que representan a los procesos reales suelen ser ecuaciones diferenciales no lineales, impidiendo su integración simbólica, de modo que es necesario recurrir a métodos numéricos, cuya solución útil será una tabla de valores o la representación gráfica de éstos. MatLab dispone de las funciones ode23, ode45, ode113, ode15s, ode23s, ode23t y ode23tb para la evaluación numérica de ecuaciones diferenciales. Las funciones más utilizadas son ode23 y ode45, basadas en el método de Runge-Kuta de 2 y 4 parámetros respectivamente, algunas de las otras funciones se utilizan para ecuaciones con rigideces, stiff. Las funciones de integración mencionadas son de paso de integración variable, evaluándose en cada iteración la solución con el paso de integración h, h/2 y 2h, si los resultados no superan una tolerancia determinada, el paso de integración se reduce a h/2, si la superan el paso de integración se incrementa, y así sucesivamente. En opciones de integración se posibilita limitar el paso de integración a topes máximo y mínimo. En ocasiones ocurre que integrando una ecuación diferencial en un determinado valor existen diferencias significativas en la solución encontrada, aún reduciéndose sucesivamente el paso de integración, que si se supera un número determinado de reducciones, el sistema se para a causa de valores singulares, pudiendo ser de utilidad disponer de una función de integración de ecuaciones diferenciales de paso fijo, definido por el usuario, para lo cual se elaboró la función odegil4, cuyo listado y aplicaciones se muestra más adelante. Los métodos numéricos evalúan únicamente ecuaciones diferenciales de primer orden, de modo que para evaluar ecuaciones diferenciales de orden superior, han de definirse variables auxiliares para componer sistemas

58


de ecuaciones diferenciales de primer orden, equivalentes a las de orden superior. En el Capítulo anterior se resolvió simbólicamente la ecuación: d2y dy +3 + 6y = 0 2 dt dt

(5.2)

para su resolución numérica se realiza la siguiente transformación: z = y' con lo que la ecuación (5.2) se transforma en dos ecuaciones diferenciales de primer grado: z' = –3z – 6y

(5.3)

y' = z

(5.4)

Se pretende integrar desde –π hasta π, para lo cual es necesario conocer los valores de y(–π) y de dy(–π), que se obtienen según: >> >> >> >>

t=sym('t'); y=('2*15ˆ(1/2)/5*sin(15ˆ(1/2)/2*t)*exp(-3/2*t)'); dy=diff(y); subs(y,t,-pi)

ans = 34.1795 >> subs(dy,t,-pi) ans = 276.0593

A continuación, se presenta el listado de la llamada a integración de las ecuaciones consideradas y su representación gráfica: %--------------------- EcudifP.m ------------------y0=[276.0593 34.1795]; t=-pi:.01:pi; [t,y]=odegil4(’Ecudif’,t,y0,0.01); plot(t,y(:,1),’k:’,t,y(:,2),’k-’) xlabel('Tiempo'), ylabel('Funcion y derivada') axis([-pi pi -100 277]) %---------------------------------------------------


59

250

Función y derivada

200 150

y’

100 y

50 0 –50 –100

–3

–2

–1

0 Tiempo

1

2

3

Figura 5.4. Gráfico del resultado de la integración numérica.

En el siguiente listado se muestra el programa de definición de las ecuaciones a integrar: function dy=Ecudif(t,y) -------------------- Ecudif.m ----------------------------% y''+3y'+6y=0 % es equivalente al sistema: % z'=-3z-6y % y'=z %----------------------------------------------------------dy(1)=-3*y(1)-6*y(2); % y(1)=z dy(2)=y(1); % y(2)=y dy=[dy(1);dy(2)]; %-----------------------------------------------------------

5.3.1. Método de Runge-Kuta El método más común de resolución numérica de ecuaciones diferenciales es el método de Runge-Kuta de cuarto orden, cuya aplicación se resume a continuación. Sean las ecuaciones diferenciales tales como:

60


d y1 = ƒ1 (t, y1 , y2 ,… yn ) dt d y2 = ƒ 2 (t, y1 , y2 ,… yn ) dt M d yn = ƒ n (t, y1 , y2 ,… yn ) dt

Su integración se realiza evaluando los siguientes parámetros: pk = h ⋅ ƒ(tk , yk ) h p qk = h ⋅ ƒ ⎛ tk + , yk + k ⎞ ⎝ 2 2⎠ h q rk = h ⋅ ƒ ⎛ tk + , yk + k ⋅⎞ ⎝ 2 2 ⎠ sk = h ⋅ ƒ(tk + h, yk + r )

Las ecuaciones integradas numéricamente se obtienen mediante: yk ,i+1 = yk ,i +

pk + 2 ⋅ qk + 2 ⋅ rk + sk 6

Para integración de ecuaciones diferenciales por el método de Runge-Kuta de cuatro parámetros se elaboró la función odegil4.m de paso de integración fijo: function [t,y]=odegil4(f,x,y0,h) %------------------------------------------------------------------% [t,y]=odegil4(f,x,y0,h), resuelve ecuaciones y'=f(t,y), % x vector en linea de la variable independiente, % y0 vector en linea de las condiciones iniciales, % h paso de integración. %------------------------------------------------------------------x=x'; t=[]; t(1)=x(1); y=y0;


61

N=(x(length(x))-x(1))/h; Ne=length(y0); for I=1:N T=t(I); Y=y(I,:); for J=1:Ne p=h*feval(f,T,Y'); q=h*feval(f,T+h/2,Y'+p/2); r=h*feval(f,T+h/2,Y'+q/2); s=h*feval(f,T+h,Y'+r); end t(I+1)=t(1)+h*I ; Inc=(p+2*q+2*r+s)/6; y(I+1,:)=y(I,:)+Inc’; end y=spline(t,y',x); y=[reshape(y,Ne,length(x))]'; t=x; %-------------------------------------------------------------

La aplicación práctica se materializa en el siguiente ejemplo, en el que se evalúa a modo de demostración la integración de 2t, 3t2, 4t3 y cos(x), cuyas soluciones son: t2, t3, t4 y seno(t). Las ecuaciones diferenciales se definen en la función ecudif.m: function dy=ecudif(t,y) dy(1)=2*t; dy(2)=3*t.ˆ2; dy(3)=4*t.ˆ3; dy(4)=cos(t); dy=[dy(1);dy(2);dy(3);dy(4)];

Para resolver el problema se lanza con las siguientes instrucciones: >> x=[0:5]; % Valores de presentación de resultados. >> y0=[0 0 0 0]; % Condiciones iniciales. >> h=0.01; % Paso de integración. >> tic % Comienzo de temporización. >> [t,y]=odegil4('ecudif',x,y0,h); % Llamada a integración. >> toc % Final de la temporización. elapsed_time = 1.2010 >> [t,y] ans = t 0 1.0000 2.0000 3.0000 4.0000 5.0000

% Presentación de resultados. tˆ2 0.0000 1.0000 4.0000 9.0000 16.0000 25.0000

tˆ3 0.0000 1.0000 8.0000 27.0000 64.0000 125.0000

tˆ4 seno(t) 0.0000 0.0000 1.0000 0.8415 16.0000 0.9093 81.0000 0.1411 256.0000 -0.7568 625.0000 -0.9589

62


5.4. MODELO DINÁMICO DE UN TANQUE Cuando un tanque se alimenta con un caudal variable de un líquido hasta que se llena, y comienza a salir el líquido por rebose, surge una lámina del líquido por encima del borde de vertido. Esta lámina adquiere un espesor dependiente del caudal de entrada y de las dimensiones del tanque. A su vez, la altura de esa lámina provoca la cuantía del caudal de salida. El caudal de entrada es una variable impuesta, entrada, en simulación se le asignarán diferentes valores para observar la respuesta del sistema. El volumen indicará el estado del sistema. La altura de lámina y caudal de salida constituyen la salida del sistema. Siguiendo la nomenclatura del Análisis de Sistemas, este sencillo proceso se muestra en la Figura 5.5.

U

Y

X V

Qi

Qe h

Figura 5.5. Esquema general de un sistema dinámico.

U representa la entrada al sistema, caudal de entrada. X reprepresenta el estado del sistema, volumen del líquido en el tanque. Y representa las salidas del sistema, altura de lámina en el tanque y caudal de salida. El esquema físico del sistema considerado se representa en la Figura 5.6.

L h Qe

Qi

V H A

Figura 5.6. Esquema de un tanque con rebosadero lateral.


63

Considerando que el líquido circulante es agua, el caudal que sale por un vertedero rectangular, [m3/s], de longitud L [m], viene dado por la expresión de Bazin, ecuación (5.5), de uso general en Francia: Q = μ ⋅ L ⋅ h1,5 ⋅ 2 ⋅ g

(5.5)

μ es un coeficiente que depende de la altura de lámina h [m], y de la profundidad del tanque H [m], según la ecuación (5.6). μ = 0, 405 +

0, 003 ⎡ h2 ⎤ ⋅ ⎢1 + 0, 55 ⋅ h ( H + h)2 ⎥⎦ ⎣

(5.6)

El caudal de salida, calculado con este procedimiento, es válido para alturas de lámina de 2,5 a 80 cm. La variación del volumen de agua en el tanque, es la diferencia de los caudales de entrada al de salida, según la ecuación (5.7). dV = Qi − Qe dt

(5.7)

La altura de la lámina será el volumen total de agua entre la sección, menos la altura del rebosadero, según la ecuación (5.8): h=

V −H A

(5.8)

Para resolver el ejemplo numérico utilizando un programa en MatLab, se dimensiona el tanque con área de la base de 150 m2, longitud de vertedero 10 m y altura 4 m. El caudal medio de entrada se fija en 100 m3/h, con variaciones de ± 50%, como es el caso de un reactor biológico, alimentado con caudales mínimos en horas nocturnas y máximo en horas diurnas. El tiempo de análisis para este ejemplo se fijó en 20 minutos, suficiente para observar esta dinámica. El programa de cálculo se estructura en tres partes, un programa y dos funciones. La primera función calcula el estado del sistema, volumen de agua en el tanque X, en función de la diferencia de caudales de entrada, impuesto, y del de salida. La segunda función calcula las salidas Y, altura de lámina y caudal de salida. En el programa principal se definen los datos de partida, desde el que se llama a la primera función, y ésta a la segunda, para finalmente realizar los cálculos requeridos y presentación gráfica de resultados.

64


El listado del programa principal es: %--------------------- Vertedero.m ----------------------% clear, clf global A H L Qe1 Qe2 Qe3 tf % % Nomenclatura %--------------------------------------------------------% % A Superficie de la base % H Altura del vertedero % L Longitud del vertedero % Qe1 Caudal del primer tramo % Qe2 Caudal del segundo tramo % Qe3 Caudal del tercer tramo % V Volumen inicial de agua en el tanque % t0, tf Tiempos inicial y final de análisis % % X Volumen variable del agua % U Caudal de entrada % Y(:,1) Caudal de entrada % Y(:,2) Caudal de salida % Y(:,3) Altura de lámina % % Parámetros y Estado Inicial %--------------------------------------------------------% A=150; H=4; L=10; tf=20*60; Qe1=50/3600; Qe2=100/3600; Qe3=150/3600; V=590; % % Integración de la ecuación de estado %--------------------------------------------------------% Opciones=odeset('RelT',1e-7,'AbsTol',1e-7); [t,X]=ode23('VerteX',[0,tf],V,Opciones); % % Cálculo de las Salidas %--------------------------------------------------------% Y=verteY(t,X); % % Representación de resultados %--------------------------------------------------------% T=t/60; tff=T(length(T)); % subplot(311); plot(T,Y(:,1)*3600,'*-',T,Y(:,2)*3600,'o--') axis([t0 tff -5 154]), ylabel('Cuadales (m3/h)') % subplot(312);plot(T,Y(:,3)*1000,'*-'),axis([t0 tff -1 14]) ylabel('Altura de Lámina (cm)')


65

% subplot(313); plot(T,X,'*-'), axis([t0 tff 590 603]) ylabel('Volumen del Agua (m3)'),xlabel('Tiempo (Minutos)') % %---------------------------------------------------------

La función que calcula el estado del sistema, volumen del tanque, es: function dX=VerteX(t,X) %--------------------- VerteX.m--------------------------% % Cálculo del volumen del tanque %--------------------------------------------------------% % global A H L Qe1 Qe2 Qe3 tf % % Cálculo del caudal de salida %--------------------------------------------------------% Y=VerteY(t,X); U=Y(1,1); % % Variación de volumen del tanque %---------------------------------------------------------

Caudales (m3 /h)

150 100 50

Volumen del Agua (m3 )

Altura de Lamina (cm)

0 0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10 12 Tiempo (Minutos)

14

16

18

20

10 5 0

600

595

590

Figura 5.7. Caudales de entrada, salida, altura de lámina y volumen de un tanque con caudal de entrada escalonada.

66


% dX=U-Y(1,2); % %---------------------------------------------------------

El listado de la función que calcula las salidas del sistema es: function Y=VerteY(t,X) %---------------------- VerteY.m ------------------------% % Cálculo de Caudal de Salida y Altura de Lámina % %--------------------------------------------------------% % global A H L Qe1 Qe2 Qe3 tf % N=length(t); Y=zeros(N,2); for i=1:N % Caudal de Entrada if t(i)
En la Figura 5.7 se representan los resultados de la simulación obtenida, caudal de salida, altura y volumen de lámina de un tanque alimentado con un caudal de entrada escalonado. Observando la Figura 5.7 se aprecia que la dinámica de este proceso no es instantánea; al aumentar el caudal de entrada, la altura de lámina aumenta de forma gradual, tardando en este ejemplo alrededor de 5 minutos en alcanzar el estado estacionario. El caudal de salida, por ser función de la altura de la lámina, lleva una dinámica análoga.


67

5.5. DETERMINACIÓN DE RETRASOS Y DERIVADAS En los procesos de flujo, como en el caso precedente, hay una variable de salida, caudal de salida, retrasado con respecto al caudal de entrada mientras hay rebose. En otros ejemplos de la industria es muy frecuente encontrar variables sometidas a retrasos, mostrándose a continuación un método aproximado de cálculo numérico de variables afectadas por retrasos. El cálculo de una variable sujeta a un retraso se calcula teniendo en cuenta que la transformación de Laplace de un retraso puro es: [ f(t – τ)] = e–s · τ · F(s) Operador Transformada de Laplace. Siendo: τ Tiempo de retraso. s Variable compleja. Las aproximaciones de retrasos más comunes son las aproximaciones de Pade de primer y segundo orden:

e

e − s⋅τ ≈

− s⋅τ

τ ⋅s 2 ≈ τ 1+ ⋅s 2 1−

τ 2 ⋅ s 2 − 6 ⋅ τ ⋅ s + 12 τ 2 ⋅ s 2 + 6 ⋅ τ ⋅ s + 12

En procesos químicos en los que ocurren tiempos muertos, como en una colunma de rectificación o retrasos hidráulicos, es más útil la aproximación menos común: e − s⋅τ ≈

1 ⎛ τ ⋅ s + 1⎞ ⎝n ⎠

n

La aproximación más sencilla de aplicar, a la vez que provoca error mínimo, salvo en los instantes iniciales, es la primera aproximación mencionada. τ y( s ) 1 − 2 ⋅ s = x(s) 1 + τ ⋅ s 2

(5.9)

68


Siendo:

x(s) y(s)

Variable de la que se realiza el retraso. Variable igual a x(s) retrasada τ.

Realizando la transformación inversa de la ecuación (5.9), se obtiene: y+

τ τ ⋅ y' = x − ⋅ x' 2 2

De donde se obtiene la derivada de la variable retrasada: y' =

2 ⋅ ( x − y) − x' τ

(5.10)

Esta ecuación integrada conduce a: y( t ) = x ( t − τ ) =

2 ∫ ( x (t ) − y(t )) ⋅ dt − x (t ) τ

El cálculo implicado en diagrama de bloques se muestra en la Figura 5.8. x(t) 2/

–

x(t- )

–1

Figura 5.8. Diagrama de bloques de un retraso de primer orden.

Para el cálculo de las variables retrasadas, haciendo uso de la integración numérica de la ecuación (5.10), se necesitan conocer las derivadas de las variables a las que se someten los retrasos. Las derivadas se calculan basándose en la propiedad de que, en el dominio de Laplace, multiplicar por s, al pasar al dominio del tiempo es derivar: d ƒ (t ) ⎤ = s ⋅ F( s) − ƒ(0) ⎡⎢ ⎣ dt ⎥⎦

f (0) se suele despreciar, ya que únicamente tiene influencia en los momentos iniciales. La aproximación para el cálculo de derivadas se muestra en la


69

Figura 5.9, en donde a representa una variable a la que se le calcula su derivada, representando b su derivada. a

s

b

1+0,01 s

Figura 5.9. Diagrama de bloques del cálculo aproximado de derivadas.

De la Figura 5.9, se obtiene: a · s = 0,01 · s · b + b Pasando al dominio del tiempo se obtiene: a' = 0,01 · b' + b Despejando b′ e integrando se obtiene:

(

b = 100 ⋅ a – ∫ b ⋅ dt

)

Este es un procedimiento aproximado para obtener la derivada de una función a través de su integral. Para visualización rápida de la aproximación del cálculo de variables retrasadas y derivadas, se elabora un programa de demostración en el que la función de entrada es el tiempo al cuadrado, y la función retrasada, el tiempo al cuadrado retrasado en tres unidades. En el programa principal, RetraLan.m se definen los parámetros de operación, tiempo, retraso, condición inicial, integrador, la representación gráfica y la tabla de resultados numéricos. %------------------RetraLan.m ---------------------t=0:.1:12; global tau tau=3; y10=0; y20=0; [t,y]=odegil4('Retraso',t,[y10 y20],0.01); Y=y(:,1); Yr=y(:,2); subplot(211),plot(t,Y), grid axis([0 12 -2 145])

ylabel('\fontsize{12} tˆ2') subplot(212),plot(t,Yr), grid

70


axis([0 12 -3 82]) ylabel('\fontsize{12} (t-3)ˆ2') xlabel('\fontsize{12} t) T=t(1):t(length(t)); Yn=round(spline(t,Y,T)); Ynr=round(spline(t,Yr,T)); [T’ Yn’ Ynr’] %---------------------------------------------------

En la función Retraso.m, se define la función de partida F, de la que se obtiene, a modo de exposición, el cálculo de la derivada numérica, dy(1)=100*(F-y(1)), de la que se obtiene la función retrasada y(2). %-----------------Retraso.m -----------------------function dy=Retraso(t,y) global tau F=t.ˆ2; % Funcion. dy(1)=100*(F-y(1)); % Derivada de Funcion. dy(2)=(y(1)-y(2))*2/tau-dy(1); % Funcion retrasada. dy=[dy(1);dy(2)]; %---------------------------------------------------

Los resultados numéricos y gráficos se obtienen según: >> RetraLan ans = t tˆ2 (t-3)ˆ2 0 0 0 1 1 -1 2 4 -1 3 9 -1 4 16 0 5 25 4 6 36 9 7 49 16 8 64 25 9 81 36 10 100 49 11 121 64 12 144 81


71

t2

100

50

0 0

2

4

6

8

10

12

2

4

6 t

8

10

12

80

(t–3)2

60 40 20 0 0

Figura 5.10. Resultado gráfico del cálculo de retrasos.

En el programa principal se hizo uso de la función spline, para obtener datos interpolados. Su sintaxis es: Yn=spline(X,Y,Xn)

Su funcionamiento se pone de manifiesto en el ejemplo siguiente: Se trata de obtener los cuadrados de los números enteros, por interpolación de los cuadrados de los intermedios correspondientes. >> X=0.5:5 X = 0.5000

1.5000

2.5000

3.5000

4.5000

2.2500

6.2500

12.2500

20.2500

>> Y= X.ˆ2 Y = 0.2500 >> Xn=1:5 Xn= 1 2 3 4 >> Yn=spline(X,Y,Xn)

5

Yn = 1

4

9

16

25

72


La interpolación así conseguida es perfecta. En el programa principal se utilizó length, para suministrar el número de elementos de un vector: >> length(Xn) ans = 5

5.6. AJUSTE DE DATOS EXPERIMENTALES A UNA RECTA El ajuste de datos a una recta es un caso particular de polyfit correspondiente al primer grado. Aquí se muestra su cálculo por varias razones: para mostrar el manejo de datos, para practicar la evaluación de sumatorios con MatLab, y porque muchas funciones no lineales pueden transformarse en lineales, siendo la recta la mejor prueba de hipótesis de correlación. Los mejores coeficientes de la recta son los que consiguen hacer mínimo el sumatorio de las desviaciones al cuadrado, de los puntos experimentales a la recta postulada. El método se resume en las ecuaciones (5.11), (5.12), (5.13) y (5.14). La ecuación de la recta buscada se expresa por la ecuación (5.11). Y=m·X+b

(5.11)

Con datos procedentes de medidas experimentales hay una desviación:

δ=Y–m·X–b Elevando al cuadrado la desviación, y extendiendo al conjunto de datos, se obtiene el sumatorio de las desviaciones al cuadrado: Σδ 2 = ΣY 2 + m2 · ΣX 2 + N · b2 – 2 · m · ΣX · Y – 2 · b · ΣY + 2 · m · b · ΣX Los coeficientes se determinan de las condiciones del mínimo: d Σδ 2 = 2 ⋅ m ⋅ ΣX 2 − 2 ⋅ ΣX ⋅ Y − 2 ⋅ b ⋅ ΣX dm d Σδ 2 = 2 ⋅ N ⋅ b – 2 ⋅ ΣY + 2 ⋅ m ⋅ ΣX db

Igualando a cero, y simplificando estas expresiones, se obtiene un sistema de dos ecuaciones con dos incógnitas: m · ΣX 2 + b · ΣX = ΣX · Y m · ΣX + N · b = ΣY


73

Resolviendo el sistema se obtienen los coeficientes de la recta: b=

ΣY ⋅ ΣX 2 − ΣX ⋅ ΣX ⋅ Y N ⋅ ΣX 2 − ( ΣX )2

(5.12)

N ⋅ ΣX ⋅ Y − ΣX ⋅ ΣY N ⋅ ΣX 2 − ( ΣX )2

(5.13)

m=

El grado de bondad del ajuste se analiza mediante el coeficiente de correlación, cuya expresión general toma la forma de: r=±

Σ(Yc − Y )2 Σ (Y − Y ) 2

El coeficiente de correlación específico para la recta es: r=±

N ⋅ ΣX ⋅ Y − ΣX ⋅ ΣY ( N ⋅ ΣX − ( ΣX )2 ) ⋅ ( N ⋅ ΣY 2 − ( ΣY )2 ) 2

(5.14)

Los datos con los que se va a trabajar se guardarán en un archivo separado, con nombre alusivo a su contenido, formando una matriz a la que se la denomina datos, en cuya primera columna están las X, y en la segunda las Y: %----------------- DatosAL.m ----------------------datos=[1 6.23 2 8.58 3 10.84 4 12.48 5 13.35 6 15.81 7 18.87 8 20.83 9 21.82]; %---------------------------------------------------

El programa principal de ajuste a una recta por mínimos cuadrados, según las ecuaciones previas, es: %----------------- Ajulineal.m --------------------DatosAL X=datos(:,1); Y=datos(:,2); N=length(X); SX=sum(X); SX2=sum(X.ˆ2);

74


SY=sum(Y); SY2=sum(Y.ˆ2); SXY=sum(X.*Y); b=(SY*SX2-SX*SXY)/(N*SX2-SXˆ2); m=(N*SXY-SX*SY)/(N*SX2-SXˆ2); r=(N*SXY-SX*SY)/((N*SX2-SXˆ2)*(N*SY2-SYˆ2))ˆ.5; Yc=m*X+b; plot(X,Y,'*',X,Yc,'-') text(1.5,15,'r='), text(2,15,num2str(r)) text(1.5,17,'b='), text(2,17,num2str(b)) text(1.5,19,'m='), text(2,19,num2str(m)) xlabel('X'), ylabel(’Y’) %---------------------------------------------------

22 20 m= 1,975 18 b= 4,4372 16 r= 0,99502 14 12 10 8 6

1

2

3

4

5

6

7

8

9

Figura 5.11. Ajuste de datos a una recta.

En la Figura 5.11 se muestra el resultado numérico y gráfico del ajuste de datos experimentales a una recta; los datos se muestran con el signo «*», y en línea continua la recta calculada. Como r se aproxima a la unidad, indica excelente grado de ajuste. 5.6.1. Ajuste de funciones no lineales por linealización Las funciones no lineales con uno o dos parámetros son susceptibles de linealización por transformaciones de variables, para lograr nuevas variables relacionadas linealmente, a las que se le aplica el ajuste lineal expuesto. Este tratamiento de datos tiene la ventaja de que siempre es más fácil elucidar si una


75

distribución de datos es lineal o no, con relación a comparar distribuciones de datos a curvas diferentes. A continuación se muestra un ejemplo de linealización del ámbito de la destilación, operación básica de la Ingeniería Química. En el equilibrio líquido-vapor de mezclas binarias, no azeotrópicas, la composición del vapor en función de la composiciórı del líquido sigue la ecuación de uso restringido: y=x+

C ⋅ R ⋅ x ⋅ (1 − x ) [ x + R ⋅ (1 − x )]2

(5.15)

en donde: x Fracción molar del componente más volátil del líquido. y Fracción molar del componente más volátil del vapor. R Constante, que indica la posición del máximo en la representación y-x. C Constante, que indica el valor del máximo en la representación y-x. Si en la fracción de la ecuación (5.15) se divide numerador y denominador por (1 – x)2, se obtiene: x 1− x y−x = 2 ⎡ x + R⎤ ⎢⎣1 − x ⎥⎦ C⋅ R⋅

(5.16)

Realizando los siguientes cambios de variables: D=y–x A=C·R U=

x 1− x

la ecuación (5.16) toma la forma de: D=

A ⋅U (U + R)2

Realizando un nuevo cambio: V2 =

U D

(5.17)

76


la ecuación (5.17) toma la forma de: (U + R)2 = A · V2 Ecuación que, reordenada, toma la forma lineal: V=

R U + A A

(5.18)

La ecuación (5.18) relaciona linealmente V con U. Si la distribución de los datos (x,y) sometidos a las transformaciones expuestas siguen una trayectoria lineal, indica que la ecuación (5.15) es consistente. De la recta obtenida se determinan su pendiente y ordenada en el origen, relacionados con R y C según: m=

1 1 = A R⋅C

b=

R = A

R R⋅C

Deshaciendo los cambios hechos, los coeficientes de el ecuación (5.15) son: 1 b⋅m b R= m

C=

Con los valores así calculados de R y C, y los valores de x, se calculan unas yc, que deberán estar muy próximos a las y de partida. De la bibliografía especializada se obtuvieron los datos de equilibrio isobárico, 760 mm de Hg. líquido-vapor de Benceno-nButanol, que se utilizan en el programa Equilibrio.m, en donde se aplican los cálculos expuestos. %----------------- Equilibrio.m ---------------datos=[.0064 .0201 .0087 .0285 .0098 .0318 .0103 .0495 .0111 .0568 .0152 .0666 .0198 .0778 .0208 .0939 .0250 .0984 .0331 .1295 .0488 .1898


77

.0553 .2170 .0568 .2199 .0620 .2356 .0753 .2824 .0864 .3152 .1120 .3840 .1209 .4082 .1553 .4881 .1779 .5361 .1992 .5693 .2390 .6421 .2718 .6712 .2852 .6836 .3440 .7493 .3978 .7874 .4528 .8161 .5295 .8601 .5753 .8689]; %----------------------------------------------x=datos(:,1); y=datos(:,2); D=y–x; U=x./(1–x); V=sqrt(U./D); est=polyfit(U,V,1); m=est(1); b=est(2); subplot(221), plot(x,y,'o-') subplot(222), plot(U,V,'o',U,m*U+b,'-') C=1/(b*m); R=b/m;

V=(U/(y–x))0,5

2

m= 1,152 b= 0,55015

R= 0,47756 C= 1,5779

1,5

1

0,5

0,2

0,4

0,6 0,8 U=x/(1–x)

1

1,2

0,8

y

0,6 0,4 0,2 0,05

0,1

0,15

0,2

0,25

0,3 x

0,35

0,4

0,45

0,5

0,55

Figura 5.12. Composición benceno en la fase vapor frente a composición de benceno en la fase líquida.

78


yc=x+C*R*x.*(1-x)./(x+R*(1-x)).ˆ2; subplot(212),plot(x,y,'o',x,yc,'-') %-----------------------------------------------

En la Figura 5.12, en la parte superior, se muestran los datos sometidos a las transformaciones indicadas, y la recta ajustada a estos datos, haciendo uso de polyfit, observándose un buen grado de concordancia. En el gráfico inferior se representan los datos (x,y) de partida, representados por «o», y en línea continua se muestra la ecuación (5.15) a partir de las x y con los datos calculados R y C. 5.7. ANÁLIS ESPECTRAL El procesado de señales digitales para análisis espectral y series temporales, se facilita con la función fft, fast fourier transform. En el análisis de señales con perturbaciones, el problema más frecuente es la determinación de las frecuencias. En el siguiente ejemplo se crea una señal compuesta, a la que se le añade un ruido aleatorio de cuantía análogo a la señal. El período de muestreo se toma en milisegundos, creando una base de tiempo de 0 a 0,5 segundos. >> t=[0:0.001:0.5]'; >> size(t) ans = 501

1

Se crea una señal sinusoidal compuesta de 50, 100 y 200 Hz. >> x=sin(2*pi*50*t)+sin(2*pi*100*t)+sin(2*pi*200*t);

Esta señal se distorsiona con un ruido aleatorio de 0 a 3. y=x+3*rand(length(t),1);

La representación gráfica de y frente a t conduce a una representación en la que la periodicidad está oculta. El objetivo de este tratamiento es encontrar las frecuencias de esta señal distorsionada, para ello se procede con estos datos distorsionados, como datos de partida, (t,Y) para calcular las frecuencias. A las y a analizar, se les aplica la transformada de Fourier. Y=fft(y);


79

Ahora se calcula la potencia espectral o energía, mediante: P=Y.*conj(Y); >> size(P) ans = 5001

1

La frecuencia se forma mediante la instrucción: >> f=1000*(1:256)/512; >> size(f) ans = 1 256 >> f=f'; >> size(f)

1

y=x+3 ⋅ rand(length(t),1)

x=sin(2 ⋅ pi ⋅ 50 ⋅ t)+sin(2 ⋅ pi ⋅ 100 ⋅ t)+sin(2 ⋅ pi ⋅ 200 ⋅ t)

ans = 256

2 1 0 –1 –2 0

Potencia Espectral

7

0,05 0,1 Tiempo (seg.)

5 4 3 2 1 0 –1 –2 0

0,05 0,1 Tiempo (seg.)

x 10 4

6 5 4 3 2 1 0

50

100

150 Frecuencia

200

250

Figura 5.13. Análisis de frecuencias.

Las representaciones gráficas se consiguen con: >> subplot(221), plot(t(1:50),x(1:50)) >> ylabel('x=sin(2·pi·50·t)+sin(2·pi·100·t)+sin(2·pi·200·t)')

80

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


xlabel('Tiempo') subplot(222), plot(t(l:50),y(1:50)) ylabel('y=x+3·rand') xlabel('Tiempo') subplot(212), plot(f,P(1:256)) axis([ 0 500 0 70000]) ylabel('Potencia espectral') xlabel('Frecuencia')

En la Figura 5.13 se muestran las variables relacionadas en los cálculos expuestos, destacando claramente las potencias espectrales en las frecuencias buscadas sobre el rizado de la línea de base. 5.8. EVITANDO LA DIVISIÓN POR CERO Y REBOSE La exactitud con que funciona MatLab es eps, cuyo valor es 2,2204 · 10–16. Al dividir una expresión por una variable x, que en un proceso de cálculo tomase el valor 0, provoca un warning Divide by zero. Este inconveniente se obvia mediante: x=x+(x==O)*eps;

Este procedimiento únicamente cambia el valor de x, cuando su valor es 0, asignándole el valor de eps. Los números más grande y más pequeño, en valor absoluto, que utiliza MatLab son realmax y realmin, de valores 1,7977 · 10308 y 2,2251 · 10–308, respectivamente. Para evitar que una variable x, no exceda de un valor límite dado de saturación o rebose, ± X, se procede según: x=x.*(abs(x) < X)+X.*(abs(x)>X);

6 Simulink

6.1. INTRODUCCIÓN A SIMULINK

Figura 6.1. Ventana de paquetes de librerías de Simulink. 81

82


Simulink es un software para simulación, análisis y modelado de sistemas dinámicos, que acompaña a MatLab, en forma de toolbox. Su interface gráfica permite «ver» los modelos.

Figura 6.2. Librerías de bloques de operadores continuos.

Simulink soporta sistemas lineales y no lineales, continuos, discretos e híbridos. Simulink viene con muchos ejemplos a modo de demos para facilitar su uso.

Figura 6.3. Librerías de bloques no continuos.

En la Figura 6.1 se presentan los menús de las libreras que suministran conjuntos de bloques operacionales, cuyos contenidos se muestran en las Figuras 6.2, a 6.8, con los que se compondrán los modelos deseados.

SIMULINK

83

Figura 6.4. Bloques generadores de señales.

La notación matemática de los bloques es la utilizada en los tratados de Control Automático, en los que se utiliza la variable compleja s, que tiene, entre otras, la propiedad de que multiplicar por s representa derivar, y dividir por s representa integrar.

Figura 6.5. Librerías de funciones y tablas.

En los diagramas de bloques lo que sale es el producto de lo que entra por el contenido del bloque.

84


Figura 6.6. Librerías de bloques matemáticos.

Figura 6.7. Librerías de señales y sistemas.

SIMULINK

85

Figura 6.8. Librerías de salidas.

6.2. CONSTRUCCIÓN DE UN MODELO MUY SENCILLO A continuación se describe paso a paso cómo se elabora el esquema Simulink del modelo más sencillo, la integración de 2 × t, con condición inicial 0, cuyo resultado es t2. 1. Desde la ventana de comandos de MatLab se ejecuta simulink, o se activa pulsando sobre su icono en la barra de tareas, obteniéndose la ventana de librerías de Simulink 2. De la ventana de librerías se consigue la ventana de construcción del modelo, pulsando File, New, Model, con lo que se obtiene una nueva ventana vacía, para construir ahí el modelo deseado. 3. Se despliega Sources de la ventana de librerías. Se arrastra o copia Clock a la ventana del modelo. Clock suministrará t, tiempo continuo. 4. Se abre la ventana Math de la ventana de librerías, y se copia Gain, ganancia, a la ventana del modelo. Gain se coloca delante de Clock. Para conectar Gain a Clock, se posiciona el ratón en el ángulo de salida de Clock, manteniendo pulsado el botón izquierdo del ratón, se desplaza a la entrada de Gain. De este modo Clock y Gain quedan unidos. Pulsando dos veces sobre Gain se seleciona su valor a 2. La salida de Gain es su entrada multiplicada por su valor seleccionado. La salida de Clock suministra t, y la salida de Gain suministra 2 × t. 5. Pueden ponerse comentarios informativos donde se desee haciendo doble click en el lugar elegido para pasar a modo texto. Para cambiar tamaños y tipos de letra, se coloca el puntero del ratón encima del texto, al pulsar el botón derecho se accede a Propiedades del texto.

86


Figura 6.9. Esquema de la integración de 2 · t.

6. Desplegando Continuous de la ventana de librerías, se copia Integrator a la ventana del modelo. Se une la salida de Gain con la entrada de Integrador. Pulsando dos veces sobre Integrador, se abre una ventana de propiedades del bloque integrator, en este caso únicamente se selecciona la condición inicial a 0, valor por defecto. 7. La salida del integrator, es la integral de su entrada, por tanto, con él se logra el objetivo buscado, t2. 8. Para ver la integración conseguida, es necesario llevar esa salida a un visualizador, X Y graph, que se halla en la librería Sinks, desde donde se copia a la ventana del modelo, cuyas entradas corresponden a las salidas del reloj, X, y del integrador, Y. 9. Los métodos y parámetros de integración, son seleccionables desde la ventana del modelo en Simulation parameters, del menú Simulation. 10. Para ejecutar la simulación se activa Start del menú Simulation de la ventana del modelo. Para ver el resultado gráfico se pulsa dos veces en X Y graph. 11. Para llevar variables del modelo de simulación al entorno MatLab, se utiliza el módulo To workspace de la librería Sinks. Este módulo se conecta en la salida del bloque del que se desea guardar sus datos; pulsando dos veces sobre este módulo, se accede a la asignación del nombre de la variable y al formato en que se desean guardar los datos. 6.3. SOLUCIÓN SIMULINK DE UNA ECUACIÓN DIFERENCIAL Para representar un modelo, ecuaciones diferenciales, se abre el espacio de representación mediante Nuevo, del menú Archivo de la ventana de librerías de Simulink, hacia donde se arrastrarán los iconos operacionales desde las librerías correspondientes, para componer el modelo mediante las conexiones y relleno con los parámetros pertinentes.

SIMULINK

87

En la Figura 6.10 se representa el diagrama de bloques de la ecuación (5.2); y su solución gráfica, en la Figura 6.11.

Figura 6.10. Esquema analógico de la resolución de la ecuación diferencial de segundo grado.

Figura 6.11. Resultado gráfico de la integración de la ecuación diferencial de segundo grado.

88


Como la ecuación diferencial es de segundo grado, es necesario integrar dos veces. El primer integrador tiene por entrada la derivada segunda, y su salida será la derivada primera, que será la entrada del segundo integrador, para suministrar la primitiva en función del tiempo. La entrada del primer integrador es el segundo miembro de la ecuación diferencial expresada como: d2y dy = –3 – 6 y 2 dt dt

6.4. SIMULACIÓN DINÁMICA DE UN ECUALIZADOR Al ecualizador, tanque de regulación de caudal y amortiguador de oscilaciones de concentración de influentes, en procesos de aguas, llega una corriente de caudal y concentración de sustancias variables, y se pretende determinar el volumen variable ocupado por el líquido en el ecualizador cuando se extrae un caudal de valor medio de las últimas 24 horas, media móvil, a la vez que determinar la concentración de las sustancias, sustrato, variables en la salida.

Qi(t) Ci(t)

Volumen de seguridad

Variación de volumen

Volumen inicial

V(t)

C(t)

Q C(t)

Figura 6.12. Esquema y variables de un ecualizador.

SIMULINK

89

La variación de volumen de líquido en el ecualizador está determinado, por lo que entra y por lo que sale en un instante: dV = Qi – Q dt

(6.1)

Si en el tanque ecualizador no hay reacción química alguna, conservación de componentes, la concentración de cada componente que se mezcla con el contenido del tanque de volumen variable sigue la ecuación: dC Qi dV ⎞ C = ⋅ Ci – ⎛ Q + ⋅ ⎝ dt V dt ⎠ V

(6.2)

La variación de volumen, puede calcularse analíticamente, si se considera que el caudal de entrada es sinusoidal, tal como: 2 ⋅π ⎞ Qi (t ) = Q + Qv ⋅ sin ⎛ ⋅t ⎝ T ⎠

(6.3)

El volumen del líquido en el ecualizador estará dado por: V = V (0) + Qv ⋅

t

⎛ 2 ⋅ π ⋅ t⎞ T ⎠

∫ sin ⎝ 0

(6.4)

Cuyo resultado es: V = V (0) +

2 ⋅π ⎞ ⎤ Qv ⋅ T ⎡ ⋅ ⎢1 – cos ⎛ ⋅t ⎝ T ⎠ ⎥⎦ 2 ⋅π ⎣

(6.5)

Ecuaciones en las que V(0), representa el contenido del tanque en el instante inicial; T, el período de una oscilación completa; y Qv, representa la amplitud de la oscilación del caudal y C la concentración. En la Figura 6.13 se muestra el esquema de cálculo Simulink para resolver la ecuación (6.1). En el bloque denonimado Qi, se suministran los datos del caudal de entrada frente al tiempo, correspondientes a un ciclo, mediante la Repeating table, que pide una tabla de valores, tiempo-valor, correspondientes a un ciclo completo, en este caso: [0 2 4 6 8 10 12 14 16 18 20 22] [211 168 155 162 220 330 460 533 475 370 290 250]

90


Figura 6.13. Esquema Simulink para el cálculo de la evolución del volumen contenido en un ecualizador con caudal variable.

El bloque Mux permite conducir por una línea varias señales. En el ejemplo, Mux se utiliza para realizar operaciones en los bloques f(u), y para representar varias señales en un gráfico. 2500

Caudales, (m3/h), Volumen (m3)

Qi Q Volumen 2000

1500

1000

500

0

0

20

40

60

80

100

120

Tiempo, (Horas)

Figura 6.14. Resultado gráfico del cálculo de la evolución del volumen contenido en un ecualizador con caudal variable.

En el bloque Qm(t-24), se realiza la operación u(l)/u(2), la señal u(1) corresponde a la integral del caudal retrasado, y u(2) corresponde al tiempo retrasado.

SIMULINK

91

El bloque Qi(t-24) suministra el caudal de entrada retrasado en 24 horas; durante las primeras 24 horas, suministra un valor de designación inicial. El bloque Vol, integra el caudal retrasado, para obtener el volumen aportado por el caudal retrasado, que dividido entre t – 24 produce el caudal medio retrasado en las últimas 24 horas. El bloque denominado Limitador, tiene por objeto limitar posibles valores de rebose, por valores iniciales muy altos o muy bajos. El bloque Volumen integra la diferencia del caudal de entrada menos el caudal medio retrasado.

Figura 6.15. Esquema Simulink para el cálculo del volumen y composición de un ecualizador, con alimentación de caudal y concentración variables.

Los resultados de la simulación de los caudales y del volumen ocupado en el tanque se presentan en la Figura 6.14, en la que se observa que el caudal medio de las últimas 24 horas, caudal de salida, en el comienzo de la simulación, tiene un valor constante hasta llegar al tiempo de 24 lıoıas, seguido de un valor constante, límite superior, y a continuación oscila, con tendencia hacia un valor constante, debido a que el caudal de entrada se repite de la misma manera. El volumen del líquido sufre una primera aportación muy grande, debido a que en el comienzo, el caudal purgado está en un valor mínimo, observándose que tiende a una oscilación constante. En la Figura 6.15 se muestra el esquema de cálculo Simulink para resolver las ecuaciones (6.1) y (6.2); cabe destacar únicamente que el integrador C va provisto de limitadores, a fin de evitar valores anómalos iniciales que enlentecen el cálculo en alcanzar el estado estacionario. El bloque Ci suministra los datos de la concentración de entrada de un supuesto componente químico o sustrato: [0 2 4 6 8 10 12 14 16 18 20 22] [150 130 89 109 165 195 398 413 364 270 215 174]

92


La Figura 6.16 muestra la reducción de la variabilidad de la concentración de salida del ecualizador, línea continua, frente a la concentración de entrada, línea discontinua. 450 C Ci

400

Conc. Entrada/salida (mg/l)

350

300

250

200

150

100

50

0

20

40

60

80 100 120 Tiempo (Horas)

140

160

180

200

Figura 6.16. Concentración de salida del ecualizador, línea continua, frente a concentración de entrada.

Las Estaciones Depuradoras de Aguas Residuales que se vienen instalando desde hace unos años, incluyen el tanque de homogeneización o ecualización, que si es de tamaño suficiente y con buen control de operación, ahorrará energía, diluirá posibles tóxicos o permitirá su «by-pass» a esta zona segura, evitará los problemas de las variaciones bruscas de caudal y corıtribuirá a la estabilidad del proceso biológico posterior.

Bibliografía

[1] Brian D. Hanhn. Essential MatLab for scientists and engineers. Buterworth-Heinemann. (2002). [2] Brian D. Hunt. A guide to MatLab for beginners and experienced users. Cambridge Univ. Press. (2001). [3] Duane Hanselman, Bruce R. Littlefield. Mastering MatLab 6. The MatLab curriculum series. Prentice Hall. (2000). [4] James B. Dabney; Thomas L. Harman. Mastering simulink 4 S/E. The MatLab curriculum series. Prentice Hall. (2001). [5] Edward B. Magrab. Engineers guide to MatLab. Prentice Hall. (2000). [6] Rudra Pratrap. Getting starting with MatLab, Version 6: Quick introduction. Oxford University Press. (2001). [7] Stephen J. Chapman. MatLab programming for engineers. Brooks-Cole Pub. Co. (2001) [8] Gil Rodríguez M. Cálculos avanzados en procesos de descontaminación de aguas. Consejo Superior de Investigaciones Científicas. (2003).

93

Índice alfabético

..., 21 LATEX, 19, 45 TEXtos en gráficos, 45

continue, 40 conv, 16 date, 20 datenum, 20 datestr, 20

Ajuste de datos, 56 Ajuste lineal, 72 angle, 12 ans, 12 area, 48 ascii, 22, 23 atan2, 13 axes, 44 axis, 43, 44

datetick, 44 datevec, 20 deconv, 16 demo, 1 Derivadas, 27, 67 de polinomios, 17 det, 15, 26 diary, 23 diff, 27, 32 double, 28 dsolve, 29, 32

bar, 41 bar3, 46 Bazin, ecuación, 63 break, 40 case, 40 ceil, 22 clear, 5 clf, 44, 54 clock, 20 Color, 47, 48 Complejos, números, 12 conj, 12, 79

else, 36, 37 elseif, 37 eps, 80 errorbar, 41 etime, 21 ezplot, 27 factor, 27 95

96


Factorial, 11 feval, 10 fft, 78 fix, 11 floor, 22 fminbnd, 52 fmins, 51 fminsearch, 53 fontsize, 45, 48 for-end, 10, 38 format, 5, 6 fplot, 41 Fracciones simples, desarrollo, 17 Funciones del tiempo, 20 function, 8 fzero, 62 global, 10, 69 grid, 8 gtext, 44 help, 1, 9 hold on/off, 44 if else end, 36 ilaplace, 30 imag, 12 inline, 11 input, 35 int, 27 integrador, 86 Integrales, 27 de polinomios, 18 interp1, 18 Interpolación cubic, 18 polinomial, 18 spline, 18, 71 inv, 15 invlaplace, 33 keywords, 4

laplace, 30, 32, 33 LaPrint, 46 leftarrow, 48 legend, 44 length, 55, 72 limit, 31 line, 44, 54 LineWidth, 44, 48 linspace, 2, 13 load, 21 loglog, 41 logspace. 13 long, 6 lookfor, 9 lsqcurvefit, 55 magic, 26 Maple V, 25, 31, 32 matlabpath, 19 Matriz, 13 max, 8 mpower, 15 mux, 90 now, 20 num2str, 57 ode113, 57 ode15s, 57 ode23, 57 ode23s, 57 ode23t, 57 ode23tb, 57 ode45, 57 odegil4, 57, 60 optimset, 55 Orden de ejecución, 19 otherwise, 40 path, 19 pie, 41 plot, 8, 41 plot3, 41

ÍNDICE ALFABÉTICO

plotmatrix, 41, 43 plotyy, 41 polar, 43 Polinomios, 16 poly, 16 polyder, 17 polyfit, 56 polyint, 18 polyval, 16 power, 15 pretty, 26 quad, 51 quadl, 51 rand, 21 rat, 5 real, 12 realmax, 80 realmin, 80 recursion, 11 reshape, 14 residue, 17 Retrasos, 67 return, 9 root, 16 round, 29, 70 Runge-Kuta, 57, 59 save, 21 scripts, 9 semilogx, 41 semilogy, 41 short, 5

Simbólico, cálculo, 25 simplify, 27 size, 4 solve, 28 spline, 71 sqrt, 6, 7, 21 stairs, 41 stem, 41 subplot, 43, 44 subs, 58 sum, 55 switch, 40 sym, 25 syms, 26 taylor, 19, 31 text, 44, 54, 55 title, 44 Transposición, 7 variables, 3 Vector, 13 Vertedero, 63 what, 21 which, 19 while-end, 38 who, 2, 3 whos, 2, 3 workspace, 2 xlabel, 44 ylabel, 44

97

Guia Rapida Para Matlab y Simulink

Recommend Documents