Programación en Ingeniería Química Con OCTAVE. Guillermo Peris Ripollés, Universitat Jaume I.

Tema 1

Introducci´ on on a Octave Guil Guil lermo lerm o Peris Peri s Ripol Ripoll´ l´ es es Objetivos Cuando finalice este tema, el alumno deberá ser capaz de: Realizar Realizar c´ alculos alculos simples con Octave, utilizando los operadores aritméticos eticos báa-sicos, sus reglas de precedencia y algunas funciones matemáticas. Asignar valores a variables y utilizarlas correctamente. Gestionar el entorno de trabajo Octave para cono c onocer cer las la s variables variable s definidas, defi nidas, as´ as´ı como gestionar ficheros y directorios. Crear ficheros-M con programas que pidan información on al usuario y proporcionan informaci´ on on por la salida estándar. andar. Utilizar la ayuda de Octave para obtener información on sobre su uso.

Aplicaci´ on on Cuando finalice este tema, el alumno deberá ser capaz de resolver problemas como el siguiente, cuya resolución on se indica a lo largo del propio tema. C´ alculo alculo de pesos moleculares Escribe un programa en Octave que pida al usuario el número umero de átomos atomos de carbono, ox´ ox´ıgeno e hidrógeno ogeno de una molécula ecula org´ anica anica (sin otro tipo de átomos) atomos) y calcule calcule su peso molecula molecularr y el porcentaje porcentaje en peso de ox´ ox´ıgeno. ıgeno. Como ejemplo, aplica el programa a las siguientes mol´ eculas: eculas: a) Aspirina Aspirina (C (C9 H6 O4 ). b) Benceno Benceno (C (C6 H6 ). c) Alcohol Et´ Et´ılico (CH3 CH2 OH). ´ o Acético d) Acido Acid eti co (CH 3 COOH). e) Acetona Acetona (CH (CH3 COCH 3 ). Pesos atómicos: omicos: Carbono, 12 g/mol; Hidr´ ogeno, ogeno, 1 g/mol; y Ox´ Ox´ıgeno, 16 g/mol.

Ingenier´ıa Qu´ımica

Programación en Octave

1-2

Introducción on a Octave

Contenidos

1.1. Introduc Introducci´ ci´ on on . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-3 1.2. Conceptos b´ asicos asicos . . . . . . . . . . . . . . . . . . . . . . . . 1-3 1.2. 1.2.1. 1. C´ alculos alculos simples . . . . . . . . . . . . . . 1.2.2. 1.2.2. Variables ariables . . . . . . . . . . . . . . . . . . 1.2.3. 1.2.3. Comentario Comentarioss y signos signos de puntuaci puntuación ón . . . 1.2.4. Formato num´ erico erico . . . . . . . . . . . . . 1.3. Funciones unciones matem´ matem´ aticas aticas . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . . . . . . . . . . . . . . . .

11 -3 11-5 11-6 11 -7 1-8

1.3.1. 1.3.1. Funciones unciones matem´ matem´ aticas aticas elementales . . . . . . 1.3.2. Funciones trigonométricas etricas . . . . . . . . . . . 1.4. 1.4. El entor entorno no de trab trabajo ajo Octave . . . . . . . . . 1.4.1. 1.4 .1. Gesti´ Gesti´ on on de variables . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . 11 -9 . . 1-9 . 1-10 . . 1-10

1.4.2. Histor 1.4.2. Historia ia de la la sesi´ sesi´ on on . . . . . . . . . . . . . . . 1.5. Fichero Ficheros-M s-M . . . . . . . . . . . . . . . . . . . . . 1.6. Instrucci Instrucciones ones de entrada entrada/salid /salida a . . . . . . . . . 1.6.1. 1.6.1. Salidas Salidas por pantalla pantalla . . . . . . . . . . . . . . . 1.6.2. 1.6 .2. Entra Entrada da de dato datoss . . . . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . 1-11 . 1-11 . 1-12 . . 1-12 . . 1-14

1.7. 1.7. Ayud Ayuda a de de Octave . . . . . . . . . . . . . . . . . . . . . . . . 1-15 1.8. Ejercicios Ejercici os pr´ acticos acticos . . . . . . . . . . . . . . . . . . . . . . . 1-17

Universitat Jaume I

Guil lermo Peris Ripol l´ es

1.1 Introducción

1.1. 1.1.

1-3

Intr Introdu oducc cci´ i´ on on

Octave es

un potente programa de cálculo alculo matem´ atico basado en software libre, origiatico nalmente desarrollado para traba jar con vectores y matrices por ingenieros qu´ qu´ımicos, aunqu aunquee sus sus ultimas u ´ ltimas versione versioness exceden exceden con mucho mucho este propósito osito inicial. inicial. Incorpora Incorpora distintas funcionalidades utiles u ´ tiles para el cálculo alculo técnico ecnico y cient´ cient´ıfico, como algebra algeb ra matricial, tricial, manipulaci´ manipulaci´ on on de polinomios, cálculo alc ulo numérico, eric o, cálculo alculo simbólico, olico, y creación on 1 y manipulación on de gráficos aficos . Además, as, incluye un lenguaje de programación on propio, que a´ un no siendo tan eficiente como los lenguajes compilados (como C/C++, Forun tran, etc), es más as fácil acil de utilizar y nos permitirá introducir en este curso conceptos b´ asicos asicos de programación on comunes con estos lenguajes. Este lenguaje es muy similar al que proporciona el software propietario Matlab, por lo que puede consultarse la bibliograf´ biblio graf´ıa ıa de este ultimo u ´ ltimo programa (mucho m´ as extensa) para aprender a utilizar as Octave. En este curso se asume que se va a ejecutar Octave en un sistema operativo basado en GNU/Linux.

1.2.

Conceptos b´ asicos asicos

Para ejecutar Octave bajo Linux abriremos un terminal bajo Linux y ejecutaremos en él el la orden orde n octave. Al hacerlo, aparecerá una ventana en la que, tras información on 2 sobre el programa, aparecerá el indicador de órdenes ordenes octave:1> . En esta ventana unicamente u ńicamente podremos escribir despu´ es es del último ultimo indicador aparecido.

1.2.1.

C´ alculos alculos simples

Para empezar a trabajar con Octave, vamos a resolver poco a poco la aplicación on práctica actica expuesta al principio de la unidad. Para ello, calcularemos inicialmente el peso molecular y el porcentaje p orcentaje en peso de ox´ ox´ıgeno del ácido acido acetilsalic acetils alic´´ılico (aspirina), (aspirin a), de fórmula ormula C9 H6 O4 , siendo los pesos atómicos omicos aproximados: C, 12 g/mol; H, 1 g/mol; y O, 16 g/mol. Para resolver este problema debemos calcular primero el peso molecular del compuesto. puesto. Podemos Podemos ayudarnos ayudarnos con Octave utilizándolo andolo como una calculadora3 :

 9*12 + 6*1 + 4* 16

ans = 178

Fijémonos emo nos en que Octave almacena el resultado con el nombre ans, abreviatura de answer (soluci´ (solución). on). Si ahora calculamos los gramos de ox´ ox´ıgeno por p or mol de producto,

 4* 16

ans = 64

podemos obtener el porcentaje en peso p eso de ox´ ox´ıgeno como

 64 / 178 * 100

ans = 35.955 35.955 1 2 3

Realmente, para la representaci´ on on gr´ afica afica Octave se apoya en otro programa libre: Gnuplot. Para simplificar simplificar esta notaci´ on, en los apuntes utilizaremos el prompt ( ). on, F´ıjate en que los espacios no son importantes en Octave.



1-4

Introducción a Octave En el ejemplo anterior hemos utilizado los operadores aritméticos suma, multiplicación y división. Los operadores aritméticos básicos que proporciona Octave son los siguientes4: Operador suma, a + b resta, a - b multiplicaci´ on, a * b divisi´ on, a b exponenciación, ab cambio de signo, -a

S´ımbolo + * /o

÷

∧

-

\

Ejemplo 3 + 22 90 - 54 3.14*0.85 56/8 2∧ 8 -(4*8)

Tabla 1.1: Operadores aritméticos básicos de Octave.

La ejecución de estas operaciones en una expresión matem´ atica se realiza siguiendo las reglas de precedencia , y de izquierda a derecha en el caso de operadores con igual precedencia. La mayor precedencia la posee el operador de exponenciación, seguido por los operadores de multiplicación y división (igual precedencia), siendo la suma y la resta los operadores con menor precedencia (igual entre ellos). Por ejemplo, la expresión 2 + 3 4∧ 2

∗

se evalua considerando primero la exponenciación (4∧ 2 = 16), por ser el operador de mayor precedencia, a continuación la multiplicació n (3 16 = 48), y por u ´ ltimo la suma (48 + 2 = 50). Este orden de precedencia puede ser alterado mediante el uso de paréntesis. As´ı, (2 + 3) 4∧ 2 = 80. En caso de duda sobre el orden de precedencia, es recomendable hacer uso de paréntesis.

∗

∗

Ejercicios 1

b

Eval´ ua las siguientes expresiones según las reglas de precedencia de Octave.



a) 2 + 3 + 4 + 2 b) 2 + 3 * 4 + 2 c) (2 + 3) * 4 + 2 d) (2 + 3) * (4 + 2) e) 1/2*( 2 +3 *4∧ 2) b

2 Traduce las siguientes expresiones matemáticas a Octave utilizando el menor número de paréntesis posible. 

a) 2 + (3 (6/2)) b)

4+6 2+3

·

c) (4/2)5 d) (4/2)5+1 e) ( 3)2

−

4

El operador \ realiza operaciones de matrices a izquierdas, aunque no lo utilizaremos en este curso.

Universitat Jaume I

Guillermo Peris Ripoll´ es

1.2 Conceptos básicos

1.2.2.

1-5

Variables

Una variable es un identificador que se utiliza para representar cierto tipo de informaci´ on dentro de una determinada parte del programa. En su forma más sencilla, una variable es un identificador que se utiliza para representar un dato individual, como por ejemplo una cantidad numérica. Este valor se puede recuperar después en el programa simplemente haciendo referencia al nombre de la variable, e incluso cambiar su contenido por otro del mismo tipo. Podemos resolver el problema anterior utilizando variables y asignándoles valores con el operador de asignación (=). De esta forma, podemos reutilizar los resultados intermedios para operaciones posteriores. En las siguientes l´ıneas se resuelve el problema definiendo 5 variables (peso_C , peso_O, peso_H, peso_molecular y porcentaje_O ).

 peso_C = 9*12;  peso_O = 4*16;  peso_H = 6*1 ;  peso_molecular = peso_C + peso_O + peso_H peso_molecular = 178  porcentaje_O = peso_O/peso_molecular*100 porcentaje_O = 35.955

En las tres primeras l´ıneas de código se observa que la utilizació n de un punto y coma (;) al final de la expresión matem´ atica evita que se muestre el resultado en pantalla. Esta opci´ o n de Octave se utiliza con frecuencia para evitar un exceso de l´ıneas de salida de programa. En las órdenes de asignaci´ on, en primer lugar se calcula la expresión a la derecha del signo =, y el valor obtenido se asigna a la variable que aparece a la izquierda del operador. As´ı, mientras la expresión

 peso_C

= 9*12;

es válida en Octave, la siguiente no lo es:

 9*12 =

peso_C ;

Adem´ as, teniendo en cuenta el orden de evaluación de las asignaciones, la siguiente secuencia de órdenes es correcta:

 a = 10;  a = 2*a a = 20

ya que en la segunda l´ınea se calcula en primer lugar 2*a y el resultado (20) se asigna a la variable a, sustituyendo el valor anterior. La utilización de nombres de variables en Octave sigue una serie de reglas que se indican a continuaci´ on: Variables que sólo se diferencian en el uso de mayúsculas/min´ usculas son distintas. Ejemplo: peso o ,peso O y Peso O se pueden utilizar como tres variables distintas. Ingenier´ıa Qu´ımica


1-6

Introducción a Octave Los nombres deben empezar por letras, y pueden contener tanto letras como n´ umeros y guiones bajos ( ). Los signos de puntuación no están permitidos. Ejemplo: Este es 1 nombre valido y 0Este,no lo es. Resulta conveniente escoger nombres que nos indiquen claramente para qu´ e se utiliza la variable. Además, hay que tener en cuenta que Octave utiliza una serie de variables predefinidas cuyo valor no conviene cambiar, aunque es posible hacerlo. Algunas de estas variables son: ans pi eps inf NaN o nan i y j e

Variable que almacena el u ´ ltimo resultado. Valor π = 3.1415.... Menor número que, sumado a 1, da un resultado mayor que 1. Infinito, p.e., 1/0. Resultado no numérico ( Not a number ), p.e., 0/0. 1. Base de los logaritmos naturales (neperianos).

√ −

Tabla 1.2: Variables predefinidas por Octave.

Ejercicios b

3 Indica cu´ ales de las siguientes expresiones son incorrectas y por qué.



a) numero bajas = 6 + 2 ; b) 8colores = 6*8; c) numero bajas = 6 + 2; d) i = 1 ; e) A1234 5678 = i ; f) A1234 5678 = A1234 5678*2 g) B-52 = 0 ; b

4

¿Qué resultado esperar´ıas encontrar para el valor de x e y tras la ejecución del código Octave siguiente? 

 x = 5;  x = 2*x;  y = x 2;  x = y/x; ∧

1.2.3.

Comentarios y signos de puntuaci´ on

Octave permite

introducir varias órdenes en una misma l´ınea separándolas por signos de puntuación. As´ı, haciendo uso de comas podemos escribir,

 peso_C =

9*12, peso_O = 4*16, peso_H = 6*1 peso_C = 108 peso_O = 64 peso_H = 6

Universitat Jaume I


1.2 Conceptos básicos

1-7

Ya hemos visto que la utilización on de un punto y coma evita que se muestre el resultado en pantalla. pantalla. Podemos introducir introducir varias órdenes ordenes en una un a misma mi sma l´ınea, terminando termina ndo con coma o fin de l´ınea aquellas cuyo resultado resultado queremos visualizar, visualizar, y con punto punto y coma aquellas que no nos interesa observar.

 peso peso_C _C =

9*12 9*12, , peso peso_O _O = 4*16 4*16; ; peso peso_H _H = 6*1 6*1 peso_C peso_C = 108 peso peso_H _H = 6

En Octave, al igual que en otros lenguajes de programación, se pueden añadir nadir comentarios comentarios a los programas, lo cual permite que otros programadores programadores entiendan entiendan más as fácilmente acilmente la estructura del código, odigo, y que incluso los propios autores de los programas comprendan sus propios trabajos al cabo del tiempo. Para ello, en Octave se utiliza el carácter acter % para iniciar comentarios comentarios,, de forma que el intérpre er prete te de Octave ignora todo lo que sigue a este signo hasta el final de la l´ınea.

 % Calcu Calculo lo de los g/mol g/mol de cada cada eleme elemento nto. . peso_C C = 9*12; 9*12; % Carbon Carbono o  peso_  peso_ peso_O O = 4*16; 4*16; % Oxigen Oxigeno o  peso peso_H _H = 6*1 6*1 ; % Hidr Hidrog ogen eno o Calculo lo del peso peso molecu molecular lar  % Calcu  peso_ peso_mol molecu ecular lar = peso_C peso_C + peso_ peso_O O + peso_ peso_H H peso_mol peso_molecul ecular ar = 178 Calculo lo del porce porcenta ntaje je de oxige oxigeno no  % Calcu  porcentaje_O porcentaje_O = peso_O/peso_molecu peso_O/peso_molecular*100 lar*100 porcenta porcentaje_O je_O = 35.9551 35.9551

