CURSO PARA INGENIEROS.pdf

Facultad de Ingeniería Química y Metalurgia

CURSO PARA INGENIEROS: Introducción a MATLAB y SIMULINK Por:

Ing° José S. Orbegoso López

CONTENIDOS:

1. Generalidades. 2. Manipulación de matrices. 3. Gráficas en 2D y 3D. 4. Programación. 5. Métodos numéricos. 6. Matemáticas simbólicas. 7. Simulink.

Huacho, Abril de 2013

MANUAL DE MATLAB BASICO I MATLAB VER 7.12.0.635 (R2011a)

ING. JOSE SAUL ORBEGOSO LOPEZ

1. Generalidades Desarrollado a finales de los setenta en las Universidades de New Mexico y Stanford. En la actualidad ampliamente difundido en diversas Universidades e Industrias. Núcleo básico del programa: Cálculo numérico basado en matrices. Matlab= Matrix Laboratory. Módulos adicionales. Simulink, MatLab Symbolic Toolbox (núcleo del programa Maple), etc. La versión última es R2013a.

Operadores

2



Tipos de datos: MATLAB es un programa preparado para trabajar con vectores y matrices. Como caso particular también trabaja con variables escalares (matrices de dimensión 1). MATLAB trabaja siempre en doble precisión, es decir guardando cada dato en 8 bytes, con unas 15 cifras decimales exactas. Ya se verá más adelante que también puede trabajar con cadenas de caracteres (strings) y, desde la versión 5.0, también con otros tipos de datos: Matrices de más dos dimensiones, matrices dispersas, vectores y matrices de celdas, estructuras y clases y objetos.

Otras formas de definir matrices: Existen en MATLAB varias funciones orientadas a definir con gran facilidad matrices de tipos particulares. Algunas de estas funciones son las siguientes: eye (4) forma la matriz unidad de tamaño (4x4) zeros(3,5) forma una matriz de ceros de tamaño (3x5) zeros(4) ídem de tamaño (4x4) ones(3)forma una matriz de unos de tamaño (3x3) ones(2,4) ídem de tamaño (2x4) linspace(x1,x2,n) genera un vector con n valores igualmente espaciados entre x1 y x2logspace(d1,d2,n) genera un vector con n valores espaciados logarítmicamente entre10^d1 y 10^d2. Si d2 es pi6, los puntos se generan entre 10^d1y pi rand(3)forma una matriz de números aleatorios entre 0 y 1,con distribución uniforme, de tamaño (3x3) rand(2,5) ídem de tamaño (2x5) randn(4)forma una matriz de números aleatorios de tamaño(4x4), con distribución normal, de valor medio 0 y varianza magic(4)crea una matriz (4x4) con los números 1, 2, ... 4*4, con la propiedad de que todas las filas y columnas suman lo mismo hilb(5)crea una matriz de Hilbert de tamaño (5x5). La matriz de Hilbert es una matriz cuyos elementos (i, j) responden a la expresión (1/(i+j-1)). Esta es una matriz especialmente difícil de manejar por los grandes errores numéricos a los que conduce.

Formación de una matriz a partir de otras A continuación se describen algunas de las funciones que crean una nueva matriz a partir de otra o de otras, comenzando por dos funciones auxiliares: [m,n]=size(A) devuelve el número de filas y de columnas de la matriz A. 3



Si la matriz es cuadrada basta recoger el primer valor de retorno N=length(x)calcula el número de elementos de un vector x zeros(size(A))forma una matriz de ceros del mismo tamaño que una matriz A previamente creada. ones(size(A))ídem con unos. A=diag(x)forma una matriz diagonal A cuyos elementos diagonales son los elementos de un vector ya existente x x=diag(A)forma un vector x a partir de los elementos de la diagonal de una matriz ya existente A diag(diag(A))crea una matriz diagonal a partir de la diagonal de la matriz A. blkdiag(A,B)crea una matriz diagonal de submatrices a partir de las matrices que se le pasan como argumentos. triu(A)forma una matriz triangular superior a partir de una matriz A (no tiene por qué ser cuadrada). tril(A)ídem con una matriz triangular inferior

OPERADOR DOS PUNTOS (:) Defínase un vector x con el siguiente comando: >> x=1:10 x= 1

2 3 4 5 6 7 8 9 10

En cierta forma se podría decir que el operador (:) representa un rango: en este caso, los números enteros entre el 1 y el 10. Por defecto el incremento es 1, pero este operador puede también utilizarse con otros valores enteros y reales, positivos o negativos. En este caso el incremento va entre el valor inferior y el superior. >> x=1:2:10 x= 1

3 5 7 9

>> x=10:-1:1 x= 10 9 8 7 6 5 4 3 2 1 4



OPERADORES RELACIONALES

OPERADORES LOGICOS

FUNCIONES MATEMÁTICAS ELEMENTALES QUE OPERAN DE MODO ESCALAR

5



6



FUNCIONES QUE ACTÚAN SOBRE LOS VECTORES

FUNCIONES QUE ACTUAN SOBRE MATRICES

7



Tabla de códigos ASCII - Formato de caracteres estándares ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0 1 2 3 4 5 6 7 8 9 A B C D E F

NUL SOH STX ETX EOT ENQ ACK BEL BS TAB LF VT FF CR SO SI

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

10 11 12 13 14 15 16 17 18 19 1A 1B 1C 1D 1E 1F

DLE DC1 DC2 DC3 DC4 NAK SYN ETB CAN EM SUB ESC FS GS RS US

32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

20 21 22 23 24 25 26 27 28 29 2A 2B 2C 2D 2E 2F

(espacio) ! " # $ % & ' ( ) * + , . /

48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

30 31 32 33 34 35 36 37 38 39 3A 3B 3C 3D 3E 3F

0 1 2 3 4 5 6 7 8 9 : ; < = > ?

ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

ASCII

Hex

Símbolo

64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

40 41 42 43 44 45 46 47 48 49 4A 4B 4C 4D 4E 4F

@ A B C D E F G H I J K L M N O

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

50 51 52 53 54 55 56 57 58 59 5A 5B 5C 5D 5E 5F

P Q R S T U V W X Y Z [ \ ] ^ _

96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111

60 61 62 63 64 65 66 67 68 69 6A 6B 6C 6D 6E 6F