Ejercicios 

5 Ejecuta el programa programa anterior anterior con Octave, comprobando que el resultado es correcto.

Por ultimo, u ´ ltimo, podemos utilizar puntos suspensivos (...) para continuar una orden larga en otra l´ınea. ınea. En este caso, no se permite la división on de nombres y continuación on de comentarios: peso_mol molecu ecular lar =  peso_

peso_C peso_C + ... ... % Conti Continua nua en sigui siguient ente e linea linea

peso_O peso_O + peso_H peso_H peso_mol peso_molecul ecular ar = 178

1.2.4.

Formato Formato num´ num´ erico erico

Al hablar del lenguaje de programación on Octave, al contrario que en otros lenguajes, no tiene demasiado sentido hablar de tipos de variables (entero, real, doble precisión, etc), ya que Octave almacena todas las variables con doble precisión (con 16 decimales), aún un cuando sean enteros simples5 . Sin embargo, s´ s´ı que podemos cambiar el formato con el que se muestran las variables. Para ello, podemos utilizar las órdenes que se muestran en la Tabla 1.3 1.3,, en la que se muestra la representación de pi tras ejecutar las ordenes órdenes respectivas. 5

Esto sorprender´ sorprender´ a a aquellos que conozcan otros lenguajes de programación, como C, en los que hay que declarar el tipo de variable. En Octave, todas las variables se presuponen que son de doble precisi´ on, por lo que no es necesaria la declaración de tipos. on,



1-8

Introducción on a Octave Orden format format format forma format t

pi short long short e long long e

3.1416 3.14159265358979 3.1416e+00 3.141 3.141592 592653 653589 58979e 79e+00 +00

Tabla 1.3: 1.3: Distintos formatos numéricos ericos de Octave.

Cabe insistir en que estas órdenes unicamente u ´ nicamente afectan a la impresión on de un número umero en pantalla, y no a la representación on interna, que siempre es de 16 d´ıgitos. Ejercicios 6 Comprueba Comprueba el efecto efecto de los formatos formatos de la Tabla 1.3 sobre la variable pi. Para ello, ejecuta ejecuta primero primero la orden format correspondiente antes de mostrar el valor de pi . 

1.3.

Funciones matem´ aticas aticas

Hasta ahora hemos estudiado las operaciones aritméticas eticas b´ asicas (suma, resta, multiasicas plicación, on, división on y exponenciación), on), que son las más as utilizadas en cálculos alculos básicos. asicos. Sin embargo, en ocasiones ocasiones son necesarias necesarias operaciones operaciones más as complejas, como cálculo alculo de raices cuadradas, operaciones o peraciones trigonométricas, etricas, valores absolutos, etc. Para ello, Octave dispone de una serie de funciones predefinidas que nos permiten realizar estas operacione ope raciones. s. As A s´ı, para calcular calcul ar la ra´ ra´ız cuadrada cuadrad a de d e un u n n´ n umero u ´ mero utilizaremos la función on sqrt, introduciendo el n´ umero umero entre paréntesis entesis después es del nombre de la función:

 a = 18;  b = sqrt sqrt(a (a) ) b = 4.24 4.2426 26  sqrt(4)

ans = 2

A los parámetros ametros (constantes, variables o expresiones) que se le pasan a una función entre par´ entesis entesis se les denomina argumentos . Una función o n puede no tener ningún un argumento, uno (sqrt), o varios argumentos separados por comas. Por ejemplo, la función on rem calcula el resto de la división on entera de dos números umeros enteros:

 rem(25,4)

ans = 1

ya que si dividimos 25/4, el resultado es 6 con un resto de 1. Por supuesto, los argumentos de una función on tambi´ en en pueden contener otras llamadas a otras funciones,

 rem(sqrt(16),4)

ans = 0

Universitat Jaume I

Guil lermo Peris Ripol l´ es

1.3 Funciones matemáticas

1.3.1.

1-9

Funciones matem´ aticas aticas elementales

En la Tabla 1.4 1.4 se se proporciona una lista de algunas funciones matemáticas aticas simples, simples, como raices raices cuadradas, cuadradas, exponenciales exponenciales y logaritmos, logaritmos, as´ as´ı como otras utilizadas utilizadas para redondear redondear n´ n umeros. u ´ meros. Función sqrt(x) exp(x) log(x) log10(x) rem(x,y) abs(x) round(x) floor(x) ceil(x)

Ejemplo Orden Resultado Ra´ız ız cuadra cua drada da de x de x 5 sqrt(25) x 2.7183 e (e = base de logaritmos neperianos) exp(1) ln( ln(x), logaritmo neperiano de x de x log(2.7183) 1.0000 log10 (x), logaritmo de x de x en base 10 log10(100) 2 Resto de la divisi´ división x/y on x/y 3 rem(17,7) Valor absoluto de x de x 3.4000 abs(-3.4) Redondea x Redondea x al entero más as cercano cercano round(4.7) 5 Redondea Redondea al entero entero inferior inferior 4 floor(4.7) Redondea Redondea al entero entero superior superior 5 ceil(4.3) Comentario

Tabla 1.4: 1.4: Funciones matemáticas aticas elementales.

1.3.2.

Funciones trigonom´ etricas etricas

En Octave, se asume que los lo s argumentos de funciones f unciones trigonométricas etricas están a n en ra6 dianes . En la Tabla 1.5 Tabla 1.5 se ejemplifica el uso de algunas funciones trigonométricas. etricas. Funci´ on sin(x) cos(x) tan(x) asin(x) acos(x) atan(x)

Comentario

Orden

Ejemplo

Resultado Seno de x sin(30*pi/180) 0.50000 Coseno Coseno de x de x cos(30*pi/180) 0.86603 Tangente de x de x tan(30*pi/180) 0.57735 Arcoseno de x de x asin(0.5)*180/pi 30.000 Arcocoseno Arcocoseno de x de x acos(0.87)*180/pi 29.541 Arcotangente de x de x atan(0.58)*180/pi 30.114 Tabla 1.5: 1.5: Funciones Funcio nes trigo t rigonom´ nométricas. etri cas.

F´ıjate en que el resultado de las tres últimas ultimas funciones se encuentra en radianes, por lo que es necesario transformarlo para obtener grados sexagesimales. Ejercicios b

7 Eval´ ua las siguientes expresiones sin utilizar Octave, y después ua es comprueba c omprueba tu respuesta re spuesta haciendo uso de él. el. Si has realizado el ejercicio anterior, ejecuta antes la orden format. 

a) round(-2.6 round(-2.6))

e) round(2.6) round(2.6)

b) floor(-2.6) floor(-2.6)

f) floor(2. floor(2.6) 6)

c) abs(-2.6) abs(-2.6)

g) abs(2.6) abs(2.6)

d) ceil(-2.6) ceil(-2.6)

h) ceil(2.6) ceil(2.6)

6

Recordemos Recordemos las expresiones expresiones para transformar transformar radianes a grados sexagesimales sexagesimales y viceversa: viceversa: radianes = radianes = grados grados ∗ ∗ pi/ pi/180 180 grados = grados = radianes radianes ∗ ∗ 180 180/pi /pi . .



1-10

Introducción a Octave

b

i) rem(5,4)

m) round(exp(2*log(3)))

j) rem(17,7)

n) ceil(4*sin(3*pi/2))

k) log10(100) + log10(0.001)

n ˜ ) abs(log10(0.01)*sqrt(25))

l) sqrt(rem(89,10))

o) sin(pi)∧ 2 + cos(pi)∧ 2 -1

8



El flujo de gas que escapa de un tanque a presi´ on P 0 y en condiciones adiabáticas es: ψ =

      − k

k

−1

P e P 0

2/k

P e P 0

(k +1)/k

done P e es la presión externa y k la constante del gas reversible adiabática. Escribe esta ecuaci´ on en notación Octave y comprueba si es correcta con los valores k = 1.4 y P e /P 0 = 0.3. (Respuesta: ψ = 0.4271.)

1.4.

El entorno de trabajo Octave

Se conoce como entorno de trabajo Octave a la memoria del programa donde residen tanto las variables utilizadas en una sesión de trabajo como el historial de las órdenes utilizadas hasta el momento. El entorno de trabajo es, por lo tanto, el que nos permite recuperar valores de variables y recordar órdenes ya ejecutadas. A continuación, vamos a estudiar algunas de las caracter´ısticas y utilidades de dicho entorno.

1.4.1.

Gesti´ on de variables

Podemos conocer el conjunto de variables definidos durante una sesión de trabajo con la orden who,

 who

*** local user variables: __nargin__ peso_C peso_molecular porcentaje_O

peso_O

peso_H

La orden whos da información más detallada acerca de las variables, como la clase y el tama˜ no que ocupan en memoria.

 whos

*** local user variables: Prot Name Size ==== ==== ==== rw- __nargin__ 1x1 rwd peso_C 1x1 rwd peso_H 1x1 rwd peso_O 1x1 rwd peso_molecular 1x1 rwd porcentaje_O 1x1

Bytes ===== 8 8 8 8 8 8

Class ===== scalar scalar scalar scalar scalar scalar

Total is 6 elements using 48 bytes

Universitat Jaume I


1.5 Ficheros-M

1-11

En este ejemplo, todas las variables definidas son escalares de dimensión 1x1. La orden clear borra variables del espacio de trabajo Octave. En el siguiente ejemplo, se elimina la variable peso_C, de forma que un intento de utilización posterior de dicha variable genera un mensaje de error.

 clear  who

peso_C

Your variables are: __nargin__ peso_O peso_H peso_molecular porcentaje_O peso_C ??? Undefined function or variable ’peso_C’.



Pero hay que ir con cuidado con esta orden, ya que si se ejecuta sin ningún nombre de variable borra todos las variables y funciones de la memoria.

1.4.2.

Historia de la sesi´ on

en Octave tambi´

memoriza las expresiones introducidas en una sesión de trabajo en el mismo orden en que fueron ejecutadas, de forma que podemos navegar por la historia de la sesi´ on repitiendo órdenes ya ejecutadas. Para ello se hace uso de los cursores: con las teclas y recordamos las expresiones introducidas hacia el principio o el final, respectivamente, y con las teclas y podemos editar la orden para modificarla convenientemente.

↑ ↓

1.5.

← →

Ficheros-M

En el caso de problemas sencillos, como los tratados hasta ahora, es útil introducir directamente las órdenes en la ventana de Octave, pero en el caso de problemas más complejos o en la repetición de las mismas órdenes con distintos valores de variables, esta técnica resulta muy incómoda. Para automatizar la ejecución de secuencias de órdenes, Octave permite utilizar ficheros-M , cuyo nombre proviene de su extensión “.m”. La secuencia de órdenes contenida en un fichero-M constituye un programa , que podremos ejecutar fácilmente cuando lo deseemos. Para crear un fichero-M utilizaremos un editor de texto cualquiera guardándolo como nombre.m, siendo nombre el nombre del programa. En el ejemplo siguiente se ha guardado el fichero-M con el nombre pmol.m. %***************************************************************** % Programa : pmol.m % Descripcion: Este programa calcula el peso molecular de una % molecula organica. %***************************************************************** % Calculo de los g/mol de cada elemento. peso_C = 12*9 ; peso_H = 1*6 ; peso_O = 16*4; % Calculo del peso molecular peso_molecular = peso_C + peso_H + peso_O % Calculo del porcentaje de oxigeno porcentaje_O = peso_O/peso_molecular*100



1-12

Introducción a Octave Para ejecutar este programa basta con escribir en la ventana de Octave su nombre sin la extensión “.m”. As´ı, ejecutar´ıamos el ejemplo anterior como se muestra a continuaci´ on7:

 pmol

peso_molecular = 178 porcentaje_O = 35.9551

Al ejecutar el programa, se ejecutan una tras otra las instrucciones contenidas en el fichero pmol.m, en el mismo orden en que aparecen en éste. Como era de esperar, sólo se muestran los resultados de las instrucciones no finalizadas con punto y coma. Las variables que se usan en el programa se quedan en el entorno de trabajo de on del programa. Por ello, debemos evitar utilizar en Octave cuando finaliza la ejecuci´ ficheros-M nombres de variables que coincidan con el nombre del propio programa, ya que en ese caso al realizar posteriores llamadas al programa se mostrar´ıa simplemente el valor de la variable. Es decir, si nuestro programa lo hemos guardado con el nombre pmol.m , debemos evitar el uso de una variable pmol. Ahora resulta evidente la utilidad de los ficheros-M a la hora de repetir conjuntos de órdenes. As´ı, si quisiéramos calcular el peso molecular de distintas moléculas, bastar´ıa con que cambi´ aramos los n´ umeros de cada tipo de átomos. Adem´ as, es aconsejable documentar los ficheros-M con comentarios, para que no olvidemos en el futuro qué pretend´ıamos al escribir el programa. Ejercicios 9 Escribe con editor de texto el programa pmol.m, y prueba a ejecutarlo.



10 Escribe y ejecuta un programa (ll´ amalo nitrogeno.m) que calcule la presión P que ejerce 1 mol de N2 en un volumen de 0.419 l. a 227C, según la ecuación de Van der Waals: 

P =

nRT V nb

2

n a − V 2 −

.

En esta ecuación, n es el n´ umero de moles de un gas, T la temperatura y V el volumen a las ·atm que se encuentra. R es la constante de los gases ideales (R = 0.082 litro mol·grad ). Las constantes del N2 son a = 1.390 l 2 atm/mol2 y b = 3.913e-02 l/mol.

1.6.

Instrucciones de entrada/salida

1.6.1.

Salidas por pantalla

En ocasiones, estaremos interesados en mostrar texto por pantalla, por ejemplo, para informar al usuario del uso del programa. Esto se puede conseguir con la orden disp(’texto’) . Si queremos mostrar el valor de una variable previamente calculada, podemos utilizar la misma orden incluyendo el nombre de la variable sin comillas:

 disp(’Este programa calcula...’) Este programa calcula... t = 5 ;  disp(t) 5

7

Tambi´ en podemos ejecutar el programa sin entrar en Octave: basta con ejecutar en un terminal Linux la orden octave -q pmol.m. ¿Por qué no pruebas a hacerlo?

Universitat Jaume I


1.6 Instrucciones de entrada/salida

1-13

La orden fprintf8 permite tener un mayor control sobre las salidas por pantalla que el que se tiene con disp, ya que permite especificar el formato con el que se van a mostrar los valores. La sintaxis de esta orden es la siguiente: fprintf (formato, expresiones)

El formato contiene el texto y las especificaciones de formato para las salidas, y va seguido de los nombres de las expresiones cuyo valor se desea visualizar separadas por comas. Dentro del formato se pueden usar los especificadores %e, %f y %g. Si se usa %e los valores se exhiben en una notación exponencial; si se usa %f se muestran en notaci´ on decimal; y si se usa %g, se usará el formato que sea más corto de los dos anteriores. La cadena n fuerza un cambio de l´ınea en la salida. Veamos un ejemplo:

\

 presion = 5.64;  temp = 245 ;  fprintf(’A %g atm. \n

la temperatura es

%f K’, presion, temp);

A 5.64 atm. la temperatura es 245.000000 K

F´ıjate en que se ha producido un cambio de linea a la mitad de la frase debido al especificador n. Además hay una variable por cada uno de los especificadores de formato: presion se corresponde con %g y temp con %f.

\

Los especificadores de formato tambi´ en pueden contener información para especificar el n´ umero de posiciones decimales que se muestran y las posiciones totales que se van a ocupar en la salida. Veamos un ejemplo para entender mejor el uso de esta funcionalidad: si ejecutamos la orden

 fprintf(’La temperatura es %6.2f K \n’, temp); lo que indicamos es que se reserven 6 posiciones para mostrar la variable temp, con 2 decimales. Las posiciones reservadas se ocupan (de derecha a izquierda) con los decimales, el punto, la parte no decimal y, en su caso, el signo. Si se reservan más posiciones de las necesarias, éstas se muestran como espacios en blanco a la izquierda. Veamos la salida de la orden anterior: La temperatura es 245.00 K

F´ıjate en qué ocurre si aumentamos el número de posiciones reservadas:

 fprintf(’La temperatura es La temperatura es

245.00 K

\

%8.2f K n’, temp);

En general, si sólo estamos interesados en indicar el número de decimales, podemos obviar la reserva de posiciones. 8

La orden original de Octave es printf, pero aqu´ı utilizaremos la orden fprintf porque, además de ser tambi´ en correcta, es compatible con Matlab.



1-14

Introducción a Octave

 fprintf(’La temperatura es %.4f K \n’, temp); La temperatura es 245.0000 K

Ejercicios b

11 Escribe las instrucciones Octave adecuadas para que se impriman las siguientes frases con el formato que se indica. Ten en cuenta que previamente se han definido las siguientes variables: a = 5 b = 48.56 c = 4.7864 d = 11111111111 

−

1. El valor de a es 5 2. El valor de a es 5.00

6. El valor de c es -4.7864 7. El valor de c es -4.8

3. El valor de b es 49 4. El valor de b es 48.56 5. El valor de b es 48.56000

1.6.2.

8. El valor de d es 1.111e+010 9. El valor de d es 1.e+010

Entrada de datos

Con lo visto hasta ahora, hemos podido realizar cálculos con distintas moléculas sin más que modificar el programa pmol.m. Sin embargo, ser´ıa más interesante que ni siquiera hubiera que modificar dicho fichero, sino que el propio programa nos “preguntara” los cambios que queremos hacer. Por ello, cualquier lenguaje de programación posee instrucciones que permiten al usuario introducir datos durante la ejecución del programa. En el caso de Octave disponemos de la orden input. Consideremos el siguiente fichero-M, modificación del archivo pmol.m y documentado con comentarios: %******************************************************************* % Programa : pmol.m % Descripcion: Este programa calcula el peso molecular de una % molecula organica. %******************************************************************* % Calculo de los g/mol de cada elemento. disp(’Programa para el calculo de pesos moleculares’); numero_C = input(’Introduce el numero de atomos de carbono: ’); numero_H = input(’Introduce el numero de atomos de hidrogeno: ’); numero_O = input(’Introduce el numero de atomos de oxigeno: ’); % Calculo de los pesos de los elementos peso_C = 12*numero_C ; peso_H = 1*numero_H ; peso_O = 16*numero_O; % Calculo del peso molecular peso_molecular = peso_C + peso_H + peso_O ; fprintf(’El peso molecular de la molecula es %g g/mol. peso_molecular);

Universitat Jaume I

\n’,

...


1.7 Ayuda de Octave

1-15

% Calculo del porcentaje de oxigeno porcentaje_O = peso_O/peso_molecular*100; fprintf(’El porcentaje en oxigeno de la molecula es %8.4f. n’, ... porcentaje_O);

\

Al ejecutar este programa, cada una de las instrucciones input mostrar´ıa en pantalla el mensaje que incluye, y esperar´ıa a que el usuario introdujera un dato y presionara la tecla Intro, almacenando dicho dato en la variable correspondiente:

 pmol

Programa para el calculo de pesos moleculares Introduce el numero de atomos de carbono: 9 Introduce el numero de atomos de hidrogeno: 6 Introduce el numero de atomos de oxigeno: 4 El peso molecular de la molecula es 178 g/mol. El porcentaje en oxigeno de la molecula es 35.955.

Ejercicios 12 Modifica el programa pmol.m para que incluya los cambios de entrada de datos e impresi´ on de resultados. Ejecuta dicho programa y comprueba que da los resultados esperados, aplic´ andolo sobre las siguientes moléculas orgánicas: 

a) Aspirina (C9 H4 O6 ). b) Benceno (C6 H6 ). c) Alcohol Et´ılico (CH3 CH2 OH). ´ d) Acido Acético (CH 3 COOH). e) Acetona (CH3 COCH3 ).

1.7.

Ayuda de Octave

Octave es

capaz de proporcionar ayuda al usuario mediante distintos medios. Por ejemplo, si queremos recordar la forma de uso de la orden rem (vista anteriormente), podemos obtener una respuesta rápida con la orden help rem:

 help rem

rem is the user-defined function from the file /usr/share/octave/2.1.40/m/general/rem.m - Mapping Function: rem (X, Y) Return the remainder of ‘X / Y’, computed using the expression x - y .* fix (x ./ y) An error message is printed if the dimensions of the arguments do not agree, or if either of the arguments is complex. See also: mod, round.



1-16

Introducción a Octave Ejercicios 13 Una funci´ on interesante que no hemos estudiado es la función diary. Busca ayuda sobre el funcionamiento de esta funci´ on. 

Universitat Jaume I


1.8 Ejercicios prácticos

1.8.

1-17

Ejercicios pr´ acticos

Es conveniente que pienses y realices los ejercicios que han aparecido a lo largo de la unidad marcados con el s´ımbolo b antes de acudir a la sesión de prácticas correspondiente. Deberás iniciar la sesión realizando los ejercicios marcados con el s´ımbolo . A continuaci´ on, deberás hacer el mayor número de los ejercicios siguientes. Ejercicios 14 Escribe un programa llamado OperBas.m que pida al usuario dos n´ umeros y calcule con ellos las operaciones básicas suma, resta, multiplicación y divisi´ on, mostrando los resultados en pantalla. 

15 Escribe un programa llamado farenheit.m que, dada una temperatura en grados Celsius, la convierta en grados Farenheit. La expresión matemática para realizar dicho cambio es: 9 F = C + 32 5 

16 La siguiente figura muestra una masa m en reposo sobre una superficie sin rozamiento. La masa está conectada a dos muros por muelles con constantes elásticas k 1 y k 2 . El periodo de este sistema viene dado por la expresión 

T = 2π

m k1 + k2



Escribe un programa Octave llamado muelles.m que pida al usuario los valores de m, k1 y k2 y que calcule y muestre el periodo T .

k1

k2 m

17 Escribe un programa llamado distancia.m que calcule la distancia entre dos puntos de un plano (x1 , y1 ) y (x2 , y2 ). 

distancia =



(x2

− x1)2 + (y2 − y1)2

Deberás pedir las coordenadas al usuario del programa. 18 Escribe un programa de nombre resistencia.m que calcule la resistencia equivalente de un circuito de tres resistencias en paralelo como el que se muestra en la siguiente figura. 

R1

R2

R3

La expresión de la resistencia equivalente viene dada por Req =


1 1 R1

+

1 R2

+

1 R3


1-18

Introducción a Octave 

19 Escribe un programa llamado GasIdeal.m que utilice la ecuación de los gases ideales P =

nRT V

de modo que pida al usuario el n´ umero de moles n de un gas, la temperatura T y volumen V ·atm a las que se encuentra, y calcule la presión P seg´ un esta ecuación (R = 0.082 litro mol·grad ). Ejemplo:

1 mol de N 2 que ocupa 0.419 l. a 227C, ejerce una presi´ on de 97.85 atm.

20 Repite el ejercicio anterior utilizando la ecuación de Van der Waals para gases imperfectos, nRT n2 a P = . V nb V 2 Aqu´ı deberás pedir al usuario, además los valores de a y b. El programa deberá llamarse VanDerWaals.m. 

− −

Ejemplo:

1 mol de N 2 que ocupa 0.419 l. a 227C, ejerce una presi´ on de 100 atm. Las constantes del N 2 son a = 1.390 l2 atm/mol2 y b = 3.913e-02 l/mol

21 La siguiente figura muestra la trayectoria de un proyectil lanzado a una velocidad v0 y un ángulo θ sobre un plano horizontal. 

y(t)

v0

Altura θ

x(t) Alcance

Suponiendo que pueden despreciarse todos los efectos de resistencia del aire, las ecuaciones del movimiento en los dos ejes vienen dadas por las ecuaciones siguientes: x(t) = v o cos(θ)t 1 2 y(t) = vo sin(θ)t gt 2

−

donde g = 9.81m/s2 . Teniendo en cuenta que la altura máxima del proyectil se corresponde con el punto donde dy/dt = 0, y que el movimiento es parabólico y, por lo tanto, simétrico, es relativamente fácil obtener la duración del vuelo del proyectil y su alcance. 2v0 sin(θ) g 2 v sin(2θ) Alcance = 0 g

Tiempo de vuelo =

Escribe un programa en Octave de nombre proyectil.m que calcule e imprima el tiempo de vuelo y alcance de un proyectil, pidiendo al usuario los valores de v0 y θ. Ejec´ utalo con distintos valores de ángulos entre 0 y 90 grados (de 5 en 5 grados, por ejemplo), comprobando que se obtiene el alcance máximo para θ = 45 y el tiempo de vuelo máximo para θ = 90. Ejemplo:

Universitat Jaume I

Para v 0 = 15m/s y θ = 45o , el tiempo de vuelo es de 2.16 s. y el alcance 22.94 m.


1.8 Ejercicios prácticos 

1-19

22 En la siguiente figura se muestra una versi´ on simplificada de un radio receptor AM. Antena

L C +

V

R

V

0

R

−

Tierra

El voltaje que atraviesa la resistencia V R var´ıa en funci´ on de la frecuencia seg´ un la ecuación V R =

R

 

R2 + ωL

−

1 ωC



2

V 0

donde ω = 2πf , y f es la frecuencia en hercios. Escribe un programa que pida al usuario los valores de R, V 0 , f , L y C , y calcule el voltaje V R . La salida del programa deberá ser lo más parecida a la que se muestra a continuación. % Ejemplo de uso Voltaje_antena Introduce el valor de la resistencia R (en ohms): 50 Introduce el valor del voltaje V0 (en voltios): 0.01 Introduce el valor de la frecuencia f (en Hertzios): 1E6 Introduce el valor de la inductancia L (en henries): 0.0001 Introduce el valor de la capacitancia C (en faradays): 0.25E-9 El valor de VR es 0.00986496 voltios.





Tema 3

Representaci´ on gr´ afica 2D y 3D Guillermo Peris Ripoll´ es Aplicaci´ on Optimizaci´ on del flujo de un canal de agua.

A la hora de dise˜ nar un canal de agua, como el que se muestra en la siguiente figura, se pretende que la velocidad de flujo sea máxima.

h

θ

θ

c

Las variables a considerar a la hora de diseñar el canal son su altura h, la anchura de la base c, y el ángulo lateral θ. Puede demostrarse que la velocidad de flujo es inversamente proporcional al per´ımetro, cuya expresión es: p = c +

2h sin(θ)

mientras que el área total de flujo es A = ch + h2 cot(θ). Considerando que el canal debe tener 200 cms. de área transversal (A = 200), la funci´ on 1/p u ńicamente dependerá de h y θ. Representa la superf´ıcie 1/p = f (h, θ), mostrando el diagrama de contornos simultáneamente. ¿Para qué valores aproximados de h y θ la velocidad de flujo es máxima? Velocidad de flujo en un canal

1/p

0.03 0.025 0.02 0.015 0.01

20 0


15 20

40 60 80 Theta (Grados) 100

10 120

140

5 160

Altura (cm)

0


3-2

Representación gráfica 2D y 3D

Contenidos

3.1. Introducci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-3 3.2. Representaciones 2D . . . . . . . . . . . . . . . . . . . . . . 3-3

3.2.1. Gráficas simples . . . . . . . . . . . . . . . . . . . . . . . . . 3-3 3.2.2. Gráficas m´ ultiples . . . . . . . . . . . . . . . . . . . . . . . . 3-4 3.2.3. Mejorando el aspecto de las gr´ aficas . . . . . . . . . . . . . . 3-6 3.3. Representaciones 3D: l´ıneas, contornos y superf´ıcies. . . . 3-7

3.3.1. L´ıneas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-7 3.3.2. Superficies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-8 3.3.3. Contornos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-10 3.4. Aplicaci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-11 3.5. Ejercicios pr´ acticos

Universitat Jaume I

. . . . . . . . . . . . . . . . . . . . . . . 3-13


3.1 Introducción

3.1.

3-3

Introducci´ on

Resulta muy habitual que los ingenieros utilicen gráficos para mostrar sus ideas de una forma m´ as clara, ya que es más sencillo identificar tendencias en una figura que una tabla de resultados. Octave dispone (junto con el paquete octave-forge) de un gran conjunto de funciones u ´ tiles para la creación de gráficos. En este tema estudiaremos algunas de ellas.

3.2.

Representaciones 2D

Todas las funciones de que dispone Octave para la creació n de gráficos utilizan el programa gnuplot. Este es un programa que podemos usar de forma independiente de Octave, aunque aqu´ı aprenderemos a utilizarlo desde la interfaz de Octave.

3.2.1.

Gr´ aficas simples

Para dibujar gráficas, Octave dispone de la orden plot(x,y), donde x e y son dos vectores de la misma dimensión que representan las coordenadas de las abscisas y ordenadas de los datos a representar, respectivamente. Supongamos que queremos representar la gráfica de sin(x) entre 0 y 2π. Entonces, deber´ıamos crear en primer lugar un vector con valores de x, y el vector sin(x). Para ello, podemos utilizar por ejemplo 100 puntos equiespaciados en el eje x y calcular su seno:  x = linspace(0,2*pi,100);  y = sin(x) ;  plot(x,y) ;

Al ejecutar este código, nos aparece una ventana similar a la siguiente: 1 line 1 0.8

0.6

0.4

0.2

0

−0.2

−0.4

−0.6

−0.8

−1 0

1

2

3

4

5

6

7

afica simple. Figura 3.1: Gr´

Pensemos que cuantos más puntos utilicemos, más fiel será la representación obtenida, pero más tiempo le costará a Octave obtener el resultado.



3-4