` a b c d e f g h i j k l m n o

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

70 71 72 73 74 75 76 77 78 79 7A 7B 7C 7D 7E 7F

p q r s t u v w x y z { | } ~ •

8



En caso de que el teclado se desconfigure se puede recurrir a ésta tabla, digitando: ALT + NUMERO DE LA COLUMNA ASCII, para obtener el carácter necesario. Ejemplo: ALT + 92 va a dar \ ALT + 91 va a dar [ ALT + 93 va a dar ] ALT + 94 va a dar ^ ALT + 58 va a dar : ALT + 59 va a dar ; Etc.

GUIA DE LA CLASE1 Ambiente de trabajo de MatLab. Workspace, Windows y Help.

9



Fig 1: Ventana Inicial de Matlab 7.12.0.635(R2011a) En la ventana inicial se sugieren ya algunos comandos para el usuario inexperto que quiere echar un vistazo a la aplicación. En dicha ventana aparece también el prompt característico de MATLAB (»). Esto quiere decir que el programa está preparado para recibir instrucciones. Puede hacerse que aparezca un saludo inicial personalizado por medio de un fichero de comandos personal que se ejecuta cada vez que se entra en el programa. Una vez que el programa comienza, se visualizan tres pequeñas ventanas: Command, Wokspace y Command History. Estas son 3 de las ocho diferentes ventanas en Matlab. El listado de las ventanas y sus propósitos se muestran en la tabla 1.

Tabla 1 : Ventanas de Matlab Ventana

Propósitos

Command Window

Ventana principal, variables enteras, corre programas.

Figure Window

Contiene información para controles gráficos.

Editor Window

Crea y elimina escritos y archivos de funciones.

Help Window

Proporciona información de ayuda

Launch Pad Window

Provee acceso a herramientas, demostraciones y documentación. 10



Command History Window

Introduce los controles principales en la ventana de control .

Workspace Window

Proporciona información de las variables que son usadas.

Current Directory Window

Muestra los archivos actuales en el directorio.

Ventana Command: Esta ventana es la ventana principal y se abre cuando comienza Matlab. Es conveniente tener la ventana command como la única ventana visible y esto puede hacerse cerrando todas las otras ventanas. Ventana Figure: La ventana figure se abre automáticamente cuando ejecutamos los comandos gráficos y contienen gráficos creados por estos comandos. Un ejemplo de la ventana Figure es mostrada en la figura 2.

{ Puede usted ingresar lo siguiente, utilizando una hoja de editor: >> x = [0:0.1:16]; >> y = x.^ (1.4).*cos(x) >> plot(x,y) >> title('Gráfica de la función y = x.^1^.^4.*cos(x)') }

FIGURA 2 : GRÁFICO DE LA FUNCIÓN COSENO(X)

11



Ventana Editor: La ventana editor es usada para escribir y editar, crear y modificar programas, tanto como ejecutarlos paso a paso para ver si contienen errores (proceso de depuración). El editor muestra con diferentes colores los diferentes tipos o elementos de los comandos (en verde, para los comentarios; en rojo, las cadenas de caracteres, etc.). El editor se preocupa de que las comillas o paréntesis que se abren, no se queden sin cerrar. Figura 3

Figura 3: Ventana Editor Ventana Help MATLAB 7.12.0.635 (R2011a) dispone de un excelente Help con el que se puede encontrar la información que se desee. La Figura 4 muestra las distintas opciones que aparecen en el menú Help. 1. Help Window, Se abre la ventana de la Figura 4, en la que se puede buscar ayuda sobre la función o el concepto que se desee. 2. Help Tips. Ofrece ideas prácticas para utilizar la ayuda. 3. Help Desk. Se abre un browser de Internet (Netscape Communicator, en el caso de la Figura 5) que permite acceder a toda la información sobre MATLAB en formato HTML. Esta información es equivalente a los manuales impresos del programa. Desde la parte inferior de esta página, mediante el enlace Online Manuals (in PDF) se puede acceder a la versión *.pdf (PortableDocument Format) de los manuales de MATLAB.

12


FIGURA 4:


VENTANA HELP

13



Este formato es menos adecuado para consultar que el HTML, pero mucho más adecuado para imprimir y revisar luego sobre papel. El formato *.pdf requiere del programa gratuito Adobe Acrobat Reader 5.0 o una versión superior. 4. Examples and Demos. Se abre una ventana que da acceso a un buen número de ejemplos resueltos con MATLAB, cuyos resultados se presentan gráficamente de diversas formas. Es bastante interesante recorrer estos ejemplos para hacerse idea de las posibilidades del programa.

FIG 5: VENTANA DEMOS DE HELP WINDOW.

FIG 6 : VENTANA TOOLBOX DE HELP

DESK

14



Además, se puede también recurrir al Help desde la línea de comandos. Se aconseja hacer prácticas al respecto. Por ejemplo, obsérvese la respuesta a los siguientes usos del comando help: » help » help lang El comando helpwin seguido de un nombre de comando muestra la información correspondiente a ese comando en la ventana Help Window (ver Figura 5), incluyendo también comandos similares sobre los que se ofrece ayuda. El comando doc seguido de un nombre de comando muestra la información correspondiente a ese comando a través de Netscape Navigator o Internet Explorer, en formato HTML. El Debugger es un programa enormemente útil para detectar y corregir errores, que hay que conocer muy bien. Para ello lo mejor es practicar. De hecho, al realizar operaciones con el Debugger gráfico van a apareciendo las correspondientes instrucciones en la línea de comandos. VISITE:

VENTANA WORKSPACE

El espacio de trabajo consiste en el sistema de variables acumuladas durante una sesión de usar el software de MATLAB y almacenadas en memoria. Usted agrega variables al espacio de trabajo usando las funciones, M-archivos corrientes, y cargando espacios de trabajo ahorrados. Por ejemplo, si usted funciona estas declaraciones, A = magic(4) R = randn(3,4,5) El espacio de trabajo incluye dos variables, A y R. Usted puede realizar workspace y operaciones relacionadas usando el browser del workspace. Cuando están disponibles, las funciones equivalentes se documentan con cada característica del browser del espacio de trabajo. Si usted tiene una conexión activa del Internet, usted puede mirar el Workspace browser video demo para una descripción de la funcionalidad principal:

15



Matemáticas escalares. Funciones matemáticas básicas. 2. Operaciones con Vectores y Matrices >> A=[1 2 3; 4 5 6; 7 8 9] La respuesta del programa es la siguiente: A=

123 456 789 Una sencilla operación con A es hallar su matriz traspuesta. En MATLAB el apóstrofe (') es el símbolo de trasposición matricial. Para calcular A' (traspuesta de A) basta teclear lo siguiente:

>> A' ans

= 14 7 258 369

La variable ans puede ser utilizada como operando en la siguiente expresión que se introduzca. También podría haberse asignado el resultado a otra matriz llamada B:

>> B = A' 16



B = 147 258 369 Ahora ya están definidas las matrices A y B, y es posible seguir operando con ellas. Por ejemplo, se puede hacer el producto B*A (deberá resultar una matriz simétrica): >> B*A ans = 66 78 90 78 93 108 90 108 126 En MATLAB se accede a los elementos de un vector poniendo el índice entre paréntesis (por ejemplo x(3) ó x(i)). Los elementos de las matrices se acceden poniendo los dos índices entre paréntesis, separados por una coma (por ejemplo A(1,2) ó A(i,j)). Las matrices se almacenan por columnas (aunque se introduzcan por filas, como se ha dicho antes), y teniendo en cuenta esto puede accederse a cualquier elemento de una matriz con un sólo subíndice. Por ejemplo, si A es una matriz (3x3) se obtiene el mismo valor escribiendo A(1,2) que escribiendo A(4). Invertir una matriz es tan fácil como transponerla. A continuación se va a definir una nueva matriz A – no singular-en la forma: >> A = [1 4 -3; 2 1 5; -2 5 3] A= 1 4 -3 2 1 5 -2 5 3 Ahora se va a calcular la inversa de A y el resultado se asignará a B. Para ello basta hacer uso de la función inv ( ) >> B = inv(A) B= 0.1803 0.1311 -0.0984

0.2213 0.0246 0.1066

-0.1885 0.0902 0.0574

Para comprobar que este resultado es correcto basta con multiplicar A por B >> B*A ans = 17


1.0000 0.0000 0.0000

0.0000 1.0000 0.0000


0.0000 0.0000 1.0000

De forma análoga a las matrices, es posible definir un vector fila x en la forma siguiente (silos tres números están separados por blancos o comas, el resultado será un vector fila): >> x=[10 20 30] % vector fila x= 10 20 30 MATLAB considera comentarios todo lo que va desde el carácter tanto por ciento (%) hasta el final de la línea. Si los números están separados por intros o puntos y coma (;) se obtendrá un vector columna: >> y=[11; 12; 13] % vector columna y= 11 12 13 MATLAB tiene en cuenta la diferencia entre vectores fila y vectores columna. Por ejemplo, si se intenta sumar los vectores x e y se obtendrá el siguiente mensaje de error: >> x+y ??? Error using ==> + Matrix dimensions must agree. (Deben estar de acuerdo las dimensiones de la matriz). Estas dificultades desaparecen si se suma x con el vector transpuesto de y: >> x + y' ans = 21 32 43

EJERCICIOS UTILIZANDO MATLAB 1. Transferencia de Calor: Un objeto con una temperatura inicial de T 0 es colocado al tiempo t = 0 dentro de una cámara que tiene una temperatura constante de Ts experimentará un cambio de temperatura de acuerdo a la ecuación: T = Ts + (T0 – Ts) e-kt Donde:T = Temperatura del objeto al tiempo. 18



k = constante. Una soda a una temperatura de 120ºF es ubicada dentro de un refrigerador donde la temperatura es 38ºF. Determine la temperatura alcanzada después de 3 horas. Asumir k=0.45. Primero define todas las variables y entonces calcula la temperatura usando un comando Matlab. Solución >> Ts= 38;

T0 = 120;

k = 0.45;

t = 3;

>>T = round (Ts +(T0 – Ts) *exp(-k*t)) T = 59 2. Identidad Trigonométrica: Una identidad trigonométrica está dado por:  x  tan x  sin x cos2    2 tan x 2

Verificar que la identidad es correcta calculando en cada lado de la ecuación y sustituyendo x = ∏/5. Solución >> x = pi/5; >>LHS = cos(x/2)^2 LHS = 0.9045 >>RHS = (tan(x) +sin(x))/(2*tan(x)) RHS = 0.9045

3. La presión de vapor de agua Escribir una función¸ en MATLAB, para el cálculo de la presión de vapor del agua. Use una función¸ es decir, s para calcular la presión de vapor de agua para las temperaturas de 30, 60 y 90°C.

Resolución: La presión de vapor de agua, pv [bar] es función de la temperatura T [K] y puede ser calculado por la expresión: 𝑝𝑣 𝑎𝑇𝑟 + 𝑏 𝑇𝑟1,5 + 𝑐 𝑇𝑟3 + 𝑑 𝑇𝑟6 ln = 𝑃𝑐 1 − 𝑇𝑟 Donde: 19



𝑇𝑟 = 1 −

𝑇 𝑇𝐶

Los coeficientes¸ a, b, c, d¸ con valores -7,76451, 1,45838, -2,7758 respectivamente 1,23303. En estas relaciones Tç es una temperatura crítica y Pc la presión crítica del agua. Esta tiene los valores de 647,3 y 221,2K¸bar. De esta relación se puede expresar la presión de vapor:

esta ecuación se utiliza para escribir la función MATLAB¸ necesario. Esta función llama un argumento que tiene en cuenta el valor de la temperatura y valor de retorno actual correspondiente presión de vapor. La función se presenta a continuación: function pv = pvapH2O(T); % Función para calcular la presión de vapor del agua %% Mod de llamar: % PV = PVAPH2O(T) % donde: PV - presión de vapor, [bar] % T - temperatura, [K] Tc = 647.3; % temperatura critica, [K] Pc = 221.2; % presión critica, [bar] a = -7.76451; b = 1.45838; c = -2.7758; d = -1.23303; Tr = 1 - T./Tc; pv = exp((a*Tr+b*Tr.^(1.5)+c*Tr.^3+d*Tr.^6)./(1-Tr)).*Pc; % end function pvapH2O

Llamado esta función puede obtener un escalar. >> pv = pvaph2o(303.15) pv = 0.0425 O a través de un vector matriz. >> pv = pvaph2o(273.15+[30 60 90]) pv = 0.0425 0.1994 0.7015

OTROS EJERCICIOS: Polinomios:

20



Coeficientes del polinomio completo y ordenado decrecientemente: 𝑃 𝑥 = 𝑥3 + 2 𝑥2 + 3 x + 4 >>A=[1 2 3 4] Comandos: roots: nos da las raíces >>raíces = roots(A) poly: nos da los coeficientes del polinomio que originan esas raíces. >>poly (raíces) >>r=[1 1+j 1-j 2]; >>poly(r) 1

-5 10 -10

4

(𝑥 4 − 5 𝑥 3 + 10 𝑥 2 − 10 𝑥 + 4 = 0)

conv: nos da el product de polinomios >>pp=[1 2 3 4 ]; >>qq=[1 0 0 2 3]; >>producto=conv(pp, qq) deconv: nos da el cociente y el residuo de una división de polinomios >>[co, re]=deconv(qq,pp) co= 1

-2

re= 0

0

1

4 11

polyval: evalá el polinomio P(x) en un valor xo, es decir P(xo) >>A=[1 2 3 4]; >>p=polyval(A,[1 2 3 ]) pl = 26 polyder: nos da los coeficientes de la derivada del polinomio >>D=polyder(A)

21



FUNCIÓN input >>n = input(„teclee el numero de ecuaciones‟); >>nombre=input(„¿Cómo te llamas?‟, „s‟) (en este caso, la respuesta se lee y se devuelve sin evaluar, con lo que se almacena en la cadena nombre). FUNCION disp: Siempre imprime en pantalla un mensaje de texto o el valor de una matriz, pero sin imprimir su nombre >>disp(„EL programa ha terminado‟) >>A=rand(4) >>disp(A) FUNCIONES Definición:

function [lista de valores de retorno]=name (lista de argumentos)

Ejm: Se creará una función que evalúe el numerador, denominador y la división de: P(x,y)=

𝑥 𝑦 + 𝑠𝑒𝑛 𝑥 ∗𝑦 − 𝑥 𝑥 2+ 𝑦 3

function[num, den, c]=division(x,y) num=x.^y + sin(x.*y) – x; den=x.^2 + y.^3; c=num./den; Guardar el archivo como division. m y lo ejecutamos

GUIA DE CLASE 2:

22



Recuerde: Los paréntesis cuadrados “[” y “]” son empleados para matrices, mientras que los paréntesis “(” y “)” se emplean para funciones. El carácter coma “,” para separar elementos [1,2,3,4] y el punto y como “;” para separar filas en las matrices. [ 1 2; 3 4]. El carácter dos puntos “:” indica números consecutivos entre dos valores, [2:5] es equivalente a [2 3 4 5]. Gráficas de funciones MATLAB tiene un gran potencial de herramientas gráficas. Se pueden dibujar los valores de un vector frente a otro (de la misma longitud): >>x=pi*(-1:0.1:1); >>y=x.*sin(x); >>plot(x,y) % Por defecto une los puntos (x(i),y(i)) mediante una poligonal.

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -4

-3

-2

-1

0

1

2

3

4

Como se ve, con pocos puntos la gráfica tiene un aspecto demasiado lineal a trozos. Para "engañar" al ojo, basta tomar más puntos. >>x=pi*(-1:0.01:1); >>y=x.*sin(x); >>plot(x,y)

23



2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -4

-3

-2

-1

0

1

2

3

4

3

4

>>hold on % Mantiene en la ventana gráfica los dibujos anteriores >>plot(x,cos(x)) % Dibuja sobre la gráfica anterior la función cos(x)

2

1.5

1

0.5

0

-0.5

-1 -4

-3

-2

-1

0

1

2

>>hold off % Con esto olvida los dibujos anteriores % y dibuja en una ventana nueva El mismo gráfico de la figura superior puede obtenerse con una única orden, dibujándose cada curva de un color distinto: 24



>>plot(x,cos(x),x,y)

2

1.5

1

0.5

0

-0.5

-1 -4

-3

-2

-1

0

1

2

3

4

Además se puede especificar el color de cada curva, el tipo de línea y el marcador, según la tabla:

Otros comandos muy útiles a la hora de realizar gráficas son: title: le ponemos un título a la gráfica xlabel: le ponemos una etiqueta al eje x ylabel: le ponemos una etiqueta al eje y legend: incluimos una legenda en la gráfica 25



Por ejemplo >> plot(x,cos(x),'ro-',x,y,'bp:')

2

1.5

1

0.5

0

-0.5

-1 -4

-3

-2

-1

0

1

2

3

4

>> xlabel('tiempo') >> ylabel('valor') >> legend('coseno','y') >> title('grafico de ejemplo')

26



grafico de ejemplo 2 coseno y 1.5

valor

1

0.5

0

-0.5

-1 -4

-3

-2

-1

0 tiempo

1

2

3

4

También se pueden realizar gráficas logarítmicas con las instrucciones: semilogx (eje x logarítmico), semilogy (eje x logarítmico) y loglog (ejes x e y logarítmicos). Además de gráficas en 3D con las instrucciones: mesh, surf, plot3... Ejercicio: Representar simultáneamente dos funciones sinusoidales de frecuencia 50 Hz desfasadas 90 º (π/2 rad) de amplitud 10 y 20, respectivamente. Siendo la primera de ellas en color verde a puntos y la segunda en azul continua y con asteriscos. Incorporar además: nombre eje x: “tiempo” nombre eje y: “señal” título gráfica: “onda seno” legenda: “onda1” y “onda 2” Función sinusoidal:

y  A sin  kx   

A: amplitud, k: factor multiplicativo del argumento, que se denomina pulsación  en el caso de que la variable sea independiente del tiempo;  el desfasaje

>>y = 10*sin(x +pi/2); >> plot(x,y) >> hold on >> y = 20*sin(x +pi/2); >> plot(x,y) 27



>>xlabel('tiempo');ylabel('señal'); title('onda seno'); legend('onda1','onda2'); onda seno 20 onda1 onda2

15 10

señal

5 0 -5 -10 -15 -20 -4

-3

-2

-1

0 tiempo

1

2

3

4

POLINOMIOS Y ECUACIONES Para resolver un polinomio, se formula el vector de los coeficientes del mismo, en forma completa y ordenado decrecientemente, siendo el objetivo aquí es la determinación de las raíces. Ejemplo: 𝑃 𝑥 = 𝑥 4 + 2𝑥 2 + 3𝑥 + 4 >>A=[1 2 3 4] A= 1

2

3

4

RAICES: se emplea el comando roots que al digitarlo en la ventana de comando nos da las raíces del polinomio. >> raices=roots(A) raices = -1.6506 28



-0.1747 + 1.5469i -0.1747 - 1.5469i poly: da los coeficientes del polinomio que origina esas raíces. >> poly(raices) ans = 1.0000 2.0000 3.0000 4.0000 >> r=[1 1+j 1-j 2]; >> poly(r) ans = 1

-5 10 -10

4

Producto de polinomios: conocidos los polinomios, ´se pueden multiplicar utilizando el comando conv: >> pp=[1 2 3 4]; >> qq=[1 0 0 2 3]; >> producto=conv(pp,qq) producto = 1

2

3

6

7

12

17

12

El Cociente y Residuo de una división de polinomios: se obtiene empleando el comando deconv: >> [co,re]=deconv(qq,pp) co = 1

-2

re = 0

0

1

4

11

Polyval: evalúa un polinomioP(x) en un valor xo, es decir P(xo) >> A=[1 2 3 4]; >> p=polyval(A,2) 29



p= 26 >> p1=polyval(A,[1 2 3]) p1 = 10 26 58 Polyder: nos da los coeficientes de la derivada del polinomio >> D=polyder(A) D= 3

4

3

FUNCION INPUT Permite imprimir un mensaje en la línea de comandos de MatLab y recuperar como valor de retorno un valor numérico o el resultado de una expresión tecleada por el usuario. Después de imprimir el mensaje, el programa espera que el usuario teclee el valor numérico o la expresión. El usuario puede teclear un vector o una matriz. En cualquier caso, la expresión introducida es evaluada cono los valores actuales de las variables de MATLAB y el resultado se devuelve como valor de retorno. >> n=input('Teclee el numero de ecuaciones') Teclee el numero de ecuaciones 4 n= 4 Otra forma: obsérvese el parámetro „s‟: >> nombre=input('¿Cómo te llamas?', 's') ¿Cómo te llamas? José nombre = José En este caso el texto tecleado como respuesta se lee y se devuelve sin evaluar, con lo que se almacena en la cadena nombre. Así, pues, en este caso, si se teclear una fórmula, se almacena como texto sin evaluarse. FUNCION DISP 30



Permite imprimir en pantalla un mensaje de texto o el valor de una matriz, pero sim imprimir su nombre. En realidad, disp. Siempre imprime vectores y/o matrices: las cadenas de caracteres son un caso particular de vectores. >> disp('El programa ha terminado') El programa ha terminado >> A=rand(4) A= 0.9501 0.2311 0.6068 0.4860

0.8913 0.7621 0.4565 0.0185

0.8214 0.4447 0.6154 0.7919

0.9218 0.7382 0.1763 0.4057

0.8913 0.7621 0.4565 0.0185

0.8214 0.4447 0.6154 0.7919

0.9218 0.7382 0.1763 0.4057

>> disp(A) 0.9501 0.2311 0.6068 0.4860

Programación con MATLAB Para escribir un programa o fichero de comandos con MATLAB habrá que crear un fichero que tenga extensión .m y contenga las instrucciones. Esto se puede hacer con cualquier editor de textos, pero tiene algunas ventajas usar el editor propio de MATLAB llamándolo desde la barra de herramientas con los menús FILE > NEW > M-FILE. Una vez estamos en el editor de funciones se puede ejecutar el código escrito con el menú DEBUG>RUN o pulsando la tecla F5, o bien escribiendo en el Command Windows el nombre de la función. MATLAB trabaja con memoria dinámica, por lo que no es necesario declarar las variables que se van a usar. Por esta misma razón, habrá que tener especial cuidado y cerciorarse de que entre las variables del espacio de trabajo no hay ninguna que se llame igual que las de nuestro programa (proveniente, por ejemplo, de un programa previamente ejecutado en la misma sesión), porque esto podría provocar conflictos. FICHEROS *.m Los ficheros con extensión (.m) son ficheros de texto sin formato (ficheros ASCII) que constituye el centro de la programación en MATLAB. Existen dos tipos de ficheros *.m, los ficheros de comandos (llamados scripts en inglés) y las funciones. Los primeros contienen simplemente un conjunto de comandos que se ejecutan sucesivamente cuando se teclea el nombre del fichero en la línea de comandos de MATLAB. 31



Las funciones permiten definir funciones enteramente análogas a las de MATLAB, con su nombre, sus argumentos y sus valores de retorno. Los ficheros *.m que definen funciones que permiten extender las posibilidades de MatLab: de hecho existen bibliotecas de ficheros *.m que se venden (toolkits) o se distribuyen gratuitamente ( a través de internet). Las funciones definidas en ficheros *.m se caracterizan por que la primera línea (que no sea comentario) comienza por la palabra función, seguida por los valores de retorno (entre corchetes [ ] y separados por comas, si hay más de uno), el signo igual ( = ) y el nombre de la función, seguido de los argumentos (entre paréntesis y separados por comas). DEFINICION DE FUNCIONES: La primera línea de un fichero llamado name.m que define una función tiene la forma: function [lista de valores de retorno] = name (lista de argumentos) donde name es el nombre de la función. Entre corchetes y separados por comas van los valores de retorno (siempre que haya más de uno) y entre paréntesis también separados por comas los argumentos. Recuerde que los argumentos son los datos de la función y los valores de retorno sus resultados. Si no hay valores de retorno se omiten los corchetes y el signo igual ( = ); si solo hay un valor de retorno no hace falta poner corchetes. Tampoco hace falta poner paréntesis si hay argumentos. Ejemplo: Se creará una función que nos evalúe el numerador, denominador y la división de: 𝑥 𝑦 + 𝑠𝑒𝑛 𝑥 ∗ 𝑦 − 𝑥 𝑃 𝑥, 𝑦 = 𝑥2 + 𝑦3 Observamos que depende de 2 entradas x e y , además tiene tres salidas: num, den y la división num/den. >> [a,b,c]=division(2,3) a= 5.7206 b= 31 c= 0.1845 Ejemplo: Crear una función de Newton („fnewton‟) que resuelva una ecuación no – lineal de una sola variable, abriendo una ventana de programación y teclear el siguiente contenido: % Función que resuelve Ecuaciones No - Lineales de una sola variable 32



% empleando el método de Newton - Raphson % _________________________________________________________________ % % [xn]=fnewton(f,xo,to1) % xn es el valor que satisface la ec. no - lineal (respuesta) % f es el nombre de la función donde se guardó la ecuación a resolver % xo es un valor inicial % to1 es el error menor a cuanto por ejemplo 0.000001 function [xn]=fnewton(f,xo,to1) x=xo; h=0.0001; error=1000; while (error>to1) y=feval(f,x); f_deriv=(feval(f,x+h)-feval(f,x))/h; xn=x-(y/f_deriv); error=abs(xn-x); x=xn; end Para correr la función („fnewton‟) necesitamos abrir uan ventana de programa y creamos una función donde guardaremos la ecuación a resolver, por ejemplo x=cos(x), y se guarda con el nombre „calcular‟: function y=calcular(x) y=x-cos(x); Y, desde el command window ejecutamos: >> x=fnewton('calcular',0,0.00001) x= 0.7391 Un programa escrito en MATLAB admite la mayoría de las estructuras de programación al uso y su sintaxis es bastante estándar. En los siguientes apartados se muestra la sintaxis de algunas de estas estructuras (if, for, while,...).

La condiciones “if” y “switch” IF

33



Una de las operaciones más habituales en programación es la ejecución de unos determinados comandos si se cumple una determinada condición lógica, es lo que se conoce como bifurcaciones. Para ello se emplea la orden “if...end” de la siguiente forma: if condición comandos end Si la bifurcación es múltiple entonces podemos combinar con los comandos else y elseif If

condicion comandos elseif condicion comandos elseif condicion comandos else comandos end Por ejemplo, crear el fichero “prueba_condicion.m” y ejecutar el siguiente código: A = 1; B = 2; if A + B == 3 disp('el resultado es 3') elseif A+B > 5 disp('el resultado es mayor que 5 ') else disp('el resultado ni es 3 ni mayor que 5 ') end (la orden “disp” se utiliza para visualizar texto en la orden de comandos) Para ejecutar la función realizada se pulsa F5 se estamos en el editor de funciones, o se escribe su nombre en la Command Window: >>prueba_condicion el resultado es 3 También se puede llamar la función recién descrita desde otra función sin más que escribir su nombre. SWITCH La sentencia switch realiza una función análoga a un conjunto de if...elseif concatenados. Su forma general es la siguiente: 34



switch expresion case case_expr1 comandos case case_expr2 comandos ... otherwise % opción por defecto comandos end Por ejemplo si queremos realizar acciones distintas en función del valor de una variable entera podemos escribir: A = 1; B = 2; switch A + B case 3 disp('el resultado es 3') case 5 disp('el resultado es 5 ') otherwise disp('el resultado ni es tres ni cinco') end Los bucles for y while FOR...END Cuando queremos ejecutar un conjunto de comandos un número predeterminado de veces utilizamos la función “for...end”, cuya forma general es: for contador = [vector de valores] comandos end A modo de ejemplo si queremos sumar los cuadrados de los números enteros impares del 1 al 11 podemos escribir: suma = 0 for k = 1:2:11 suma = suma + k^2; end suma Aunque produciría el mismo resultado la instrucción: >> suma = sum([1:2:11].^2) Ejercicio: 35



Crear un fichero llamado “primos.m” (FILE>NEW>M-FILE) para calcular la suma de los 100 primeros términos de la sucesión 1, 2x, 3x2, 4x3, ... WHILE...END A diferencia del comando for, el comando while se utiliza cuando el número de veces que se ejecuta un determinado conjunto de sentencias depende de una determinada condición lógica. Su utilización genérica es la que sigue: while condicion comandos end Por ejemplo si queremos sumar el cuadrado de números pares mientras no se supere un determinado valor, por ejemplo, 500, entonces: suma = 0; num = 2; while suma < 500 suma = suma + num^2; num = num + 2; end suma Una instrucción muy útil es el comando “break” que se utiliza para salir de un bucle. Por ejemplo, si desea salir del bucle en el caso de que alguna suma parcial sea exactamente igual a 100, entonces: suma = 0; num = 2; while suma < 500 suma = suma + num^2; num = num + 2; if suma == 100 break end end suma

Ejercicio: Dada la función y = e-x calcular el valor x que haga que se cumpla la igualdad x = e-x con un error de 10-6. Tomar como valor inicial 0.5. Poner una condición que limite a 10 el número de iteraciones para llegar a la solución (utilizar break). Además, probar el programa realizado suponiendo como condición inicial x=10. 36



Ejercicio: Se desea programar una función que calcule, mediante el algoritmo de Euclides, el máximo común divisor de dos números naturales, basta escribir un fichero euclides.m cuyo contenido sea: function m=euclides(a,b) % Cálculo del máximo común divisor de dos números naturales % mediante el algoritmo de Euclides if a0 c=rem(a,b); a=b; b=c; end m=a; Si, una vez escrito el fichero anterior, en el espacio de trabajo o en un programa se escribe la instrucción mcd=euclides(33,121) en la variable mcd se almacenará el valor 11. Las variables de una función son siempre locales. Por tanto, aunque en el seno de la función se modifiquen los argumentos de entrada, el valor de las variables correspondientes queda inalterado. Por ejemplo, en la función euclides.m se modifica el valor de los argumentos de entrada, pero, sin embargo: >>x=15; >>mcd=euclides(x,3); >>x x= 15 Si se pretende que las modificaciones de un argumento de entrada afecten a la variable correspondiente, deberá situarse dicho argumento, además, en la lista de argumentos de salida. 37



Ejercicio: Crear una función llamada “miexp” que dado un valor de x devuelva como los dos resultados de las ecuaciones: 𝑌1 = 𝑥 + 𝑒 −𝑥 𝑌2 = 𝑒 −𝑥 Utilizar dicha función para resolver la ecuación no lineal 𝑥 = 𝑥 + 𝑒 − 𝑥 mediante NewtonRaphson con un error de 10− 6 limitando el número de iteraciones a 10.

GUIA DE CLASE N° 3: GRAFICOS BI Y TRIDIMENSIONALES: Para realizar gráficos 2 – D hay cuatro funciones básicas:  PLOT( ): permite crear un gráfico a partir de vectores y/o columnas de matrices con escalas lineales sobre ambos ejes (x e y).  Loglog( ): igual que el anterior, pero con escala logarítmica.  Semilogx( ): igual que el primero con escala lineal en el eje de las ordenadas y logarítmica en el eje de las abscisas.  Semilogy( ): igual que el primero con escala lineal en el eje de las abscisas y logarítmica en el eje de las ordenadas. Se toma como referencia la primera y las demás, como consecuencia. Para agregar títulos, nombres a cada eje, cuadrículas, textos, etc., existen otras funciones: 

title (‘Gráfica P vs T’)

 añade un título al gráfico

 xlabel(‘Temp(°C)’)

 añade etiqueta al eje x (con xlabel off desaparece).

 ylabel(Presión(KPa)’)



 text (x,y,’texto’)

 introduce „texto‟ en el lugar especificado por las coordenadas x e y Si x e y son vectores, el texto se repite pro cada par de elementos. Si el texto es también un vector d cadenas de texto de la misma dimensión, cada elemento se escribe en las coordenadas correspondientes.

 gtext(„texto‟)

 introduce texto con ayuda del ratón: el cursor cambia de forma y se espera un clic para introducir el texto en esa posición.

 legend ( )

 define rótulos para las distintas líneas o ejes usados en la figura.

añade etiqueta al eje y (con ylabel off desaparece).

38


 grid on


 activa la inclusión de una cuadrícula en el dibujo. Con grid off desaparece la cuadrícula.

Ejemplo: Realizar una gráfica de un coseno multiplicado por una función exponencial Recordar dar espaciamiento entre puntos pequeños (0.01 o 0.00 son suficientes) y que los productos o divisiones se tienen que hacer elemento a elementos (.* ./ ). >> x= -5: 0.01:5; >> y= cos(x).*exp(-0.05* x); >>h=figure(1); >>plot(x,y,‟r‟); % grafica x vs y de color rojo. >>xlabel(„Eje X‟); %etiqueta del eje x >>yalbel(„Eje Y‟); % etiqueta del eje y >>title(„COSENO EXPONENCIAL‟); % título >>text(1.5, 0.75,‟cos(x)*exp(-0,05 * x)‟); % texto en la posición x,y, >>grid on; % rejilla >>legend(„Coseno Exponencial‟); %rótulo >>figure (gcf); % para que la figura pase a primer plano.

COSENO EXPONENCIAL 1.5 Coseno Exponencial 1 cos(x)*exp(-0,05 * x)

Eje Y

0.5

0

-0.5

-1

-1.5 -5

-4

-3

-2

-1

0 Eje X

1

2

3

4

5

COMANDOS hold on y hold off: Se usan en las gráficas para ciertas condiciones. El primero de ellos hace que los gráficos sucesivos respeten los ya dibujados en la figura. Es posible que haya modificación en la escala de los ejes. El comando hold off deshace el efecto de hold on 39



GUIA DE CLASE N° 4: SISTEMAS DE ECUACIONES TRIDIMENSIONALES.

Y

FUNCIONES

Interpretación Gráfica: Una ecuación lineal con dos variables, como 2x – y =3, define una recta (y=mx+b). Se puede escribir como y=2x -3. Si hay dos ecuaciones, estas pueden ser lineales y cruzarse o no. Si una ecuación tiene tres variables, x, y, z, representan un plano en un espacio tridimensional. Si hay dos ecuaciones con tres variables, pueden representar dos planos que se cruzan o no, etc.. Hiperplano: es el conjunto de puntos definido por una ecuación con más de tres variables. Y en general, un conjunto de M ecuaciones lineales con N incógnitas, donde cada ecuación define un hiperplano único que no es idéntico a ninguno otro del sistema. Si MN, el sistema está sobre especificado y no existe solución única. El conjunto de ecuaciones también se denomina sistema de ecuaciones. Un sistema de ecuaciones con una solución única es no singular, y uno que no tiene solución única se llama singular. Ejemplo: Para el sistema:

x+y=2; 2x –y = 1; x + 2y = 3;

M=3>N=2, de modo que la única solución es x =1, y=1. Si M>N el sistema puede no tener solución, tener solución única o infinidad de soluciones. Ejercicio: considere el siguiente sistema de tres ecuaciones con tres incógnitas: 3x + 2y – z =10 x + 3y + 2 z = 5 x - y - z = -1 Utilizando matrices, se puede escribir: 𝐴=

3 2 −1 −1 3 2 ; 1 −1 −1

𝑥 𝑋= 𝑦 ; 𝑧

10 𝐵= 5 −1

La forma matricial: 40



AX = B Se puede multiplicar y se obtiene el sistema original Pero, si se tiene más de tres variables, la notación se puede hacer torpe, x1, x2, x3, etc. El conjunto de ecuaciones se representa como vectores fila B y X: Ejemplo: 3 𝑥1 + 2 𝑥2

− 𝑥3 = 10

− 𝑥1 + 3 𝑥2 + 2 𝑥3 = 5 𝑥1 − 𝑥2 − 𝑥3 = −1

X=[ x1 x2 x3]

𝐴=

3 −1 1 2 3 −1 −1 2 −1

B=[ 10 5 -1]t

PROGRAMACION BASICA:

Limitaciones computacionales. Opciones de exhibición. Exactitud y precisión. Administración de archivos. Ficheros M-Files. depuración. Rutas de búsqueda y arranque en Matlab.

Errores

y

• Ejemplos prácticos (visión general) Reacción de primer orden Raíz de una función usando un m-file y fzero Usando fzero desde el Command Window Solución de sistemas de ecuaciones lineales Integración numérica y simbólica Ecuaciones diferenciales de primer y segundo orden Sistemas de ecuaciones diferenciales ordinarias 41



Reacciones múltiples CLASE N° 5: APLICACIONES DE MATLAB A LA RESOLUCIÓN DE PROBLEMAS DE INGENIERÍA TRADUCIR CON GOOGLE. DE RUMANO A ESPAÑOL Resolver: "En cálculos de Ingeniería, el calor específico de una solución de NaOH (expresado en [J /(kg * K)]), la Cp, de la solución se calcula a un valor aproximado: Cp, sol= 4190 • (1 − xNaOH) + Cp, NaOH • xNaOH donde Cp, NaOH representa el calor específico del hidróxido de sodio anhidro y xNaOH es la fracción de la masa del hidróxido en solución. El error respecto al verdadero valor de soluto concentrado local es válido para una concentración (xNaOH > 0.2). La solución inorgánica diluida en (xNaOH 6-0, 2) se toma mediante la relación: Cp, sol = 4190 • (1 − xNaOH) El Calor específico del hidróxido de sodio anhidro (expresada en [J /(kg * K)]) como trabajo, es decir temperatura, se calcula con la relación: CpNAOH = 0,1835 + 3.125 T + 0,347 T2 La función en MATLAB que puede calcular el calor específico del NaOH es: function cp = cpSuelo(x, T) % función para calcular el calor específico de una solución de NaOH %% Modo de aplicarse: % CP = CPSOL(X,T) % donde: CP – calor específico de un soluto, [J/(kg.K)] % X - fracción masica del NaOH en la solución % T - temperatura, [K] if x<=0.2 cpSol = 4190*(1 - x); else cpNaOH = 0.1835 + 3.125*T + 0.347/T.^2; cp = 4190*(1 - x) + cpNaOH*x; end % end function cpSol Para calcular el calor específico de una solución de Na OH a una fracción másica de 0,25, la temperatura de 50°C =323,15°K,el comando necesario es: >> cpSol(0.25,323.15) ans = 3.3950e+003 Densidad de una mezcla de gases: Ejercicio: Se desea escribir una función en MatLab para determinar la densidad de una mezcla de gases (CO 2, H2 y CO). Utilice una función que determine la densidad de la mezcla expresada en porcentajes molares: 25% de CO2, 60% de H2, y 15% de CO, a la temperatura de 150°C y una presión de 20 bar Solución: La densidad de una mezcla de gases reales, mez, se calculas mediante la expresión: mez = Mmez P Zmez RT Donde: Mmez representa la masa molecular de la mezcla, Zmez es el coeficiente de compresibilidad de la mezcla, T es la temperatura en [°K] y R es la constante universal de los gases igual a 8,3143 [J/(mol·K)]. La Masa molecular de la mezcla se calcula por la expresión: Mmez = 𝑛𝑗=1 𝑀𝑖 𝑋𝑖 y los coeficienes de comprensibilidad mediante la relación:

42



Zam = 𝑛𝑖=1 𝑍𝑖 𝑋𝑖 Para el caso en que la mezcla se compone de n componentes la relación de representación de una maza molecular con el coeficiente de compresibilidad Zi, y xi representa una fracción molar de cada componente en la mezcla tomada. Los datos necesarios para los calculos son presentados en una tabla:. Tabla. Datos de los componentes de la mezcla de gases Componente Masa Coeficiente Molecular de compresibilidad CO2 44,010 0,274 H2 2,016 0,303 CO 28,010 0,295 La función en MATLAB está diseñada para medir los tres argumentos siendo el primero un vector de tres componentes para determinar el contenido de CO2, H2 ¸y CO en fracciones molares, y la temperatura en °K, la presión, en bar, aplicada a la mezcla de gases. Tenemos la siguiente función en MATLAB: function ro = densMez(x,T,P) %Función para calcular la densidad de una mezcla de CO2 H2 y CO %% Modo de aplicarse: % RO = DENSMEZ(X,T,P) % donde: RO – densidad de mezcla, [kg/mc] % X - composición en fracción másica en orden CO2 H2 CO % T - temperatura, [K] % P - presión, [bar] R = 8.3143; % const.univ. de los gases, [J/(mol.K)] % masa molecular % CO2 H2 CO M = [44.01 2.016 28.01]; % coeficiente de compresibilidad % CO2 H2 CO z = [0.274 0.303 0.295]; Zmez = sum(x.*z); Mmez = sum(M./x); ro = Mmez.*P./(Zmez.*T.*R); % end function densMez Llamando a la función con nuestros datos: >> densMez([0.25 0.6 0.15],423,20) ans = 7.0688 Calcular el grosor de un cilindro Un recipiente cilíndrico es utilizado para almacenar un líquido corrosivo de una densidad media de 1.200 kg/m3. El vaso ha sido construido para trabajar con una presión máxima de 8 bar ¸y una temperatura de 60°C con un diámetro de 1.100 mm ¸y una altura de 2.100 mm, con fondo y tapa plana. Construir un programa en MATLAB para calcular el grosor del cilindro Solución: El grosor de la pared de un tanque en mm se calcula de la relación: _p = + cex + cr donde _0 es el espesor de la resistencia de la pared¸˘el elemento, cex está agregado para Llevar a Condit¸ establecido operativo y cr representan el redondeo aplicado ˘ ˘ para llevar el valor de espesor diseñar un espesor normalizado ˘ un Espesor resistente ¸ ˘ un _0 en [mm] para un elemento cilíndrico se da un tiempo relativamente ¸ ˘ tomar: _0 = pc · D 2 ' _ta − pc (10.10) unde: pc este presiunea de calcul a recipientului ˆın [MPa], D este diametrul recipientului ˆın [mm], _ta

43



este rezistent¸a admisibil˘a a materialului de construct¸ie al vasului la temperatura de calcul exprimat˘a ˆın [MPa] iar ' este coeficientul de rezistent¸˘a al ˆımbin˘arilor sudate. Adaosul pentru condit¸iile de exploatare, cex ˆın [mm] este grosimea suplimentar˘a care se adaug˘a la grosimea de rezistent¸˘a a peretelui vasului ce este susceptibil de a se subt¸ia datorit˘a coroziunii/eroziunii ˆın a¸sa fel ˆıncât s˘a asigure funct¸ionarea sigur˘a pe durata de funct¸ionare a vasului. Acest adaos este dat de relat¸ia: cex = (vce + vci + ve) · _ (10.11) unde: vce este viteza de coroziune pe fat¸a exterioar˘a a vasului ˆın [mm/an], vci este viteza de coroziune pe fat¸a interioar˘a a vasului, ˆın [mm/an], ve este viteza de eroziune ˆın [mm/an] iar _ reprezint˘a durata de funct¸ionare a aparatului exprimat˘a ˆın [ani]. Ultima component˘a a calculului grosimii de proiectare este rotunjirea cr, care t¸ine seama de valoarea grosimii standardizate imediat superioare sumei _0+cex a tablei din care se confect¸ioneaz˘a vasul ¸si de abaterea negativ˘a la grosime, cab, care se stabile¸ste din standarde funct¸ie de grosimea standardizat˘a. Calcul de rezistent¸˘a prezentat este valabil doar pentru vase cu peret¸i subt¸iri, condit¸ia care se cere a fi ˆındeplinit˘a pentru aceasta fiind: __ 0D 6 0,1 (10.12) unde __ 0 = _p − cex − cab este grosimea de rezistent¸˘a efectiv˘a. Presiunea de calcul, pc este egal˘a cu presiunea maxim˘a admisibil˘a de lucru a vasului, pm la care se adaug˘a presiunea hidrostatic˘a a coloanei de lichid din vas: pc = pm + ph (10.13) Presiunea hidrostatic˘a, ph reprezint˘a presiunea datorat˘a coloanei de lichid din recipient ¸si se calculeaz˘a cu relat¸ia: ph = 9,81 · 10−9 hl · _ (10.14) unde hl reprezint˘a ˆın˘alt¸imea coloanei de lichid iar _ densitatea lichidului. Presiunea de ˆıncercare hidraulic˘a a recipientului pph este dat˘a de relat¸ia: pph = 1,25 · pm · _20 a

_ta (10.15) unde _20 a este rezistent¸a admisibil˘a a materialului de construct¸ie al vasului la temperatura de 20°C. Presiunea de ˆıncercare hidraulic˘a calculat˘a, pphc este dat˘a de relat¸ia: pphc = pph + pha (10.16) unde pha este presiunea hidrostatic˘a a lichidului de ˆıncercare care este apa. Rezistent¸a admisibil˘a a materialului _20 a ¸si _ta se calculeaz˘a prin relat¸iile: _20 a = min __20 c

1,5 ; _20 r

2,4 _ (10.17) _ta = min _ _tc 1,5; _tr 2,4

44



_ (10.18) unde: _c ¸si _r reprezint˘a limita de curgere respectiv rezistent¸a la rupere la tract¸iune ˆın [MPa]. Rezistent¸a efectiv˘a a peretelui vasului la ˆıncercarea de presiune hidraulic˘a se calculeaz ˘a cu relat¸ia: _ef = pphc(D + s_ 0) 2 ' __ 0

(10.19) Verificarea vasului la presiune hidraulic˘a implic˘a respectarea condit¸iei: _ef 6 0,9 _20 c (10.20) Calculul de rezistent¸˘a al peretelui unui vas cilindric implic˘a parcurgerea urm˘atorului algoritm: 1. se calculeaz˘a grosimea de rezistent¸˘a _0 cu relat¸ia 10.10; 2. se calculeaz˘a adaosul pentru condit¸iile de exploatare, cex cu relat¸ia 10.11; 3. se identific˘a valoare rotunjirii cr prin alegerea grosimii standardizate imediat superioare valorii calculate _0 + cex; 4. se verific˘a dac˘a vasul este cu peret¸i subt¸iri conform condit¸iei 10.12. Dac˘a acest lucru nu se verific˘a, algoritmul este abandonat; 5. se calculeaz˘a rezistent¸˘a efectiv˘a cu relat¸ia 10.19; 6. se verific˘a condit¸ia 10.20. Dac˘a aceast˘a condit¸ie se verific˘a grosimea de proiectare a peretelui vasului, _p este corect˘a ¸si algoritmul furnizeaz˘a aceast˘a valoare. Dac˘a condit¸ia nu de respect˘a se incrementeaz˘a cu 1 mm grosimea de rezistent¸˘a efectiv˘a a vasului, __ 0 ¸si se revine la pasul 2 al algoritmului. Caracteristicile materialului de construct¸ie al vasului sunt ar˘atate ˆın tabelul 10.2. Tabelul 10.2. Caracteristici OL42 (conform SR EN 10207/95). Caracteristici Valoare _20 c , [MPa] 260 _tc , [MPa] la 100°C 225 la 150°C 216 la 200°C 206 la 250°C 186 _20 r , [MPa] 410 _tr , [MPa] la 100°C 395 la 150°C 380 la 200°C 360 la 250°C 300 grosime tabl˘a, [mm] 4, 5 (abatere la grosime: − 0,5 mm) 6, 7 (abatere la grosime: − 0,6 mm) 8, 9, 10, 12, 14, 15, 16, 18, 20, 22, 25 (abatere la grosime: −0,8 mm) 28, 30 (abatere la grosime: − 0,9 mm) 32, 35, 38, 40 (abatere la grosime: − 1,0 mm) Utilizând toate aceste informat¸ii, se scrie urm˘atorul program MATLAB: %Program de calcul al grosimii de proiectare % pentru un vas cilindric vertical

45



clear all;clc false = 0; true = 1; p_m = 0.8; % presiunea maxima admisibila, [MPa] D = 1100; % diametrul vasului, [mm] H = 2100; % inaltimea vasului, [mm] t = 60; % temperatura maxima de lucru, [grd.C] ro = 1200; % densitatea lichid stocat, [kg/mc] an = 10; % durata de exploatare, [ani] v_ex = 0.1; % viteza procesului de coroziune externa, [mm/an] v_in = 0; % viteza procesului de coroziune interna, [mm/an] v_e = 0.01; % viteza procesului de eroziune, [mm/an] fi = 0.9; % coef. suduri ro_apa = 1000; % densitate lichid de incercare, [kg/mc] % caracteristici material sig_c_20 = 260; sig_r_20 = 410; t_date = [100 150 200 250]; sig_c = [225 216 206 186]; sig_r = [395 380 360 300]; % grosime si abatere standard a tablei ds = [4 5 6 7 8 9 10 12 14 15 16 ... 18 20 22 25 28 30 32 35 38 40]; ab = -[0.5 0.5 0.6 0.6 0.8 0.8 0.8 0.8 0.8 0.8 0.8 ... 0.8 0.8 0.8 0.8 0.9 0.9 1 1 1 1]; % calcul rezistenta material sig_a_20 = min([sig_c_20/1.5 sig_r_20/2.4]); sig_c_t = interp1([20 t_date],[sig_c_20 sig_c],t); sig_r_t = interp1([20 t_date],[sig_r_20 sig_r],t); sig_a_t = min([sig_c_t/1.5 sig_r_t/2.4]); % presiunea hidrostatica, [MPa] p_h = 9.81e-8*H*ro; % presiunea de calcul, [MPa] p_c = p_m + p_h; % grosimea de rezistenta, [mm] d0 = p_c*D / (2*fi*sig_a_t - p_c); % adaos de exploatare, [mm] c_ex = (v_ex + v_in + v_e)*an; % grosimea de rezistenta calculata, [mm] d0c = d0 + c_ex; test = false; while ~test % grosimea standardizata, [mm] i = 1; while (d0c-ab(i))>ds(i) i = i + 1; end c_r = ds(i) - (d0 + c_ex - ab(i)); d_p = d0 + c_ex + c_r - ab(i); % presiunea de incercare hidraulica, [MPa] p_ph = 1.25*p_m*sig_a_20/sig_a_t; % presiunea hidrostatica apa, [MPa] p_ha = 9.81e-8*H*ro_apa; % presiunea de incercare hidraulica calculata, [MPa] p_phc = p_ph + p_ha; if (d0c/D)>0.1 error(’Vasul este cu pereti grosi!!!’); end

46



sig_ef = p_phc*(D + d0c)/(2*fi*d0c); test = sig_ef<=(0.9*sig_c_t); if ~test d0c = d0c + 1; end end; % raport disp(’ Calculul grosimii de proiectare’); disp(’ =========================================’); disp(sprintf(’ diametrul, [mm] %4d’,D)); disp(sprintf(’ inaltimea, [mm] %4d’,H)); disp(sprintf(’ presiunea, [MPa] %5.1f’,p_m)); disp(sprintf(’ temperatura, [grd.C] %3d’,t)); disp(’ -----------------------------------------’); disp(sprintf(’ rezistenta mat. 20C, [MPa] %5.1f’,sig_a_20)); disp(sprintf(’ rezistenta mat.%3dC, [MPa] %5.1f’,t,sig_a_t)); disp(’ -----------------------------------------’); disp(sprintf(’ presiunea de calcul, [MPa] %7.3f’,p_c)); disp(sprintf(’ pres.de incercare hidr., [MPa] %7.3f’,p_phc)); disp(’ -----------------------------------------’); disp(sprintf(’ grosimea de proiectare, [mm] %2d’,d_p)); disp(’ =========================================’); % end script Dup˘a rularea programului obt¸inem urm˘atorul raport: Calculul grosimii de proiectare ========================================= diametrul, [mm] 1100 inaltimea, [mm] 2100 presiunea, [MPa] 0.8 temperatura, [grd.C] 60 ----------------------------------------rezistenta mat. 20C, [MPa] 170.8 rezistenta mat. 60C, [MPa] 161.7 ----------------------------------------presiunea de calcul, [MPa] 1.047 pres.de incercare hidr., [MPa] 1.263 ----------------------------------------grosimea de proiectare, [mm] 6 ========================================= 10.2.2. Calculul puterii necesare la antrenarea unui amestec˘ator Exercit¸iul 10.5. Se cunosc urm˘atoarele informat¸ii despre un reactor prev˘azut cu un sistem de amestecare cu elice: • presiunea maxim˘a de lucru din vas: 0,8 MPa; • densitatea mediului de lucru: 1.200 kg/m3; • vâscozitatea dinamic˘a a mediului de lucru: 5,5 cP; • diametrul elicei: 0,4 m; • turat¸ia: 300 min −1; • diametrul arborelui: 60 mm; • l˘at¸imea unui inel de etan¸sare: 20 mm. S˘a se scrie un program MATLAB pentru calculul puterii electromotorului de act¸ionare al sistemului de amestecare. Rezolvare: Puterea efectiv˘a de act¸ionare a unui amestec˘ator, N, se stabile¸ste funct¸ie de mai multe elemente: puterea necesar˘a amestec˘arii, Nam, care exprim˘a puterea efectiv ˘a transmis˘a mediului de lucru, puterea pierdut˘a prin frecare ˆın cutia de etan¸sare, Nf ¸si puterea pierdut˘a ˆın vas datorit˘a amenaj˘arilor interioare (teac˘a de m˘asurare a temperaturii, serpentin˘a, sp˘arg˘ator de vârtej, rugozitatea peret¸ilor, ¸s.a.m.d.). Relat¸ia de calcul utilizat˘a este [20]:

47



N=k _t · h Nam _ 1+ X 'i _ + Nf i (10.21) unde: k reprezint˘a coeficientul de majorare a puterii care t¸ine seama de suprasarcinile care apar la pornire, umplere sau alimentare, _t este randamentul total al sistemului de act¸ionare iar 'i reprezint˘a factori de majorare datorate prezent¸ei unor amenaj˘ari interioare. În cazul P'i = 0,4. nostru, valorile acestor coeficien¸ti sunt: k = 1,25, _t = 0,95 ¸si Puterea necesar˘a amestec˘arii, ˆın [kW], pentru alte amestec˘atoare decât cele de tipul impeller, se calculeaz˘a cu relat¸ia: Nam = np _n 60 _3 d5 · _ · 10−3 (10.22) unde: np este num˘arul caracteristic de putere, n reprezint˘a turat¸ia amestec˘atorului [min−1], d este diametrul amestec˘atorului, [m], iar _ este densitatea lichidului antrenat. Num˘arul caracteristic de putere, np, se stabile¸ste funct¸ie de criteriul Reynolds conform datelor prezentate ˆın tabelul 10.3. Valoarea criteriului Re se calculeaz˘a cu expresia: Re = 16,25 · n · d2 · _ _ (10.23) unde _ reprezint˘a vâscozitatea mediului de lucru, ˆın [cP]. Tabelul 10.3. Num˘arul caracteristic de putere pentru amestec˘atoare cu elice. Valoarea Num˘arul caracteristic criteriului Re de putere, np 102 1 103 0,5 104 0,3 105 0,27 106 0,22 Puterea consumat˘a prin frecare ˆın cutia de etan¸sare a sistemului de amestecare, Nf ˆın [kW], se calculeaz˘a cu relat¸ia: Nf = 2,255 · 10−7 · d2 a · m · b · μ · pm · n (10.24) unde: da este diametrul arborelui amestec˘atorului ˆın dreptul cutiei de etan¸sare, ˆın [mm], m reprezint˘a num˘arul de inele de etan¸sare, b esteˆın˘alt¸imea unui inel de etan¸sare, ˆın [mm], μ este coeficientul de frecare ¸si pm reprezint˘a presiunea maxim˘a admisibil˘a de lucru ˆın vas, [MPa]. Valorile utilizate ˆın acest caz sunt: m = 4, b = 20 [mm] ¸si μ = 0,075. Puterea calculat˘a a electromotorului, conform relat¸iei 10.21, este folosit˘a pentru alegerea puterii electromotorului dintr-o list˘a de puteri disponibile. Puterile electromotoarelor disponibile pe piat¸˘a formeaz˘a urm˘atorul ¸sir de valori, ˆın [kW]: 0,8; 1,1; 1,5; 2,2; 3; 4; 5,5; 7,5; 10; 13; 17; 22; 30; 40; 55; 70. Utilizând aceste informat¸ii, se construie¸ste urm˘atorul program MATLAB: %Program de calcul al puterii electromotorului

48



% de antrenare al unui amestecator cu elice clear all;clc false = 0; true = 1; p_m = 0.8; % presiunea maxima admisibila, [MPa] d = 0.4; % diametrul elicei, [mm] d_a = 60; % diametrul arbore, [mm] n = 300; % turatia elicei, [1/min] ro = 1200; % densitatea mediului de lucru, [kg/mc] eta = 5.5; % vascozitatea mediului de lucru, [cP] k = 1.25; % coef. majorare putere rt = 0.95; % randament totalal sistemului de actionare sfi = 0.4; % suma amenajari interioare m = 4; % numar inele de etansare b = 20; % inaltimea inelului de etansare, [mm] c_f = 0.075; % coeficient de frecare Re_t = [1e2 1e3 1e4 1e5 1e6]; np_t = [1 0.5 0.3 0.27 0.22]; P_disp = [0.8 1.1 1.5 2.2 3 4 5.5 7.5 ... 10 13 17 22 30 40 55 70]; % calculul coeficientului de putere Re = 16.25*n*d^2*ro/eta; np = interp1(log(Re_t), np_t, log(Re)); % calculul puterii necesare pentru amestecare N_am = np * (n/60)^3 * d^5 * ro * 1e-3; % calculul puterii consumate prin frecare N_f = 2.255e-7 * d_a^2 * m * b * c_f * p_m * n; % calcul puterii efective de actionare N = k/rt*(N_am*(1 + sfi) + N_f); disp(sprintf(’Puterea efectiva %6.2f kW’,N)); % alegerea puterii electromotorului i = 1; while N>P_disp(i) i = i + 1; end; P = P_disp(i); disp(sprintf(’Putere electromotor %6.2f kW’,P)); % end script Apelând acest program, obt¸inem urm˘atorul rezultat: Puterea efectiva 2.27 kW Putere electromotor 3.00 Kw

FUNCION plot3( ) La función plot3 es análoga a su homóloga bidimensional plot. Su forma más sencilla es la siguiente: >> plot3(x,y,z) que dibuja una línea que une los puntos (x(1), y(1), z(1)), (x(2), y(2), z(2)), etc. y la proyecta sobre un plano para poderla representar en la pantalla. Ejemplo: 49



x=cost(t); y=cost(t); z=t;

DIBUJO DE MALLADOS: FUNCIONES meshgrid(), mesh(), surf() Meshgrid

50



51



52



DIBUJO DE CURVAS DE NIVEL: FUNCIONES contour, contour3, surf.

Siguiendo con el ejemplo anterior:

53



54



Otros tipos de ploteos tridimensionales de Funciones: La table de abajo mestra todas la variedad de graficaciones de funciones y volúmenes 3D. Inclusive funciones que generan data 3-D (cilindro, elipsoide, esfera), ademas otros arreglos de data o funciones. Para informaciónacerca de los ejemplos del uso de funciones gráficas 3-D, ver “Creación de Grafics 3 – D en la documentación Visualización 3 – D.

Líneas

mallas

GRAFICOS DE: áreas superficies

direcciones

volumétricos

55



Sistemas de ecuaciones no lineales. Comando fsolve()

mifuncion.m

56



En la ventana de comandos:

Integración quad(@f, a, b)Integra por cuadraturas la función f, desde ahasta b.

Diferenciación numérica: diff(x) Si x = [x(1) x(2) …x(n)] entonces diff(x) = [x(2)-x(1) x(3)-x(2) …x(n)-x(n-1)] EDOs de primer orden

57



Límites

58



Derivación e integración

59



60



61



62



Series de Taylor

63



Ecuaciones diferenciales ordinarias

64



Simulación de un reactor no isotérmico

65



66



67



68



69



70



71

CURSO PARA INGENIEROS.pdf

Recommend Documents