Representación gráfica 2D y 3D No es demasiado complicado mejorar esta gráfica. Por ejemplo, con la orden title se a˜ nade un t´ıtulo, con las o´rdenes xlabel e ylabel se a˜ naden, respectivamente, etiquetas a los ejes de abscisas y ordenadas; además, la orden grid(’on’) a˜ nade una rejilla a la gráfica. Estas o´rdenes no surten efecto hasta que no se actualiza la gráfica con la orden replot. Observa el uso de estas órdenes en el siguiente ejemplo, cuyo resultado se muestra en la Figura 3.2:  title(’Representacion de sin(x)’), xlabel(’x’), ...

ylabel(’sin(x)’), grid(’on’), replot

F´ıjate en que los textos deben escribirse entre comillas simples. Representacion de sin(x) 1

line 1

0.8 0.6 0.4 0.2 ) x ( n i

0

s

−0.2 −0.4 −0.6 −0.8 −1 0

1

2

3

4

5

6

7

x

Figura 3.2: Grafica simple con t´ıtulo, rejilla y etiquetas.

3.2.2.

Gr´ aficas m´ ultiples

Si estamos interesados en dibujar varias curvas en una misma gráfica, podemos hacerlo de forma simple utilizando varios pares de vectores. Por ejemplo, la orden plot(x,y,w,z) dibujar´ıa en una misma gr´ afica las curvas de y en funció n de x, y w en funci´ on de z. F´ıjate en el siguiente código, que dibuja las curvas sin(x) y cos(x) en la misma figura (Figura 3.3):  x = linspace(0,2*pi,100);  y = sin(x) ; z = cos(x);  plot(x,y,x,z),title(’Grafica del seno y coseno de x’), ...

xlabel(’x’), grid(’on’) , replot Octave distingue

las dos l´ıneas utilizando colores, lo cual no se puede apreciar en impresiones en blanco y negro (como la que tienes entre tus manos ;) ). Ser´ıa interesante poder distinguir las curvas bien por el formato de las l´ıneas, bien por el formato de los puntos. Para ello, debemos indicar en la orden plot tripletes de datos, a saber: el vector x, el vector y y un formato de lineas y puntos, teniendo en cuenta las opciones que se muestran en la Tabla 3.1. Si antecedemos un guión al indicador del Universitat Jaume I


3.2 Representaciones 2D

3-5 Grafica del seno y coseno de x

1

line 1 line 2

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

1

2

3

4

5

6

7

x

afica m´ ultiple. Figura 3.3: Gr´

punto, éstos se unen con una l´ınea cont´ınua 1 . Por ejemplo, una mejor representación de la curva anterior se har´ıa con la siguiente orden (Figura 3.4):  plot(x,y,’-o’,x,z,’+’),...

title(’Gr´ afica del seno y coseno de x’),xlabel(’x’), grid(’on’) , replot

Grafica del seno y coseno de x 1

line 1 line 2

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

1

2

3

4

5

6

7

x

afica con marcas. Figura 3.4: Gr´

Hemos comprobado que, al crear una nueva figura, esta aparece en una ventana de figuras, desapareciendo la figura anterior (si la hab´ıa). Si estamos interesados en que aparezcan distintas figuras (en distintas ventanas) debemos ejecutar previamente la orden figure. Por u ´ ltimo, la orden clg limpia la ventana de gráficos. Esta orden resulta conveniente utilizarla antes de crear una nueva figura. 1

Consulta el manual de Octave para obtener m´ as informaci´ on sobre el dise˜ n o de gr´ aficas.



3-6

Representación gráfica 2D y 3D Punto punto más estrella c´ırculo marca

Indicador . + * o x

Tabla 3.1: Opciones de marcas.

Ejercicios 1 Ejecuta las l´ıneas de código necesarias en Octave para obtener las gráficas de las figuras anteriores. 

on y = sin(x) + x − x · cos(x) en el intervalo x ∈]0, 30[. A˜ nade un t´ıtulo 2 Dibuja la funci´ y las etiquetas que creas oportunas. 

afica, es importante el muestreo que se realice de la función. 3 A la hora de realizar una gr´ Para entender mejor ésto, compara los siguientes ejemplos y piensa en cual es la causa de que la segunda gráfica presente un mejor aspecto que la primera. 

 n = 5;

 n = 25;

 x = 0:1/n:3;

 x = 0:1/n:3;

 y = sin(5*x) ;

 y = sin(5*x) ;

 plot(x,y)

 plot(x,y)

3.2.3.

Mejorando el aspecto de las gr´ aficas

En este apartado vamos a introducir algunas órdenes más para cambiar el aspecto de los gráficos. Por ejemplo, podemos estar interesados en personalizar la leyenda explicativa de las distintas l´ınea de una gráfica. Esta leyenda normalmente aparece con los textos line 1, line 2, etc. Para cambiar este texto, debemos escribirlo entre punto y comas (;) del formato de l´ınea y punto de la orden plot. Esto se entiende mejor con el anterior ejemplo de la representación de seno y coseno:  plot(x,y,’-o;seno(x);’,x,z,’.;coseno(x);’)  title(’Gr´ a fica del seno y coseno de x’),xlabel(’x’)  grid(’on’) , replot

En ocasiones resulta interesante tener distintos gráficos en una misma ventana. Esto es lo que se conoce como subventanas. Para ello, podemos dividir la ventana de gráficos en dos subventanas (apiladas o contiguas) o cuatro subventanas, utilizando la orden subplot. Esta orden acepta 3 argumentos, donde los dos primeros indican las divisiones en filas y columnas de la ventana, y el tercero a la posición en la que se va a situar la siguiente gráfica (las ventanas se numeran de izquierda a derecha, y de arriba y abajo). La utilización de subgráficas la entenderás mejor con el siguiente ejercicio.

Universitat Jaume I


3.3 Representaciones 3D: l´ıneas, contornos y superf´ıcies.

3-7

Gráfica del seno y coseno de x 1

seno(x) coseno(x)

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0

1

2

3

4

5

6

7

x

afica con marcas. Figura 3.5: Gr´

Ejercicios 

4

Ejecuta las siguientes lineas de c´ odigo y trata de entender por ti mismo la orden subplot.

 x = linspace(0,2*pi,100);  y = sin(x) ; z=cos(x);  subplot(2,1,1)  plot(x,y,’;seno(x);’)  subplot(2,1,2)  plot(x,z,’;coseno(x);’)

Hasta ahora hemos observado que Octave fija automáticamente la escala de los ejes, ajustándola a los valores de los datos. Podemos cambiar estos valores ejecutando la orden axis([xmin xmax ymin ymax]), donde se indican respectivamente los valores m´ınimos y máximo de x, y m´ınimo y máximo de y.

3.3.

Representaciones 3D: l´ıneas, contornos y superf´ıcies.

Hasta ahora hemos obtenido gráficas en dos dimensiones a partir de funciones dependientes de una sola variables. Si las funciones dependen de dos variables (o más) se necesitan representaciones de orden superior. Para ello, Octave incluye distintas funciones para la representaci´ on tridimensional. A continuación vamos a explicar las funciones b´ asicas para crear tres tipos de gráficas: l´ıneas 3D, superficies y contornos.

3.3.1.

L´ıneas

Se pueden representar l´ıneas en el espacio tridimensional con la orden plot3(x,y,z) de Octave, donde x, y y z son funciones de un parámetro. Por ejemplo, las siguientes Ingenier´ıa Qu´ımica


3-8

Representación gráfica 2D y 3D ecuaciones generan una curva en tres dimensiones a medida que var´ıa el parámetro t: x = e−0.05tsin(t) y = e−0.05tcos(t) z = t Esta curva se puede representar con las siguientes órdenes Octave, dando lugar a la figura que se muestra a continuación:  t = [0:pi/50:10*pi];  plot3(exp(-0.05*t).*sin(t), exp(-0.05*t).*cos(t), t)  xlabel(’x’), ylabel(’y’), zlabel(’z’), grid(’on’), replot

line 1

z

35 30 25 20 15 10 5 0

−0.8

−0.6

−0.4

−0.2 x

0

0.2

0.4

0.6

0.8

1 −1

0 −0.2 −0.4 −0.6 −0.8

1 0.8 0.6 0.4 0.2 y

Figura 3.6: Curva 3D.

Ejercicios 

5

Ejecuta las l´ıneas de código anteriores en Octave para obtener la gráfica de la Figura 3.6.

3.3.2.

Superficies

La función z = f (x, y) representa una superficie en un sistema de coordenadas xyz. Antes de realizar la representación, es necesario crear una malla de puntos en el plano xy para calcula el valor de z en cada uno de ellos. Para ello, Octave dispone de la funci´ on meshgrid. La sintaxis de esta orden es: [X, Y] = meshgrid(x,y)

donde x e y son vectores con los valores de esta variables. X es una matriz en la que el vector x se copia en cada una de sus filas, e Y es una matriz en la que el vector y se copia en cada una de sus columnas. De esta forma, podemos trabajar con las matrices X e Y para obtener una matriz Z en t´ erminos de la función representada. El uso de meshgrid se ilustra con el siguiente ejemplo sencillo: Universitat Jaume I


3.3 Representaciones 3D: l´ıneas, contornos y superf´ıcies.

3-9

 x = [1:5] ; y = [6:10];  [X, Y] = meshgrid(x,y); X

X 1 1 1 1 1

= 2 3 2 3 2 3 2 3 2 3 Y Y = 6 6 6 7 7 7 8 8 8 9 9 9 10 10

4 4 4 4 4

5 5 5 5 5

6 6 7 7 8 8 9 9 10 10 10

Una vez calculada la malla de valores X e Y , y con ellos la matriz Z , utilizaremos la funci´ on mesh(X,Y,Z) para crear la superficie. Veamos un ejemplo en el que se genera la superficie de la función z = xe −[(x−y ) +y ] para −2 ≤ x ≤ 2 y −2 ≤ y ≤ 2 con un espaciado de 0.1. 2 2

2

 x = [-2:0.1:2]; y = x;  [X, Y] = meshgrid(x,y);  Z = X.*exp(-((X-Y.∧ 2).∧ 2+Y.∧ 2));  mesh(X,Y,Z), xlabel(’x’), ylabel(’y’), zlabel(’z’) line 1

z

0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 2 1.5 1 −2

0.5 −1.5

−1

−0.5 x

0 −0.5 0

0.5

1

1.5

y

−1 −1.5 2 −2

Figura 3.7: Superficie 3D.

Ejercicios 

6

Ejecuta las l´ıneas de código anteriores en Octave para obtener la gráfica de la Figura 3.7.



3-10


3.3.3.

Contornos

Las representaciones topográficas muestran los contornos terrestres mediante l´ıneas de altura constante. Estas lineas se conocen como lineas de contorno. Si caminamos a lo largo de una de estas lineas, la altura se mantiene constante. Octave dispone de la función contour(X,Y,X), que se utiliza de una foram similar a mesh:  x = [-2:0.1:2]; y = x;  [X, Y] = meshgrid(x,y);  Z = X.*exp(-((X-Y.∧ 2).∧ 2+Y.∧ 2));  contour(X,Y,Z), xlabel(’x’), ylabel(’y’), zlabel(’z’) 0.4 0.3 0.2 0.1 −2.78e−17 2 −0.1 −0.2 1.5 −0.3 −0.4 1 0.5 0

y

−0.5 −1 −1.5 −2 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

x

Figura 3.8: Contorno.

Las representaciones de contorno y superficie pueden combinarse con la función meshc(X,Y,Z).  meshc(X,Y,Z), xlabel(’x’), ylabel(’y’), zlabel(’z’)

z

0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5 2 1.5 1 −2

0.5 −1.5

−1

−0.5 x

0 −0.5 0

0.5

1

1.5

y

−1 −1.5 2 −2

Figura 3.9: Contorno+superficie.

Ejercicios odigo necesarias en Octave para obtener las gráficas de las Figu7 Ejecuta las l´ıneas de c´ ras 3.8 y 3.9. 

Universitat Jaume I


3.4 Aplicación

3.4.

3-11

Aplicaci´ on

A continuaci´ on se muestra el programa que resuelve la aplicación expuesta al principio del tema. Las únicas novedades que introduce son el uso de dos nuevas funciones: axis: Esta funci´ on recibe como argumento un vector de 4 o 6 componentes, en

funci´ on de que la gráfica sea bidimensional o tridimensional. Por ejemplo, en el presente caso los 6 elementos del vector ser´ıan los valores m´ınimos y máximos de x, seguidos de los valores m´ınimos y máximos de y , y por último los valores m´ınimos y máximos de z. La gráfica utilizará estos valores para los l´ımites de cada uno de los ejes. view: Esta funci´ on permite cambiar el punto de vista de una función 3D. Para entender mejor su uso, ejecuta la ayuda de Octave ( help view).

Por lo demás, no presenta demasiadas complicaciones, as´ı que trata de entenderla por ti mismo, y ejecútala en Octave para comprobar que obtienes el mismo resultado. % Limpiamos variables y graficos (si existen). clear clg % Definicion de constantes y vectores iniciales. A=200; Theta = 20:10:150; h = 2:20; % Creamos la malla de valores (x,y) y la funcion p. [alt,Th] = meshgrid(h, Theta); c = (A-alt.∧2./tan(Th*pi/180))./alt; p = c + (2*alt./sin(Th*pi/180)); % Realizamos la representacion grafica. meshc(Th,alt, 1./p) title(’Velocidad de flujo en un canal’) xlabel(’Theta (Grados)’), ylabel(’Altura (cm)’), zlabel(’1/p’) axis([0 160 0 20 0.008 0.03]) view(30,15), replot



3-12


Velocidad de flujo en un canal

1/p

0.03 0.025 0.02 0.015 0.01

0

20

40 60 80 Theta (Grados) 100

120

140

160

0

5

10

15

20

Altura (cm)

on del flujo de un canal de agua. Figura 3.10: Optimizaci´

Universitat Jaume I



3.5.

3-13


Es conveniente que pienses y realices los ejercicios que han aparecido a lo largo de la unidad marcados con el s´ımbolo b antes de acudir a la sesión de prácticas correspondiente. Deberás iniciar la sesión realizando los ejercicios marcados con el s´ımbolo . A continuaci´ on, deberás hacer el mayor número de los ejercicios siguientes. Ejercicios 

8

Consideremos de nuevo la ecuaci´ on de Van der Waals para un mol de gas: P =

RT V − b

−

a V 2

Escribe un programa isotermas.m que pida al usuario los valores de a y b, y represente una gr´ afica con las isotermas del gas (P en función de V, para un valor constante de T) a 100, 200, 300 y 400 grados cent´ıgrados. Recuerda que en la ecuación de Van der Waals la temperatura debe expresarse en Kelvin. Cada curva debe ir con un trazo diferenciado, con el texto que indique la isoterma que se ha representado, as´ı como el t´ıtulo de la gráfica y la etiqueta de los ejes. Comprueba el funcionamiento correcto del programa con el benceno, para el cual a = 18.78 atm · l/mol2 y b = 0.1208 l/mol. Ecuacion de Van der Waals: Isotermas 25 T=100ºC T=200ºC T=300ºC T=400ºC 20

15

) m t a ( n o i s e r P

10

5

0 0


10

20

30

40

50 Volumen (l)

60

70

80

90

100


3-14

Representación gráfica 2D y 3D on y se someten a cambios de presión 9 Se consideran 4.003 g. de helio y 39.944 g. de arg´ a la temperatura de 0 grados Celsius, obteniéndose las siguientes tablas de valores: 

´ a 0o ARGON P(atm) V(l) 1.000 22.4 11.10 2.000 32.79 0.667 43.34 0.500 53.68 0.400 63.68 0.333

HELIO a 0o P(atm) V(l) 1.002 22.37 0.8067 27.78 0.6847 32.73 0.5387 41.61 0.3550 63.10 0.1937 115.65

Representa gráficamente el volumen frente a la presión en ambos gases. Comprueba que el helio es un gas que verifica la ley de Boyle-Mariotte: P × V = constante, pero que el argón no cumple exactamente esta ley. Para ello, deberás representar el producto P V para cada gas. El programa se llamará mariotte.m. 25

120

Argon a 0ºC

Helio a 0ºC 110

20

100

90

15

80

) l ( n e m u l o V

) l ( n e m 70 u l o V

10

60

50

5

40

30

20 0.1

0.2

0.3

0.4

0.5

0.6 0.7 Presion(atm)

0.8

0.9

1

1.1

0 0

20

40

60 Presion(atm)

80

100

120

40 Argon a 0ºC: P*V

22.416 Helio a 0ºC: P*V

38 22.414

36 22.412

34

22.41

32 V * 30 P

V *22.408 P

28 22.406

26 22.404

24 22.402

22

22.4 1

Universitat Jaume I

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

20

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6



3-15

on de la molécula de CO 2 consiste en la extensión 10 Uno de los modos normales de vibraci´ lineal de los enlaces carbono-ox´ıgeno. 

A

B

C

O

C

O

R

R

AB

BC

Podemos obtener una aproximación a la energ´ıa de vibración de este modo normal mediante el potencial de Morse, según la siguiente ecuación: E vib = V AB (RAB ) + V BC (RBC ) , donde



−2AXY (R−BXY )

V XY (R) = d XY e

−AXY (R−BXY )

− 2e



.

En esta ecuación R XY es la distancia entre los átomos X e Y , y el resto de parámetros para ˚ el CO2 son: dAB = d BC = 7.65eV , b AB = b BC = 1.162A. Variando los valores de RAB y RBC , obtén la superf´ıcie de energ´ıa potencial para este modo normal de vibraci´ on del CO2 , y la superf´ıcie de contorno asociado. Superficie de energia potencial del CO2

Energia (u.a.)

35 30 25 20 15 10 5 0 -5 -10 -15 -20 2.6 2.4 2.2

0.6 0.8


1

1.2 1.4 1.6 1.8 RAB (Angstrom)

2

1 0.8 2.2 2.4 0.6 2.6

2 1.8 1.6 1.4 RBC (Angstrom) 1.2


Tema 5

Instrucciones repetitivas Guillermo Peris Ripoll´ es

Objetivos Cuando finalice este tema, el alumno deberá ser capaz de: Entender el funcionamiento de los bucles for, while y do-until, utilizándolos en las situaciones que corresponda de forma adecuada. Comprender el anidamiento de bucles, aplicándolo en los problemas en que resulten convenientes.

Aplicaci´ on Cuando finalice este tema, el alumno deberá ser capaz de resolver problemas como el siguiente, cuya resolución se indica a lo largo del propio tema. C´ alculo de la energ´ıa de repulsi´ on nuclear en una mol´ ecula

La energ´ıa de repulsión nuclear entre dos átomos i y j de una molécula viene dada por la expresión Z i Z j E ij = , Rij

·

donde Z i es la carga nuclear del átomo i y Rij es la distancia que separa los átomos i y j . 1. Escribe un programa en octave que pida al usuario un vector de cargas nucleares Z y una matriz de distancias R ij , y con ellos calcule la matriz E ij de energ´ıas nucleares según la ecuación anterior. Deberás tener en cuenta que los elementos de la diagonal de la matriz E ij deben ser infinito (ya que la distancia entre un átomo y s´ı mismo es 0). 2. Añade al programa anterior las l´ıneas necesarias para calcular la energ´ıa total de repulsi´ on. 3. Añade al programa anterior las l´ıneas necesarias para decidir qué átomos experimentan una mayor repulsión.



5-2

Instrucciones repetitivas

Contenidos

5.1. Introducci´ on 5.2. Bucles for 5.3. Bucles while

. . . . . . . . . . . . . . . . . . . . . . . . . . . 5-3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-6

5.4. Bucles do-until

. . . . . . . . . . . . . . . . . . . . . . . . . 5-7

5.5. Selecci´ on del tipo de bucle . . . . . . . . . . . . . . . . . . . 5-8 5.6. Bucles anidados . . . . . . . . . . . . . . . . . . . . . . . . . 5-9 5.7. Aplicaci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5-11 5.8. Ejercicios pr´ acticos

Universitat Jaume I

. . . . . . . . . . . . . . . . . . . . . . . 5-14


5.1 Introducción

5.1.

5-3

Introducci´ on

Un bucle es una estructura de programación que permite la repetición controlada de un conjunto de instrucciones. Este tipo de estructuras, y en particular las instrucciones for y while , son utilizadas de forma generalizada en la inmensa mayor´ıa de lenguajes de programaci´ on, lo cual hace interesante su estudio. Sin embargo, en muchos casos los bucles pueden (y deben) evitarse en Octave, ya que hay estructuras t´ıpicas de este lenguaje m´ as eficaces que los propios bucles. No obstante, s´ı que existen ocasiones en las que es necesario el uso de bucles, por lo que estudiaremos estas instrucciones con detalle en esta unidad.

5.2.

Bucles for

Los bucles for se utilizan cuando nos interesa repetir un bloque de instrucciones un n´ umero predeterminado de veces. La estructura general del bucle for es la siguiente 1 : for i = vector instrucciones end

El conjunto de instrucciones se repite para cada elemento de vector , denominándose iteraci´ on cada una de estas repeticiones. En cada iteración, i toma ordenadamente el valor de cada elemento de vector . Al igual que ocurr´ıa con las instrucciones if , resulta conveniente realizar un sangrado de l´ınea en las instrucciones del interior del bucle. Veamos un ejemplo sencillo:

 for i = [1:5] j = 2*i; disp(j)

end 2 4 6 8 10

En este ejemplo, observamos como el valor de i toma en cada iteración, y de forma ordenada, los valores del vector [1 2 3 4 5], y dichos valores se utilizan para realizar cálculos en el interior del bucle. El operador ’:’ para la creación de listas es ampliamente utilizado en los bucles for. En este caso, los bucles tendrán la siguiente forma particular: for k = [inicial:incremento:final] instrucciones end

aunque los corchetes no son realmente necesarios. Por ejemplo, supongamos que queremos calcular los cuadrados de los primeros 5 números impares. Podr´ıamos hacerlo utilizando una instrucción for... 1

En lugar de recorrer un vector en cada iteració n del bucle, podemos recorrer una matriz. En ese caso, en cada ciclo i es una de las columnas de la matriz.



5-4


 for i = 1:2:9 j = i∧ 2; disp(j)

end 1 9 25 49 81

...aunque en este caso ser´ıa más eficiente no hacerlo: ∧

 j = [1:2:9]. 2;  disp(j’) 1 9 25 49 81

Dentro de un bucle podemos modificar el valor de la variable que recorre el vector, aunque al terminar la iteración seguirá tomando el siguiente valor que le corresponda. En el siguiente ejemplo, intentamos terminar el bucle cambiando el valor de la variable i, sin ning´ un efecto en el resultado final:

 for i = 1:2:9 j = i∧ 2; disp(j) i = 10;

end 1 9 25 49 81

De todas formas, para evitar posibles errores no conviene cambiar el valor de la variable de iteraci´ on en el interior del bucle. Veamos un par de ejemplos de bucles for para entender mejor el funcionamiento

de estas estructuras. Supongamos que queremos escribir un programa que sume los elementos de un vector (lo cual hace la función sum, pero nos olvidaremos momentáneamente de eso). Para ello, tomamos el vector y sumamos a una variable suma cada uno de sus elementos:

 v = 1:2:9;  for x = v

suma = suma + x;

end error: ’suma’ undefined...

Universitat Jaume I


5.2 Bucles

5-5

for

¿Qu´ e ha pasado? El c´ odigo parece correcto: sumamos todos los elementos a la variable suma, pero aparece un mensaje de error. La causa es que no hemos inicializado la variable suma: en la primera iteración, cuando v = 1, tenemos la instrucción suma = suma + 1 pero, ¿cu´ al es el valor inicial de suma al que sumarle 1? Para obtener el resultado esperado, debemos definir suma con un valor inicial de 0, como se muestra en el código correcto:

 v = 1:2:9;  suma = 0 ;  for x = v

suma = suma + x;

end disp(suma); 25



Supongamos ahora que nos interesa calcular la media de los valores de un vector. En ese caso, debemos dividir la suma de los elementos del vector por el número de elementos. Para ello, debemos contar el n´ umero de elementos, para lo cual utilizaremos un contador , que inicializaremos a 0 antes del bucle (aún no hemos tomado ning´ un valor del vector) e incrementaremos en 1 cada vez que encontremos un nuevo elemento. Finalmente, dividiremos la suma de los elementos por el número de éstos para obtener as´ı la media. Fijémonos en que esta división se realiza al terminar el bucle, y no en su interior. El código Octave se reproduce a continuación:

 v = 1:2:9;  suma = 0 ;  contador =  for x = v

0 ;

suma = suma + x; contador = contador + 1;

end media = suma/contador; disp(media); 5

 

Esta técnica general se utiliza en muchos algoritmos en los que se necesita sumar, contar elementos, o calcular medias. Octave puede realizar este cálculo m´ as eficientemente. De hecho, en este caso, no hace falta utilizar un contador, ya que la orden umero de elementos de v y, por lo tanto, el n´ umero de itelength(v) nos dar´ıa el n´ raciones. Má s a´ un, la función mean calcula con una sola instrucció n la media de un vector:

 mean(1:2:9) ans = 5

Por supuesto, internamente la función mean recorre los elementos del vector sumándolos, pero lo hace a mayor velocidad que el bucle for.



5-6

Instrucciones repetitivas Ejercicios b

b

Escribe un programa que calcule el producto de los elementos del vector [1:2:9] haciendo uso de un bucle for y sin hacer uso de las funciones length y prod. 1



Modifica el programa anterior para que calcule la media geom´ etrica de los elementos de un vector utilizando bucles y contadores, es decir, sin hacer uso de funciones predefinidas como length, prod, etc. Recuerda que la media geométrica de una serie de n valores X 1 X 2 X n se calcula como Media geom. = X 1 X 2 X n y que x x1/n . 2



√ n

5.3.

· ··

√ ≡ n

· ··

Bucles while

El bucle while es una estructura que se utiliza para repetir un conjunto de instrucciones mientras se cumpla una condició n lógica determinada. La estructura general de este bucle es la siguiente: while condici´ on instrucciones end

Mientras (en inglés, while ) la condici´ on es verdadera, se ejecutan las instrucciones, tras lo cual se vuelve a comprobar la condici´ on . En el momento en que ésta es falsa, se termina el bucle. Fijémonos en que alguna de las variables de condici´ on debe cambiar durante las instrucciones, de lo contrario el valor de condici´ on ser´ıa siempre el mismo y el bucle entrar´ıa en un ciclo infinito. Veamos un ejemplo sencillo de bucle while: i = 0 ; while i < 3 disp(i) ; i = i + 1; end disp(’Terminado’);

La ejecución de este código da lugar a la siguiente salida: 0 1 2 Terminado

Analicemos este ejemplo con un poco más de detalle: en primer lugar, se inicializa el valor de i a 0. Como se cumple la condición del bucle (i<3) se ejecutan sus instrucciones, mostrándose el valor 0 y aumentando en uno el valor de i. Se vuelve a comprobar la expresión del bucle, que vuelve a ser verdadera, por lo que se muestra el valor 1 e i pasa a valer 2. Se repite este último paso, se muestra el valor 2 e i pasa a valer 3. Pero ahora, al no cumplirse la condición del bucle, se termina el ciclo, pasando a ejecutarse la instrucción que imprime la palabra Terminado. Para las condiciones de los bucles while se pueden utilizar todos los operadores y funciones lógicas vistas en la unidad anterior. Fijémonos en que los bucle for se van a utilizar cuando se quiera repetir un conjunto de intrucciones un número predeterminado de veces, mientras que en el bucle while se busca el cumplimiento de una condición para la finalización del ciclo. Universitat Jaume I


5.4 Bucles do-until

5-7

Ejercicios b



Determina la salida de los siguientes c´ odigos:

3

a)

i = 0 ; while i <= 3 disp(i) ; i = i + 1; end disp(’Terminado’);

b)

i = 0 ; while i < 10 disp(i) ; i = i + 2; end disp(’Terminado’);

5.4.

c)

i = 3 ; while i < 10 disp(i) ; i = i + 2; end disp(’Terminado’);

d)

i = 1 ; while i < 100 i = i * 2; disp(i) ; end disp(’Terminado’);

Bucles do-until

Los bucles do-until son muy similares a los bucles while, en lo que respecta a que la finalizaci´ on del bucle está ligada al cumplimiento de una condición. Sin embargo, existen dos diferencias importantes: en los bucles do-until la condición se comprueba al final de la estructura, y se ejecutan las órdenes hasta que se cumple la condición (y no mientras, como ocurre con while). La estructura general de este bucle es la siguiente: do instrucciones until condici´ on

Fij´ emonos en una diferencia importante que surge como consecuencia de evaluar la condición al final: las instrucciones del interior del bucle se ejecutan como m´ınimo una vez . Veamos un ejemplo equivalente al visto anteriormente para while, en el que se imprim´ıan los números del 0 al 2: i = 0 ; do disp(i) ; i = i + 1; until i > 2 disp(’Terminado’);

Ejercicios b

Reescribe las l´ıneas de código del ejercicio anterior con bucles do-until, de forma que muestren la misma salida que aquellos. 

4



5-8


5.5.

Selecci´ on del tipo de bucle

En muchas ocasiones, intuimos que debemos utilizar algún tipo de bucle, pero no sabemos con certeza cuál de ellos es el más adecuado. Antes de ver ejemplos que ilustran las diferencias prácticas entre los tres bucles estudiados, recordemos brevemente sus caracter´ısticas más significativas: Bucle for: Repetición de un conjunto de instrucciones un n´ umero predeterminado de veces. Bucle while: Incumplimiento de una condición al inicio del bucle para la finalización del ciclo. Las instrucciones del bucle se ejecutan 0 o más veces. Bucle do-until: Cumplimiento de una condición al final del bucle para la finalizaci´ on del ciclo. Las instrucciones del bucle se ejecutan 1 o más veces. Vamos a aplicar los tres bucles a la implementación del desarrollo en serie de funciones. Por ejemplo, supongamos que queremos calcular 1 /(1 x) con un desarrollo de Taylor de la función:

−

1 1

−x

∞

x =1 + x + x + x + x + ... = i

2

3

4

.

i=0

Si quisiéramos calcular un número determinado de términos de la serie (por ejemplo 8) utilizar´ıamos un bucle for, como se muestra en el siguiente código con x =0.6:

 x = 0.6;  suma = 0 ;  for i = 0:7

suma = suma + x∧ i;

end disp(suma) ; suma = 2.4580 disp(1/(1-x)) 2.5000

 

Pero, ¿y si nos interesa obtener el resultado con una cierta precisió n? En ese caso, debemos comprobar que la diferencia entre los dos últimos valores calculados sea menor que un cierto valor pequeño, digamos 0.00001. En ese caso resulta más conveniente la utilización de un bucle while o do-until: x

= 0.6;  precision = 0.00001 ;  suma = 1 ;  termino = 1 ; % Inicializacion  while abs(termino) > precision termino = termino*x; suma = suma + termino; endwhile  suma suma = 2.5000

x

= 0.6;  precision = 0.00001 ;  suma = 1 ;  termino = 1 ; % Inicializacion  do termino = termino*x; suma = suma + termino; until abs(termino) < precision  suma suma = 2.5000

Varios comentarios al respecto de este código. En primer lugar, utilizamos un producto en lugar de una potencia para el cálculo de cada término 2 de la serie, 2

¿Por qu´ e se toma el valor absoluto del valor de termino en la condici´ on del while?

Universitat Jaume I


5.6 Bucles anidados

5-9

calculándolo multiplicando por x el término anterior. Utilizamos un contador i para comprobar el n´ umero de iteraciones necesarias para converger a la precisión requerida. Por u ´ltimo, utilizamos una variable suma para acumular los términos obtenidos hasta el momento. Ejercicios 

b

5

Escribe y ejecuta el programa anterior para el cálculo de 1/(1

− x) en x =0.6.

Escribe un programa que calcule 1/(1+ x2 ), siendo x un valor introducido por el usuario, mediante el siguiente desarrollo en serie de Taylor: 

6

1 = 1 + x2

∞

(−1) x i

i=0

2i

=1

2

−x

+ x4

−x

6

+ ...

.

La precisión también deberá ser introducida por el usuario. Fij´ emonos en que el signo de cada término de la serie cambia alternativamente entre + y . Calcula el desarrollo para x = 0.01.

−

5.6.

Bucles anidados

Los bucles pueden anidarse, es decir, podemos introducir un bucle dentro de otro, de forma que para cada iteración del bucle externo se ejecutan todas las iteraciones del bucle interno. Veamos un ejemplo. Supongamos que queremos escribir un programa que, tras leer un número entero N , diga qué números entre 1 y N son primos. Una implementación burda de la condición de primo implica el análisis de todos los posibles divisores desde 2 hasta N/2. As´ı pues, en este caso podemos utilizar dos bucles for: uno que recorra todos los enteros i desde 1 hasta N , y otro que (para cada i) recorra todos los posibles divisores desde 2 hasta i/2. N = input(’Introduce el valor de un numero entero: ’); for i = 1:N primo = 1; for div = 2:i/2 if rem(i,div) == 0 primo = 0; end end if primo == 1 disp(i); end end

Inicialmente, decidimos que todos los enteros son primos ( primo = 1), pero si encontramos un numero que divide al entero con un resto 0 ( if rem(i,div) == 0) entonces el n´ umero ya no es primo ( primo = 0). Fijémonos en que aunque encontremos un divisor para el número analizado, el bucle for interno sigue ejecutándose. En este caso resulta útil la sentencia break, que suspende la ejecución de un bucle cuando se cumple una cierta condición. De esta forma, una versión mejorada del programa anterior ser´ıa la siguiente: Ingenier´ıa Qu´ımica


5-10


N = input("Introduce el valor de un n´ umero entero: "); for i = 1:N primo = 1; for div = 2:i/2 if rem(i,div) == 0 primo = 0; break; % Terminamos el bucle for interno end end if primo == 1 disp(i); end end

No obstante, podemos evitar el uso de la orden break mediante un bucle condicional: N = input("Introduce el valor de un n´ umero entero: "); for i = 1:N div = 2; while (div*div <= i) & (rem(i,div) = 0) div++; end if div*div >i disp(i); end endfor

∼

En el siguiente ejemplo, se utilizan dos bucles for para crear una matriz A en la que el elemento (i, j) se calcula como la suma de los cuadrados de los números de fila y columna, o sea, A ij = i 2 + j 2 . Fijémonos en que se inicializa la matriz (creando, por ejemplo, una matriz de ceros) antes de poder utilizarla 3 .

 A = zeros(5) ;  for i = 1:5

for j = 1:5 A(i,j) = i∧ 2+ j∧ 2; end

end disp(A); 2 5 10 5 8 13 10 13 18 17 20 25 26 29 34



17 20 25 32 41

26 29 34 41 50

3

Realmente, en Octave no es necesaria esta inicializaci´ on, ya que se crean las matrices y se modifica su tama˜ no cuando es necesario.

Universitat Jaume I


5.7 Aplicación

5-11

Ejercicios 

7

Escribe y ejecuta el programa para el cálculo de n´ umeros primos con N = 100.

Dados los vectores x = [ 4 1 6 - 1 - 2 2 ] e y = [ 6 2 - 7 1 5 - 1 ] , escribe el código ıa matrices según las siguientes fórmulas: Octave que calcular´ 

8

a) aij = y j /xi b) bij = x i /(2 + xi + yj ) c) cij = 1/max(xi , yj ) b

9 Escribe un programa que transponga una matriz cuadrada (crea una aleatoria con on (’). rand(5). Comprueba que funciona comparando el resultado con el operador transposici´ 

5.7.

Aplicaci´ on

En la aplicaci´ on de esta unidad vamos a utilizar bucles para calcular la energ´ıa de repulsión nuclear entre los átomos de una molécula. Para calcular esta energ´ıa, utilizaremos la ecuación Z i Z j E ij = Rij

·

donde Z i es la carga nuclear del átomo i y Rij es la distancia que separa los átomos i y j. Por lo tanto, necesitamos pedir al usuario que introduzca un vector con las cargas de los N átomos, y una matriz con las distancias internucleares entre átomos. Esta matriz de distancias tiene una serie de caracter´ısticas especiales: por un lado, es simétrica (la distancia entre los átomos i y j es independiente del orden de las etiquetas) y presenta ceros en la diagonal (la distancia entre un átomo y s´ı mismo es nula). Z = input(’Introduce el vector de cargas nucleares: ’); R = input(’Introduce la matriz de distancias: ’);

A continuaci´ on, vamos a calcular cada uno de los elementos de la matriz E (i, j). Si tenemos en cuenta que los elementos de la diagonal van a ser todos infinito (recuerda: la diagonal de la matriz de distancias se comopone de ceros), podemos evitarnos el cálculo de la diagonal (y de paso, algunos posibles errores o avisos). Podemos inicializar la matriz E de la siguiente forma (piénsalo): N = length(Z); E = diag(zeros(1,N)*Inf);

Como tenemos que calcular la matriz E (i, j) para todos los valores posibles de i y j , debemos utilizar un bucle doble: for i = 1:N for j = 1:N E(i,j) = Z(i)*Z(j)/R(i,j); end end



5-12

Instrucciones repetitivas Pero si tenemos en cuenta que no necesitamos calcular la diagonal, y que la matriz es simétrica, basta con que calculemos uno de los triangulos de la matriz. Para ello, variamos el ´ındice j desde i + 1 hasta N , calculamos con la ecuación el elemento E (i, j), e igualamos E ( j, i) al mismo valor : for i = 1:N for j = i+1:N E(i,j) = Z(i)*Z(j)/R(i,j); E(j,i) = E(i,j); end end

En el segundo apartado se nos pide la energ´ıa total de repulsión. Para obtenerla, debemos sumar todas las energ´ıas de repulsi´ onentre átomos sin repeticiones, es decir, no podemos sumar la energ´ıa de repulsión E (i, j) y la energ´ıa E ( j, i), puesto que representan la misma interacción. Esto se consigue sumando sólo una de las energ´ıas. Además, debemos inicializar la variable que acumula la suma antes del bucle: Etotal = 0; for i = 1:N for j = i+1:N E(i,j) = Z(i)*Z(j)/R(i,j); E(j,i) = E(i,j); Etotal = Etotal + E(i,j); end end

Por u ´ ltimo, la b´ usqueda del máximo se consigue almacenando el valor mayor calculado en cada momento, y los ´ındices que lo han generado. Basta con a˜ nadir unas l´ıneas al código anteriormente escrito. A continuación, se muestra el código completo, con comentarios (f´ıjate en la inicializaci´ o n de Etotal a para la b´ usqueda del máximo). Despu´ es se muestra un ejemplo de uso para la molécula de agua.

−∞

Universitat Jaume I


5.7 Aplicación

5-13

%****************************************************** % Programa : repulsion.m %****************************************************** % Pedimos datos al usuario Z = input(’Introduce el vector de cargas nucleares: ’); R = input(’Introduce la matriz de distancias: ’); % Inicializamos E con infinitos en la diagonal % y ceros en el resto de elementos (esto no es necesario). N = length(Z); E = diag(zeros(1,N)*Inf); Etotal = 0; imin = -1; jmin = -1; Emax = -Inf; for i = 1:N for j = i+1:N E(i,j) = Z(i)*Z(j)/R(i,j); Etotal = Etotal + E(i,j); E(j,i) = E(i,j); % Buscamos m´ aximo if E(i,j) > Emax Emax = E(i,j); imin = i; jmin = j; end end end fprintf(’La energ´ ıa total de repulsi´ on es %g\n’, Etotal’) fprintf(’Los atomos que mas se repelen son %g y %g.\n’, imin,jmin)

 octave

-q repulsion.m Introduce el vector de cargas nucleares: [16 1 1] Introduce la matriz de distancias: [0.0000 0.9895 0.9895; 0.9895 0.0000 1.5163; 0.9895 1.5163 0.0000] La energ´ ı a total de repulsi´ on es 32.9991 Los atomos que mas se repelen son los atomos 1 y 2.

Ejercicios 

10

Escribe y ejecuta el programa anterior.



5-14


5.8.


Es conveniente que pienses y realices los ejercicios que han aparecido a lo largo de la unidad marcados con el s´ımbolo b antes de acudir a la sesión de prácticas correspondiente. Deberás iniciar la sesión realizando los ejercicios marcados con el s´ımbolo . A continuaci´ on, deberás hacer el mayor número de los ejercicios siguientes. Ejercicios Escribe un programa de nombre pi1.m que calcule el valor de π utilizando la siguiente serie matemática: pi2 8 1 = . 2 2 16 (2n 1) (2n + 1) n=1  11

∞

−



−

¿Cu´ antas iteraciones son necesarias para obtener π de forma que la precisión asociada al sumatorio sea 1e-12? 

12

Otra forma de calcular pi se basa en el siguiente m´ etodo:

√

Inicializaci´ on: a = 1, b = 1/ 2, t = 1/4 y x = 1. Repite las siguientes órdenes hasta que la diferencia entre a y b se encuentre por debajo de una cierta precisión: 1. y = a 2. a = (a + b)/2 3. b =

 (b ∗ y)

4. t = t 5.

2

− x(y − a) x = 2 ∗ x

Con los valores resultantes de a, b y t, se calcula la estimación de π como: (a + b)2 π = 4t Implementa este programa ( pi2.m) con Octave y calcula el n´ umero de iteraciones necesarias para obtener una precisión de 1e-12. Compara este resultado con el del ejercicio anterior. 

al es el valor de la variable ires tras la ejecución de este código? 13 ¿Cu´ ires = 0 ; for index1 = 10:-2:4 for index2 = 2:2:index1 if index2 == 6 break end ires = ires + index2; end end

Universitat Jaume I



5-15

Podemos representar el conjunto de enlaces de una mol´ ecula de N a´tomos mediante una matriz D de N N , en la que el elemento D ij es un 1 si los átomos i y j están enlazados y 0 en caso contrario. Veámoslo con un ejemplo. 

14

×

0 0 0 0 0

2

1

3 4 6

7

5

0 0

1

1

0

0 0 0 0

1 0 1 0 0 1 1 0

0 1

0

0 0 0

0 0 0 0 0

0 0

    

0 0 1 0 0 0 1 0 1 0 1 0 1 0

F´ıjate en que la matriz es simétrica (si un átomo está enlazado con otro, el otro lo está a su vez con el primero), y que se ha tomado el convenio de poner ceros en la diagonal, porque no hace falta indicar que un átomo está enlazado consigo mismo. Escribe un programa que pida al usuario una matriz de enlaces y le diga qué pares de átomos están enlazados. Para ello, deber´ as tener en cuenta lo siguiente: El programa deberá verificar que la matriz es simétrica. Recuerda que una matriz es simétrica si es igual a su transpuesta. ola vez , y no una por cada ´ S´ olo debe indicar cada par de átomos enlazados una s´ atomo del par.

Para facilitarte la escritura del programa, a continuación se muestra un ejemplo de su uso. % Ejemplo de uso: enlaces Introduce la matriz de enlaces: [0 0 1 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 1 Los atomos 1 y 3 estan enlazados. Los atomos 2 y 3 estan enlazados. Los atomos 3 y 4 estan enlazados. Los atomos 3 y 6 estan enlazados. Los atomos 5 y 6 estan enlazados. Los atomos 6 y 7 estan enlazados.



0 0 0 0 1 0]

Queremos calcular y dibujar la curva de valoraci´ o n de un ácido fuerte con una base fuerte. Para ello, vamos a escribir un Programa que pida al usuario los siguientes valores: 

15

V a : Volumen inicial de ácido. C a : Concentración del ácido. C b : Concentración de la base. V b : Un vector con distintos valores del volumen de la base a˜ nadida. Con estos valores, el Programa deberá calcular el pH para cada valor de V b , utilizando para ello las siguientes ecuaciones:

  ∆  − log V + V + 10  pH = 14 + log  −∆  V + V −7

10

a

b

10

a


si ∆

≥0

si ∆ < 0

b


5-16

Instrucciones repetitivas donde ∆ se define como ∆ = C a V a

− C bV b

.

Adem´ as, el Programa deberá representar el pH en función del volumen de base a˜ nadido V B . on para calcular el logaritmo decimal en Octave es log10. Ayuda: La funci´ % Ejemplo pH Introduce Introduce Introduce Introduce



de uso: el la la un

volumen inicial de a ćido (ml): 25 concentracion inicial de a ćido (M): 0.1 concentracion inicial de base (M): 0.1 vector de vol´ umenes de base (ml): [0:50]

Curva de valoracion acido fuerte-base fuerte 14

12

10

8 H p

6

4

2

0 0

5

10

15

20

25

30

35

40

45

50

Vb(ml)

Se tienen N disoluciones numeradas de 1 a N y se mide el pH y la temperatura de cada disolución. 

16

pH 1 T1

pH

2 T2

pH 3 T 3

pH 4 T4

a) Queremos escribir un programa que pida al usuario el número de disoluciones N , y los valores de temperatura y pH de las N disoluciones, y muestre en pantalla: La temperatura media de las disoluciones claramente ácidas (pH < 6.5). La temperatura media de las disoluciones claramente básicas (pH > 7.5). La temperatura media de las disoluciones neutras (6.5

≤ pH ≤ 7.5).

A continuación tienes un ejemplo de cómo deber´ıa usarse el programa: % Ejemplo de uso: pH ......... La temperatura media de las disoluciones a ćidas es: 276.33 K La temperatura media de las disoluciones neutras es: 238.00 K La temperatura media de las disoluciones b´ asicas es: 280.50 K



Puedes pedirle los datos al usuario como creas conveniente.

Universitat Jaume I



5-17

b) A˜ nade al programa anterior las l´ıneas que consideres oportunas para que eval´ ue cuál es la disoluci´ on de mayor temperatura, y el pH de ésta. c) Puede que al escribir el programa anterior no hayas tenido en cuenta que pudiera no haber disoluciones en uno de los grupos de pH (por ejemplo, podr´ıa no haber disoluciones ácidas). En ese caso, al ejecutar el programa se generar´ıa un error al calcular la media (por una divisi´ on por cero). Si no lo has hecho ya, incluye las modificaciones necesarias para que, en caso de que no haya disoluciones de un grupo, el programa se lo indique al usuario. Por ejemplo, % Ejemplo de uso: pH ......... La temperatura media de las disoluciones acidas es: 276.33 K No hay disoluciones neutras. La temperatura media de las disoluciones basicas es: 280.50 K





Tema 6

Funciones Guillermo Peris Ripoll´ es Objetivos Cuando finalice este tema, el alumno deberá ser capaz de: Definir y utilizar funciones en Octave.

Aplicaci´ on Cuando finalice este tema, el alumno deberá ser capaz de resolver problemas como el siguiente, cuya resolución se indica a lo largo del propio tema. Trayectoria de un proyectil

La siguiente figura muestra la trayectoria de un proyectil lanzado a una velocidad v0 y un ángulo θ sobre un plano horizontal.

y(t)

v0

Altura θ

x(t) Alcance

Suponiendo que pueden despreciarse todos los efectos de resistencia del aire, las ecuaciones del movimiento en los dos ejes vienen dadas por las ecuaciones siguientes: x(t) = v o cos(θ)t 1 y(t) = vo sin(θ)t − gt 2 2 donde g = 9.81m/s2 . Igualando y (t) = 0 podemos obtener el tiempo máximo de vuelo, que es el momento en el que el proyectil toca el suelo. Tiempo de vuelo =


2v0 sin(θ) g


6-2

Funciones

A partir de estas ecuaciones, escribe: 1. Una funci´ on que calcule un vector de posiciones (x), en funci´ on del tiempo para una velocidad y ángulo inicial determinado. 2. Una funci´ on que calcule un vector de alturas (y), en funci´ on del tiempo para una velocidad y ángulo inicial determinado. 3. Una funci´ on que calcule el tiempo de vuelo para una velocidad y ángulo inicial determinado. 4. Una funci´ on que represente la trayectoria del proyectil a partir de los vectores de posici´ on y altura. 5. Un programa que pida al usuario la velocidad y ángulo inicial del proyectil y represente su trayectoria, utilizando todas las funciones anteriores.

Contenidos

6.1. Definici´ on y uso de funciones . . . . . . . . . . . . . . . . . 6-3 6.2. Dise˜ no de programas con funciones . . . . . . . . . . . . . .

6-5

6.3. Aplicaci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6-7

6.4. Ejercicios pr´ acticos

Universitat Jaume I

. . . . . . . . . . . . . . . . . . . . . . . 6-14


6.1 Definición y uso de funciones

6-3

En la unidad de introducción a Octave vimos cómo era posible ejecutar directamente órdenes en la ventana de Octave, pero que en problemas más complejos resultaba más conveniente utilizar programas. Estos ficheros son útiles cuando deseamos repetir un conjunto de órdenes, cambiando (o no) algunos de los valores de las variables. Sin embargo, no son más que una mecanización de la introducción directa de órdenes. Asimismo, tambi´ en hemos utilizado funciones Octave, como por ejemplo sqrt, que recib´ıan (o no) una serie de parámetros de entrada y devolv´ıan (o no) un conjunto de parámetros de salida. En esta unidad, aprenderemos a definir nuestras propias funciones Octave, lo cual nos permitirá escribir programas de una forma más estructurada y reutilizar código.

6.1.

Definici´ on y uso de funciones

Antes de detallar las reglas que deben seguir las funciones vamos a considerar un ejemplo sencillo: la función estad(x), que se define a continuación, calcula la media y desviación t´ıpica de los datos contenidos en un vector: function [media, std] = estad(x) % ESTAD(x) Estadistica simple % Calcula la media y desviacion tipica de un vector x. n = length(x); media = sum(x)/n; v = x - media; std = sqrt(v*v’/(n-1)) ; endfunction

La función estad(x) recibe un argumento de entrada que es un vector x y devuelve un argumento de salida que tambi´ en es un vector [media, std] que contiene los valores de la media y desviación t´ıpica. La sintaxis general para la definici´ o n de funciones es la siguiente: function [out1 out2 ... outM] = nombre_funcion (in1, in2, ...inN)

A esta l´ınea se la conoce como declaración de la función. Fijémonos en que pueden existir varios argumentos de salida como elementos de un vector, por lo que deberán aparecer entre corchetes. En el caso de un solo argumento de salida los corchetes no son necesarios, y tambi´ en es posible que no haya ningún argumento de salida. Del mismo modo, puede pasarse a la función una lista de argumentos de entrada separados por comas, e incluso no pasar ningún argumento. Cada uno de estos argumentos puede ser un escalar, vector o matriz. El fichero que contenga a la función debe tener su mismo nombre, por lo que la funci´ on estad debe guardarse en un fichero de nombre estad.m. La primera l´ınea del fichero debe ser necesariamente la declaración de la función. Tras la declaración de la función, aparecen un conjunto de comentarios que documentan el uso de la función. Estas l´ıneas son las que aparecerán cuando se pida ayuda sobre la función a Octave:



6-4

Funciones

 help

estad estad is the user-defined function from the file /home/user/estad.m ESTAD(x) Estadistica simple Calcula la media y desviacion t´ ıpica de un vector x.

Por u ´ ltimo, el resto de la función se denomina el cuerpo de la funci´ emonos on . Fij´ en que, en algún momento, deben definirse los valores de los argumentos de salida, en el ejemplo las variables media y std, utilizando los argumentos de entrada, en este caso el vector x. La orden endfunction indica que la ejecución de la función ha terminado 1 .

Llamada a funciones Una vez definida la función, para utilizarla debemos ejecutar ordenes ´ de llamada a la funci´ on , de la misma forma que lo hac´ıamos con las funciones predefinidas por ı, a continuación se muestra un ejemplo de llamada a la función estad: Octave. As´  vec

= rand(1,5);  [a b] = estad(vec) a = 0.24375 b = 0.20014 F´ ıjate en que los nombres de las variables que se le pasan a la funci´ on no tienen por qu´ e coincidir con la definici´ on de la funci´ on. As´ı, mientras que en la definici´ o n de estad el argumento de entrada era x, aqu´ ı hemos utilizado la variable vec para realizar la llamada, y lo mismo ocurre con los argumentos de salida.

Por u ´ ltimo, recuerda que para que Octave localice una de tus funciones, debes ejecutar el programa en el directorio en el que se encuentre la función (podemos evitar ésto, cambiando el valor de una variable de entorno de Octave. Si te interesa saber cómo se hace, preg´ untale a tu profesor). Ejercicios Escribe la función estad, gu´ ardala con el nombre estad.m en tu directorio de trabajo, y comprueba su funcionamiento con vectores de valores aleatorios (utiliza la función rand(1,n), cambiando n para vectores de distinta longitud).  1

b



2

Comprueba que obtienes las l´ıneas de ayuda con la orden help estad.



3

El a´rea de un triángulo de lados a, b y c viene dada por la ecuación area =



s(s − a)(s − b)(s − c) ,

donde s = (a + b + c)/2. Escribe una función que acepte a, b y c como argumentos de entrada y devuelva el área del triángulo como salida.

1

Mucho ojo para los usuarios de Matlab!!! Debeis usar return en lugar de endfunction .

Universitat Jaume I


6.2 Diseño de programas con funciones

6.2.

6-5

Dise˜ no de programas con funciones 2

Ya hemos aprendido a escribir funciones, pero no está claro qué ventajas nos reporta trabajar con ellas. A grandes rasgos, el uso de funciones facilita la lectura del código y su propia programación. Veámoslo con un ejemplo. Supongamos que en un programa deseamos leer dos números enteros y asegurarnos de que sean positivos. Podemos proceder repitiendo el bucle correspondiente dos veces: % Leemos a do a = input(’Introduce if a < 0 disp(’Error: el end until a > 0 % Leemos b do b = input(’Introduce if b < 0 disp(’Error: el end until b > 0 end

un numero positivo: ’); numero debe ser positivo’);

un numero positivo: ’); numero debe ser positivo’);

o dise˜ nar una función que lea un número asegurando que sea positivo, function numero = lee_numero_positivo % lee_numero_positivo: Pide un numero al usuario % y comprueba que sea positivo do numero = input(’Introduce un numero positivo’); if numero < 0 disp(’Error: el numero debe ser positivo’); end until numero > 0 endfunction

y realizar dos llamadas a la misma, a = lee_numero_positivo; b = lee_numero_positivo;

Un programa que utiliza funciones es, por regla general, más legible que uno que inserta los procedimientos de cálculo directamente donde se utilizan. . . siempre que se escojan nombres de función que describan bien qué hacen éstas. F´ıjémonos en que este u ´ ltimo programa es más fácil de leer que el original. 2

Esta secci´ on y la siguiente han sido adaptadas de los apuntes para Metodolog´ıa y Tecnolog´ıa de la Programaci´ on de primer curso de Ingenier´ıa Inform´ atica, escritos por Andr´ es Marzal.



6-6

Funciones Las funciones son un elemento fundamental de los programas. Ahora ya sabemos c´ omo construir funciones, pero no cu´ ando conviene construirlas. Lo cierto es que no podemos decirlo: no es una ciencia exacta, sino una habilidad que se va adquiriendo con la práctica. De todos modos, s´ı podemos dar algunos consejos. 1. Por una parte, todos los fragmentos de programa que se utilizan en m´ as de una ocasi´ on son buenos candidatos a definirse como funciones , pues de ese modo evitamos tener que copiarlos en varios lugares. Evitar esas copias no sólo resulta má s cómodo: tambi´ en reduce considerablemente la probabilidad de que cometamos errores, pues acabamos escribiendo menos texto. Además, si acabamos cometiendo errores y hemos de corregirlos o si hemos de modificar el programa para ampliar su funcionalidad, siempre será mejor que el mismo texto no aparezca en varios lugares, sino una sola vez en una función. 2. No conviene que las funciones que definamos sean muy largas. En general, una funci´ on deber´ıa ocupar menos de 30 o 40 l´ıneas (aunque siempre hay excepciones). Una funci´ on no s´ olo deber´ıa ser breve, adem´ as deber´ıa hacer una ´ unica cosa. . . y hacerla bien. Deber´ıamos ser capaces de describir con una sola frase lo que hace cada una de nuestras funciones. Si una función hace tantas cosas que explicarlas todas cuesta mucho, probablemente har´ıamos bien en dividir la funci´ on en otras funciones más peque˜ nas y simples. Recordemos que siempre podemos llamar a una función desde otra. El proceso de identificar acciones complejas y dividirlas en acciones más sencillas se conoce como estrategia de dise˜ no descendente (en inglés, top-down ). La forma de proceder es ésta: analizar primero qu´ e debe hacer el programa y hacer un diagrama que represente las diferentes acciones que debe efectuar, pero sin entrar en los detalles de cada una de ellas; descomponer el programa en una posible función por cada una de esas acciones; analizar entonces cada una de esas acciones y ver si aún son demasiado complejas; si es as´ı, aplicar el mismo método hasta que obtengamos funciones más peque˜ nas y simples. Una estrategia de dise˜ no alternativa recibe el calificativo de ascendente (en inglés, bottom-up ) y consiste en lo contrario: detectar algunas de las acciones más simples que necesitaremos en el programa y escribir pequeñas funciones que las implementen; combinar estas acciones en otras más sofisticadas y crear nuevas funciones para ellas; seguir hasta llegar a una o unas pocas funciones que resuelvan el problema planteado. Cuando se empieza a programar resulta dif´ıcil anticiparse y detectar a simple vista qué peque˜ nas funciones irán haciendo falta y cómo combinarlas apropiadamente. Ser´ a más efectivo empezar siguiendo la metodolog´ıa descendente: dividir cada problema en subproblemas má s y más sencillos que, al final, se combinan para dar solució n al problema original. Cuando tengamos mucha más experiencia, descubriremos que al programar seguimos una estrategia h´ıbrida, ascendente y descendente a la vez. Todo llega. Paciencia.

Universitat Jaume I


6.3 Aplicación

6.3.

6-7

Aplicaci´ on

A continuaci´ on se muestran las funciones que se piden en la aplicación, as´ı como el programa que las utiliza. A estas alturas del curso, no deber´ıas de necesitar ninguna explicaci´ on para comprender este código. function x = posicion(v0, theta, t) % Funcion: x = posicion(v0, theta, t) % Calcula la posici´ on de un objeto % lanzado con velocidad inicial v0 % y angulo theta, para distintos % valores de tiempo t. x = v0*cos(theta*pi/180)*t; endfunction

function y = altura(v0, theta, t) % Funcion: y = altura(v0, theta, t) % Calcula la altura de un objeto % lanzado con velocidad inicial v0 % y angulo theta, para distintos % valores de tiempo t. g =9.81 ; % Constante gravedad en m/s2 y = v0*sin(theta*pi/180)*t - 0.5*g*t.∧2; endfunction

function t = tvuelo(v0, theta) % Funcion: t = tvuelo(v0, theta) % Calcula el tiempo de vuelo de un objeto % lanzado con velocidad inicial v0 % y angulo theta. g = 9.81 ; %Constante gravedad en m/s2. t = 2*v0*sin(theta*pi/180)/g; endfunction

function trayectoria(x,y) % Funcion: trayectoria(x,y) % Dibuja la trayectoria de un objeto % a partir de vectores de posici´ on x e y. clg plot(x,y,’;;’), title(’Trayectoria’) xlabel(’x(m)’), ylabel(’y(m)’) replot endfunction



6-8

Funciones

%======================================================= % Programa proyectil.m: Pide al usuario los valores % iniciales de velocidad y angulo ´ % de un proyectil, y representa % su trayectoria. %======================================================= % Pedimos los valores al usuario v0 = input(’Introduce el valor de velocidad inicial en m/s: ’); theta = input(’Introduce el valor del angulo ´ inicial en grad: ’); % Calculamos el vector de tiempos t = linspace(0,tvuelo(v0,theta),50); % Hacemos la representacion trayectoria(posicion(v0,theta,t), altura(v0,theta,t) )

Universitat Jaume I


6.3 Aplicación

6-9

Ap´ endice: Recursi´ on Desde una función podemos llamar a otras funciones. Ya lo hemos hecho en los ejemplos que hemos estudiado, pero ¿qué ocurrir´ıa si una función llamara a otra y ésta, a su vez, llamara a la primera? O de modo más inmediato, ¿qué pasar´ıa si una función se llamara a s´ı misma? La capacidad de que una función se llame a s´ı misma, directa o indirectamente, se denomina recursi´ on es un potente concepto con el que se pueden expresar on . La recursi´ procedimientos de cálculo muy elegantemente. No obstante, al principio cuesta un poco entender las funciones recursivas. . . y un poco más dise˜ nar nuestras propias funciones recursivas. La recursión es un concepto dif´ıcil cuando estamos aprendiendo a programar. Por ello, no debes asustarte si este material se te resiste más que el resto.

C´ alculo recursivo del factorial Empezaremos por presentar y estudiar una función recursiva: el cálculo recursivo del factorial de un n´ umero. Partiremos de la siguiente fórmula matem´ atica para el cálculo del factorial: n · (n − 1)!, si n > 1. n! = 1, si n = 1.



Es una definición de factorial un tanto curiosa, pues ¡se define en términos de s´ı misma! El primero de sus dos casos dice que para conocer el factorial de n hay que conocer el factorial de n − 1 y multiplicarlo por n. Entonces, ¿cómo calculamos el factorial de n − 1? En principio, conociendo antes el valor del factorial de n − 2 y multiplicando ese valor por n − 1. ¿Y el de n − 2? Pues del mismo modo. . . y as´ı hasta que acabemos por preguntarnos cuánto vale el factorial de 1. En ese momento no necesitaremos hacer más cálculos: el segundo caso de la fórmula nos dice que el factorial de 1 vale 1. Vamos a plasmar este mismo procedimiento en una función Octave: function f = factorial(n) % FACTORIAL(n): Calcula el factorial de n de forma recursiva. if n > 1 f = n * factorial(n-1) ; else f = 1 ; end endfunction

Compara la f´ ormula matem´ atica y la función Octave. No son tan diferentes. Ocas allá tave nos fuerza a decir lo mismo de otro modo, es decir, con otra sintaxis . M´ de las diferencias de forma, ambas definiciones son idénticas en el fondo. Para entender la recursión, nada mejor que verla en funcionamiento. La Figura 6.1 nos muestra paso a paso qu´ e ocurre si solicitamos el cálculo del factorial de 5. Con el anidamiento de cada uno de los pasos pretendemos ilustrar que el cálculo de cada uno de los factoriales tiene lugar mientras el anterior aún está pendiente de completarse. En el nivel más interno, factorial(5) está pendiente de que acabe factorial(4), que a su vez est´ a pendiente de que acabe factorial(3), que a su vez está pendiente de que acabe factorial(2), que a su vez está pendiente de que acabe factorial(1). Cuando factorial(1) acaba, pasa el valor 1 a factorial(2), que a su vez pasa el Ingenier´ıa Qu´ımica


6-10

Funciones Empezamos invocando factorial(5). Se ejecuta, pues, la l´ınea 3 y como n es mayor que uno, pasamos a ejecutar la l´ınea 4. Hemos de calcular el producto de n por algo cuyo valor es aún desconocido: factorial(4). El resultado de ese producto se almacenará en la variable resultado, pero antes hay que calcularlo, as´ı que hemos de invocar a factorial(4). Invocamos ahora factorial(4). Se ejecuta la l´ınea 3 y como n, que ahora vale 4, es mayor que uno, pasamos a ejecutar la l´ınea 4. Hemos de calcular el producto de n por algo cuyo valor es a´ un desconocido: factorial(3). El resultado de ese producto se almacenará en la variable resultado, pero antes hay que calcularlo, as´ı que hemos de invocar a factorial(3). Invocamos ahora factorial(3). Se ejecuta la l´ınea 3 y como n, que ahora vale 3, es mayor que uno, pasamos a ejecutar la l´ınea 4. Hemos de calcular el producto de n por algo cuyo valor es a´ un desconocido: factorial(2). El resultado de ese producto se almacenar´ a en la variable resultado, pero antes hay que calcularlo, as´ı que hemos de invocar a factorial(2). Invocamos ahora factorial(2). Se ejecuta la l´ınea 3 y como n , que ahora vale 2, es mayor que uno, pasamos a ejecutar la l´ınea 4. Hemos de calcular el producto de n por algo cuyo valor es a´ un desconocido: factorial(1). El resultado de ese producto se almacenará en la variable resultado, pero antes hay que calcularlo, as´ı que hemos de invocar a factorial(1). Invocamos ahora factorial(1). Se ejecuta la l´ınea 3 y como n no es mayor que 1, pasamos a la l´ınea 6. En ella se dice que resultado vale 1, y en la l´ınea 8 ése será el valor devuelto como resultado de calcular factorial(1). Ahora que sabemos que el valor de factorial(1) es 1, lo multiplicamos por 2 y almacenamos el valor resultante, 2, en resultado. Al ejecutar la l´ınea 8, ´ ese será el valor devuelto. Ahora que sabemos que el valor de factorial(2) es 2, lo multiplicamos por 3 y almacenamos el valor resultante, 6, en resultado. Al ejecutar la l´ınea 8, ése será el valor devuelto. Ahora que sabemos que el valor de factorial(3) es 6, lo multiplicamos por 4 y almacenamos el valor resultante, 24, en resultado. Al ejecutar la l´ınea 8, ése será el valor devuelto. Ahora que sabemos que el valor de factorial(4) es 24, lo multiplicamos por 5 y almacenamos el valor resultante, 120, en resultado. Al ejecutar la l´ınea 8, ése será el valor devuelto. alculo recursivo de factorial(5). Figura 6.1: Traza del c´

valor 2 a factorial(3), que a su vez pasa el valor 6 a factorial(4), que a su vez pasa el valor 24 a factorial(5), que a su vez devuelve el valor 120. De acuerdo, la figura 6.1 describe con mucho detalle lo que ocurre, pero es dif´ıcil de seguir y entender. Veamos si la figura 6.2 nos es de más ayuda. En esa figura también se describe paso a paso lo que ocurre al calcular el factorial de 5, sólo que con la ayuda de unos mu˜ necos. ´ no sabe En el paso 1, le encargamos a Amadeo que calcule el factorial de 5. El calcular el factorial de 5, a menos que alguien le diga el valor del factorial de 4. En el paso 2, Amadeo llama a un hermano clónico suyo, Benito, y le pide que calcule el factorial de 4. Mientras Benito intenta resolver el problema, Amadeo se echa a dormir (paso 3). Universitat Jaume I


6.3 Aplicación

6-11

1)

10)

5!

2)

5!

4!

3!

5!

4!

3!

5!

4!

3!

5!

4!

5!

4!

2!

1! = 1

11)

5!

→

4!

3)

2!

←

1! = 1

12)

5!

4!

4)

2! = 2 1 ·

13)

5!

4!

→

3!

5)

3!

←

2! = 2

14)

5!

4!

5!

4!

3!

6)

3! = 3 2 ·

15)

3!

→

2!

5!

7)

4!

←

3! = 6

16)

5!

4!

3!

2!

5!

8)

4! = 4 6 ·

17)

5!

4!

3!

5!

4!

3!

2!

→

1!

5!

9)

←

4! = 24

18)

2!

1!

5! = 5 24 ·

omic explicativo del cálculo recursivo del factorial de 5. Figura 6.2: C´

Benito tampoco sabe resolver directamente factoriales tan complicados, as´ı que llama a su clon Ceferino en el paso 4 y le pide que calcule el valor del factorial de 3. Mientras, Benito se echa a dormir (paso 5). La cosa sigue igual un ratillo: Ceferino llama al clon David y David a Eduardo. As´ı llegamos al paso 9 en el que Amadeo, Benito, Ceferino y David est´ an durmiendo y Eduardo se pregunta cuánto valdr´ a el factorial de 1. En el paso 10 vemos que Eduardo cae en la cuenta de que el factorial de 1 es muy fácil de calcular: vale 1. En el paso 11 Eduardo despierta a David y le comunica lo que ha averiguado: el factorial de 1! vale 1. En el paso 12 Eduardo nos ha abandonado: él ya cumpli´ o con su deber. Ahora es David el que resuelve el problema que le hab´ıan encargado: 2! se puede calcular multiplicando 2 por lo que valga 1!, y Eduardo le dijo que 1! vale 1. Ingenier´ıa Qu´ımica


6-12

Funciones En el paso 13 David despierta a Ceferino para comunicarle que 2! vale 2. En el paso 14 Ceferino averigua que 3! vale 6, pues resulta de multiplicar 3 por el valor que David le ha comunicado. Y as´ı sucesivamente hasta llegar al paso 17, momento en el que Benito despierta a Amadeo y le dice que 4! vale 24. En el paso 18 sólo queda Amadeo y descubre que 5! vale 120, pues es el resultado de multiplicar por 5 el valor de 4!, que según Benito es 24. Ejercicios Escribe la funci´ on factorial en Octave y util´ızala para calcular los factoriales de varios n´ umeros enteros. b



4



5

También podemos formular recursivamente la suma de los n primeros n´ umero naturales:

   n

i =

i=1

1, n +

n

−1

i=1

si n = 1; i, si n > 1.

Dise˜ na una función Octave recursiva que calcule el sumatorio de los n primeros n´ umeros naturales. Comprueba que funciona correctamente comparando los resultados con los obtenidos con la orden sum(1:n).

Hemos propuesto una solución recursiva para el cálculo del factorial, pero también podemos escribir una función que realice este cálculo iterativamente , utilizando un bucle for: function f = factorial(n) % factorial: Calcula el factorial de un entero n. f = 1; for i = 2:n f = f * i; end endfunction

Pues bien, para toda función recursiva podemos encontrar otra que haga el mismo cálculo de modo iterativo. Ocurre que no siempre es fácil hacer esa conversión o que, en ocasiones, la versión recursiva es más elegante y legible que la iterativa (al menos se parece más a la definición matem´ atica). Por otra parte, las versiones iterativas suelen ser más eficientes que las recursivas, pues cada llamada a una función supone pagar una peque˜ na penalizaci´ on en tiempo de cálculo y espacio de memoria. Ejercicios b

Los n´ umeros de Fibonacci son una secuencia de números enteros bastante curiosa, ya que la naturaleza está llena de situaciones en las que éstos aparecen. A continuación, se indican los primeros n´ umeros de esta serie: 

6

1 1 2 3

5 8 13 21 34 55 89

Los dos primeros n´ umeros de la secuencia valen 1 y cada número a partir del tercero se obtiene sumando los dos anteriores. Podemos expresar esta definición matemáticamente as´ı: F n =

Universitat Jaume I



F n−1 + F n−2 , si n > 2. 1, si n = 1 o n = 2. Guillermo Peris Ripoll´ es

6.3 Aplicación

6-13

a) Escribe una funci´ on que acepte como argumento un entero n y devuelva F n . b) Escribe una funci´ on que lea una lista de enteros y calcule recursivamente el número de Fibonacci F n de cada uno de ellos. Haz uso de la función definida anteriormente. c) Calcula los 20 primeros n´ umeros de Fibonacci, usando la función del apartado anterior. d) Para los primeros 20 n´ umeros de Fibonacci, calcula la relación F n /F n−1 y represéntala en √ 1+ 5 funci´ on de n. Esta relación tiende a 2 . ¿Qué muestran tus resultados? e) Escribe una funci´ on que calcule los primeros 20 números de Fibonacci sin recursión, utilizando bucles.



6-14

Funciones

6.4.


Es conveniente que pienses y realices los ejercicios que han aparecido a lo largo de la unidad marcados con el s´ımbolo b antes de acudir a la sesión de prácticas correspondiente. Deberás iniciar la sesión realizando los ejercicios marcados con el s´ımbolo . A continuaci´ on, deberás hacer el mayor número de los ejercicios siguientes. Ejercicios on Octave que calcule la presión de un gas utilizando la ley de los 7 1) Escribe una funci´ gases ideales (LGI), RT P = ˆ V 

y otra función que calcule la presión con la ecuación de Soave-Redlich-Kwong (SRK), P =

RT ˆ −b V

−

αa , ˆ (V ˆ + b) V

donde ˆ V a

= V /n = 0.42747R2 T c2 /P c

b = 0.08664RT c /P c atm · l R = 0.082 mol · K La constante α se calcula a partir del factor ac´ entrico de Pitzer ω utilizando dos ecuaciones adicionales: m = 0.48508 + 1.55171ω − 0.1561ω2 α =

    1 + m 1 −

2

T /T c

Al escribir las funciones debes documentarlas convenientemente, de forma que despu´ es cualquier usuario pueda obtener ayuda sobre ellas con la orden help. Adem´ as, puedes definir las subfunciones que consideres necesarias para facilitar la programación. 2) Utiliza estas funciones para hacer una gráfica de la presión en función de la temperatura para un cilindro de gas de 820 l. que contiene 100 moles de CO2 . El rango de temperaturas debe variar entre 20 y 400o C . La gráfica debe incluir las curvas predecidas por los modelos LGI y SRK. Para el CO2 , los valores de las constantes requeridos para la ecuación SRK son: ω = 0.225 T c = 304.2K P c

Universitat Jaume I

= 72.9 atm



6-15

Un flu´ıdo compuesto de moléculas con momento dipolar µ, a una temperatura T , y sometido a un campo eléctrico E , presenta un momento dipolar medio (que se puede medir experimentalmente) dado por la expresión 

8

µmedio = µΛ, on de Langevin , cuya ecuaci´ donde Λ es la funci´ on es:

Λ = coth(x) − siendo x =

1 , x

µE , y coth(x) la cotangente hiperbólica de x, kT ex + e−x coth(x) = x . e − e−x

Por u ´ ltimo, k = 1.38 · 10−23 J/K es la constante de Boltzmann. ´ n coth que reciba como argumento de entrada un vector de valores 1. Escribe una Funcio y devuelva como argumento de salida un vector con las cotangentes hiperbólicas del vector de entrada. on que calcula e x es exp(x). Ayuda: Recuerda que la funci´

´ n Langevin que reciba como argumento de entrada un vector de 2. Escribe una Funcio valores x y devuelva como argumento de salida un vector con los valores de la función de Langevin Λ del vector de entrada.

3. Escribe un Programa que pida al usuario : el momento dipolar molecular µ, el valor del campo eléctrico externo E , y el valor de la temperatura T en Kelvin. y calcule el momento dipolar medio µ medio .

on lineal de masas puntuales respecto a un 9 El momento de inercia I C de una distribuci´ eje de rotación (ver figura) viene dado por la siguiente ecuación: 

IC

I0

CM

m

m 1

m 2

I C = m

3

m 4

m 5

d

6



x2i mi .

i

donde x i es la distancia de la part´ıcula i al eje de rotación, y m i su masa. Si conocemos el momento de inercia respecto al eje que pasa por el centro de masas, basta con que apliquemos el teorema de Steiner, según el cual I C = I 0 + M · d2



donde M = i mi es la masa total del sistema, y d la distancia del eje C al centro de masas. En este ejercicio se pretende que escribas un conjunto de funciones y un programa que permita el c´ alculo del momento de inercia de un sistema de part´ıculas lineales. Para ello, sigue ordenadamente los siguientes pasos: Ingenier´ıa Qu´ımica


Programación en Ingeniería Química Con OCTAVE. Guillermo Peris Ripollés, Universitat Jaume I.

Recommend Documents