Manual de introducción a GNU-Octave¶

Arranque de la consola de Octave¶

Partimos de un ordenador con Octave instalado. Abriremos un terminal para poder ejecutar instrucciones del sistema operativo. En Linux simplemente hay que teclear ‘octave‘, en Windows el programa que tenemos que ejecutar para abrir la consola de Octave es un programa llamado octave.exe y que está en el directorio ‘bin’ de la instalación de Octave. Si todo ha ido bien, nos aparecerá la consola de Octave con el prompt indicando algo parecido a octave:1>.

Para salir de la consola de Octave hay que teclear ‘exit‘ o ‘quit‘

Cálculos elementales¶

Octave se puede utilizar como una calculadora. Se dispone de las operaciones aritméticas suma, resta, multiplicación, división y exponenciación. Pruebe a teclear las siguientes operaciones en la consola de Octave:

1+1
2-1
2*3.5
3.0/2.0
2^5

Tras cada operación hay que pulsar la tecla intro. La consola nos mostrará la palabra ans, el signo igual y el resultado de la operación. La palabra ans quiere decir answer, respuesta.

Octave dispone de funciones para calcular raices cuadradas, logaritmos naturales, logaritmos decimales y funciones trigonométricas. Pruebe las siguientes operaciones en la consola de Octave:

sqrt(25)
log(10)
log10(10)
sin(90*pi/180)
cos(0)
tan(45*pi/180)
asin(1)*180/pi
acos(0)
atan(-1)*180/pi

Observamos un par de cosas en las expresiones anteriores:

• La utilización de la constante predefinida pi. Octave tiene varias constantes con valores predefinidos, una de ellas es el número pi=3.14....
• Las funciones trigonométricas trabajan con ángulos en radianes.

Podemos calcular exponenciales del número , por ejemplo, la ecuación de Euler:

la podríamos comprobar tecleando algo así:

exp(pi*i)+1

Vemos que Octave entiende la variable i como el número complejo . Podemos ver algunas de las variables que tiene predefinidas:

e
i
j
pi
ans

Cada vez que hacemos un cálculo, el resultado se guarda en una variable llamada ans que podemos utilizar en el siguiente cálculo.

Asignación de variables¶

Podemos almacenar valores en memoria mediante la asignación de dichos valores a nombres válidos de variables. El símbolo utilizado para la asignación es el símbolo igual ‘=’.

En el ejemplo siguiente se asigna valor a las variables x e y. Estos valores quedan almacenados en la memoria de Octave y se pueden utilizar en operaciones posteriores utilizando en las expresiones el nombre de variable elegido. Se puede reasignar el valor de una variable en cualquier momento. El valor en memoria será el último asignado a la variable.

x = 2.5
y = 3
x+y
% El resultado será: ans = 5.5
x = 4
x+y
% El resultado será: ans = 7

Vemos que el esquema de la instrucción de asignación en *Octave es situar el nombre de variable, a continuación el símbolo igual y a la derecha del símbolo igual el valor que queremos asignar a la variable.

A la derecha del símbolo igual puede aparecer un valor numérico tecleado explicitamente o cualquier expresión válida de Octave. En el siguiente ejemplo se asignan distintos valores a variables:

x = pi/2;
y = sin(x)
% El resultado será: y = 1

Nombres de variables¶

Los nombres de variables pueden contener letras, números y caracteres underscore (guión bajo), pero el primer carácter tiene que ser letra.

Ejemplos de nombres válidos de variables:

x
x2
XX_2
ultimoValor
n_factorial

Ejemplos de nombres no válidos de variables:

6x
end
n!

El número máximo de caracteres que puede tener el nombre de una variable se puede consultar con la función namelengthmax():

namelengthmax()
% El resultado en mi consola: ans = 63

Variables predefinidas en Octave¶

Vectores¶

Octave trabaja con vectores de datos, tanto con vectores fila como con vectores columna.

Para definir un vector fila por teclado solo hay que teclear los elementos separados por un espacio o por una coma y encerrados entre corchetes.

Analogamente, un vector columna se define tecleando la lista de elementos separados por punto y coma ‘;‘ o por la pulsación de la tecla intro:

Rangos¶

En Octave podemos definir un tipo de datos especial, llamado Rango, consistente en una colección ordenada (una serie) de números equiespaciados.

Para definir un rango debemos indicar a Octave el valor inicial del rango, el incremento entre valores sucesivos y el máximo valor que pueden alcanzar los valores de la serie. La construcción para definir un rango coloca los tres valores separados por ‘:‘ :

% rango = inicio : incremento : máximo_valor
rango = 1:2:10
% El resultados será: 1 3 5 7 9

Cuando se omite el incremento, Octave asume que el incremento es la unidad:

% rango = inicio : máximo_valor
rango = 1:4
% El resultado será: 1 2 3 4

Se pueden establecer rangos de valores decrecientes. Para ello el valor de inicio deberá ser mayor que el valor final y el incremento un valor negativo:

rango = 4:-1:0
% El resultado será: 4 3 2 1 0

Para especificar los valores inicio, incremento o máximo de un rango se pueden utilizar funciones, o cualquier otra expresión de Octave que devuelva un resultado válido, por ejemplo:

rango = 0:intmax('uint8')
% El resultado serán los primeros 255 números naturales: rango = 0 1 2 ... 255
rango2 = 0:2:2^3
% El resultado será: rango2 = 0 2 4 6 8

Definición de vectores a partir de rangos¶

Una de las utilizaciones de los rangos es la definición de las componentes de un vector. Veamos un ejemplo. Supongamos que queremos un vector denominado az cuyas componentes sean los ángulos 0, 90, 180, 270. Podríamos definirlo de la siguiente manera:

rango = 0:90:270
az = [rango]
% El resultado será: az=[0 90 180 270]

De hecho, no es necesario definir primero el rango y luego el vector. La manera habitual de hacerlo es definiendo el rango directamente dentro de los corchetes:

az = [0:90:270]
% El resultado será: az=[0 90 180 270]

La función linspace()¶

La función linspace() es también muy útil para definir vectores. Devuelve un vector con n números equiespaciados entre unos números de inicio y final. Los argumentos que se le pasan a la función son el número de inicio de la serie, el número final de la serie y el número total de elementos que tendrá la serie:

% v = linspace( inicio, fin, num_elementos)
v = linspace(1,10,10)
% El resultado será: v=[1 2 3 4 5 6 7 8 9 10]

La función linspace() permite también realizar series descendentes, y también podemos utilizar expresiones válidas de Octave para los parámetros:

w = linspace(10,5,5)
% El resultado será: w=[10 9 8 7 6]
z = linspace(0,pi/2,4)
% El resultado será: w=[0.00000 0.52360 1.04720 1.57080]

Lectura de las componentes de un vector¶

Cuando queramos referirnos a una componente concreta de un vector, utilizaremos el nombre del vector seguido del índice de la componente encerrado entre paréntesis. En Octave, el índice de la primera componente de un vector es el 1. Veamos un ejemplo refiriéndonos al vector definido en el epígrafe anterior az = [0 90 180 270] :

az(1)
% El resultado será: 0
az(2)
% El resultado será: 90
az(4)
% El resultado será: 270

Octave proporciona un truco para extraer la última componente de un vector, sin tener que explicitar el índice, y es utilizar la palabra end. En el ejemplo anterior es equivalente escribir az(4) que az(end):

az(4)
% El resultado será: 270
az(end)
% El resultado será: 270

También podemos definir las componentes a extraer mediante un rango. En el siguiente ejemplo extraemos la segunda y tercera componente del vector az. El resultado de la operación es un vector de dos componentes, [90 180], que se lo asignamos a un nuevo vector que hemos llamado es:

es = az(2:3)
% El resultado será: es=[90 180]

Otro ejemplo de extracción de un rango de componentes de un vector es el siguiente, en el que primeramente definimos un vector v y a continuación extraemos algunas de sus componentes a un nuevo vector W:

v = [1:2:10] % Se define un vector y se asigna a la variable v
% El resultado será: v=[1 3 5 7 9]
w = v(1:2:5) % Se extráen ciertas componentes del vector v que se asignan a w
% El resultado será: w=[1 5 9]

Concatenación de vectores¶

Podemos unir dos vectores, esto es, poner las componentes de un vector a continuación de las de otro para formar un único vector. Para ello se definen entre corchetes los elementos a concatenar separados por coma o por espacio. Veamos un ejemplo:

a = [1 2 3];
b = [4 5 6];
c = [a b]
% El resultado será: c = [1 2 3 4 5 6]

Podemos hacer otro tipo de concatenaciones, por ejemplo, concatenar un vector con un rango, o un número con un vector, o un vector con un número. De hecho, cundo definimos un vector tecleando sus componentes, lo que hacemos es concatenar varios números:

a = [1 2 3];
c = [a 5]
% El resultado será: c = [1 2 3 5]
d = [-5:0 a 4:5]
% El resultado será: d = [-5 -4 -3 -2 -1 0 1 2 3 4 5]
d = [-5 -4 -3 -2 -1 0 a 4 5]
% El resultado será: d = [-5 -4 -3 -2 -1 0 1 2 3 4 5]

Se pueden concatenar de manera similar vectores columna. En este caso las distintas subcomponentes se separarán por punto y coma.

Funciones utilitarias para trabajar con vectores¶

Hay varias funciones útiles a la hora de trabajar con vectores. Siendo v un vector:

• length(v) Devuelve el número de componentes del vector v
• max(v), min(v) Devuelve el valor máximo/mínimo de entre las componentes de v
• sum(v), prod(v) Devuelve la suma/producto de las componentes de v
• norm(v) Devuelve el módulo del vector v
• dot(v, w) Devuelve el producto escalar de los vectores v y w
• cross(v, w) Devuelve el vector producto vectorial de v y w (Las dimensiones de v y w deben ser congruentes como máximo tres)
• sort(v) Devuelve un vector con las componentes de v ordenadas de menor a mayor

Funciones para cadenas de caracteres¶

• strrep : Reemplaza una subcadena dentro de una cadena

cad = '3,78'
newcad = strrep(cad, ',','.');
% result : newcad = 3.78

–strtrim(cad): Quita los caracteres espacio del inicio y final de una cadena

• C = strsplit(str, delimiter): Devuelve un CellArray con las subcadenas que resultan de dividir la cadena original str con el separador delimiter.

cad = '3,78'
c = strsplit(cad, ',');
pentera = str2num(c{1}) % result 3 double
pdec = str2num(c{2}) % result 78 double

• str2num(cad): Convierte en número double la cadena cad

La función rand()¶

La función rand() se utiliza para generar numeros pseudoaleatorios. Se pueden generar scalares o matrices de números pseudoaleatorios. Si llamamos a la función sin argumentos nos devolverá un numeros pseudoaleatorios uniformemente distribuidos entre 0 y 1.

Si pasamos un número entero n como parámetro, la función nos devolverá una matriz cuadrada de dimensión n de numeros pseudoaleatorios comprendidos entre 0 y 1. Si le pasamos dos números enteros m y n la matriz resultante tendrá m filas y n columnas.

El numero que está utilizando Octave como semilla generadora de los números pseudoaleatorios se puede ver con rand(‘seed’). También podemos modificar la semilla con rand(‘seed’, x), donde x es la nueva semilla que queremos utilizar.

Los números generados por rand() son doubles comprendidos entre 0 y 1. Si queremos generar números comprendidos en un determinado intervalo debemos operar los números obtenidos de rand() para adaptarlos al intervalo buscado:

M = 1
rand() + M
% El resultado estará comprendido entre M y M+1 (en este caso entre 1 y 2)
%
M = 10
rand() * M
% El resultado estará comprendido entre 0 y M
%
A = 10
B = 20
A + (B-A)*rand()
% El resultado estará entre A y B

Para generar numeros enteros en un determinado intervalo se puede utilizar la función randi():

randi(10)
% Genera un número entero entre 1 y 10

También se pueden conseguir enteros mediante la función rand():

A = 10
B = 20
round(A + (B-A)*rand())
% El resultado será un entero (double con parte decimal cero) entre A y B
int32(round(A + (B-A)*rand()))
% El resultado será un entero int32 entre A y B

Utilidades para medir el tiempo¶

La función clock() nos devolverá un vector con el año, mes, día, hora, minutos y segundos:

Las función date() nos devuelve la fecha como una cadena de texto.

La función time() nos devuelve el número de segundos transcurridos desde las 00:00 horas UTC del día 1 de enero de 1970.

La función etime() (elapsed time) nos devolverá el tiempo transcurrido entre dos fechas:

t0 = clock();
# unas cuantas computaciones más tarde...
tiempo_transcurrido = etime(clock(), t0);

La función is_leap_year() nos devolverá 1 (verdadero) si el año pasado como argumento es bisiesto.

Scripts¶

Podemos escribir una serie de instrucciones de Octave en un fichero con extensión .m y luego ejecutarlo desde la consola. Al fichero con las instrucciones de Octave se le llama un script o programa. Para ejecutarlo tendremos que situarnos en el directorio donde esté el fichero y teclear su nombre sin extensión en la consola de Octave.

El lenguaje m es un lenguaje interpretado, esto es, Octave lee los ficheros m y los va ejecutando línea a línea, sin una compilación previa.

Para escribir los programas utilizaremos un editor de texto. La versión actual de Octave trae un editor incorporado y también copia en el directorio de instalación una versión del editor Notepad++.

El programa Hola mundo¶

Es tradicional comenzar el aprendizaje de cualquier lenguaje de programación con el típico programa Hola Mundo. Se trata de escribir un programa que muestre por pantalla la frase Hola mundo. El ejercicio permite comprobar que uno sabe crear el fichero de instrucciones, escribir alguna instrucción en él y hacer que el ordenador ejecute el programa.

En nuestro caso vamos a crear un fichero en el mismo directorio de trabajo en el que nos encontremos, llamado ‘hello.m’, y vamos a escribir en él una única instrucción:

disp('Hola, mundo');

La función disp() nos permite mostrar una variable por pantalla, en este caso una cadena de texto. Tendremos que grabar el archivo en el disco con el nombre ‘hello.m’ y después, desde la consola, teclear hello. Si todo va bien, el programa se ejecutará, nos mostrará el mensaje ‘Hola mundo’ y finalizará, con lo que volveremos a tener el prompt de la consola a disposición para seguir tecleando instrucciones.

Desde el editor integrado en la versión gráfica de Octave podemos crear un nuevo fichero, escribir en él la instrucción y grabarlo con el nombre elegido.

Líneas de Comentarios¶

Las líneas de comentarios son líneas que no serán ejecutadas por Octave al correr el programa. Sirven para documentar las rutinas.

Los comentarios pueden ser de una línea. Para ello se deben preceder del símbolo ‘#’ o del símbolo ‘%’. El comentario será desde el símbolo comentario hasta el final de la línea. Puede ocupar toda la línea, si el símbolo comentario se pone al principio de la misma, o comentar una línea concreta, desde el símbolo comentario hasta el final de la línea. Veamos un ejemplo:

# Cuenta atrás para cohetes
disp (3);
disp (2);
disp (1);
disp ("Blast Off!");  # El cohete abandona la plataforma

Se pueden comentar bloques de código. Para ello hay que encerrar el bloque a comentar entre los símbolos ‘#{‘ y ‘#}’. Estos símbolos tienen que aparecer cómo únicos caracteres en su línea. ‘#{‘ aparecerá como línea previa al código a comentar, y ‘#}’ como linea posterior. Una utilidad de lo anterior es ‘desactivar’ bloques de código durante el proceso de depuración de los programas, por ejemplo.

# Cuenta atrás para cohetes
disp (3);
#{
disp (2);
disp (1);
#}
disp ("Blast Off!");  # El cohete abandona la plataforma

En esta versión del programa las líneas de disp(2) y disp(1) no se ejecutarán, y se pasará directamente del disp(3) al lanzamiento.

Comentar una línea de código es habitual durante la fase de creación y depuración de los programas. Ello nos permite deshabilitar ciertas partes del código para estudiar el funcionamiento del resto.

La instrucción input¶

La instrucción input() nos permite solicitar al usuario algunos datos que se necesitan para los cálculos ulteriores. Por ejemplo, si queremos hacer un programa que devuelva el valor de elevar un número al cubo podríamos hacer lo siguiente:

x = input('Teclée un número: ');
y = x^3;
disp(y)

La función input() recibe una cadena de texto a modo de mensaje que muestra por consola, situando el cursor a continuación a la espera de que el usuario teclée el valor solicitado. Una vez el usuario pulsa la tecla intro, el valor tecleado será asignado a la variable elegida, en el ejemplo anterior la variable x. El programa, entonces, se seguirá ejecutando en la línea siguiente a la de la sentencia input y podrá utilizar el valor de x en los cálculos posteriores. En el ejemplo, se crea una nueva variable y que almacena el valor de x^3 y a continuación se muestra el valor de y por pantalla mediante la sentencia disp().

Ejercicio 1 (Función size())¶

Escribir una función que devuelva logical 1 si la variable que recibe como argumento es un vector fila o columna, y logical 0 en otro caso ( argumento escalar o argumento matriz).

function is = isVector(v)
% Devuelve logical 1 si la variable recibida
% es un vector y logical 0 en otro caso.

   % Comprobación del número de argumentos
   if (nargin ~= 1)
      is = logical(0);
      return;
   end

   if (xor(size(v,1)==1, size(v,2)==1))
      % Es un vector
      is = logical(1);
   else
      % No es un vector
      is = logical(0);
      return;
   end

end

Ejercicio 2 (Bucles for)¶

Escribir una función que calcule la suma s de los elementos de un vector que recibe como argumento. La función comprobará que el argumento recibido es un escalar o un vector. En caso de recibir un argumento no válido, la función devolverá NaN.

function s = sumvector(v)
% Suma las compnentes de un vector

   if(nargin ~= 1)
      s = NaN;
      return;
   end

   if( (size(v,1) ~= 1) & (size(v,2) ~= 1) )
      s = NaN;
      return;
   end

   s = 0;
   for index = 1: length(v)
      s = s + v(index);
   end

end

Ejercicio 3 (Bucles for. Algoritmo de la suma)¶

Escribir una función que devuelva la suma de todos los elementos de un vector o matriz que recibe como argumento. En el caso de recibir un escalar o un vector la función pedida sería equivalente a la función sum() de la librería. En el caso de recibir una matriz como argumento, el resultado de la función pedida sería equivalente a la expresión sum(sum()). En la solución no se utilizarán funciones de librería, solamente operadores aritmeticos sencillos.

function s = sumall(A)
   % Suma todos los elementos de un vector o una matriz

   filas = size(A,1);
   cols = size(A,2);

   s = 0;
   for ind1 = 1: filas
      for ind2 = 1: cols
         s = s + A(ind1,ind2);
      end
   end
end

Ejercicio 4 (Bucles for. . Algoritmo del producto)¶

Escribir una función que devuelva el producto de todos los elementos de un vector o matriz que recibe como argumento. En el caso de recibir un escalar o un vector la función pedida sería equivalente a la función prod() de la librería. En el caso de recibir una matriz como argumento, el resultado de la función pedida sería equivalente a la expresión prod(prod()). En la solución no se utilizarán funciones de librería, solamente operadores aritmeticos sencillos.

function p = prodall(A)
% Multiplica todos los elementos de un vector o una matriz

   filas = size(A,1);
   cols = size(A,2);

   p = 1;
   for ind1 = 1: filas
      for ind2 = 1: cols
         p = p * A(ind1,ind2);
      end
   end
end

Ejercicio 5 (Bucles for)¶

Escribir una función que realice el puntoproducto .* de dos vectores que recibe como argumento. La función comprobará que el número de argumentos recibidos es correcto, y que se trata de dos vectores. Si los argumentos no son correctos, la función devolverá NaN. En la solución solo se permitirá la utilización de bucles y operadores aritméticos sencillos. La función devolverá el resultado correcto independientemente de que los vectores recibidos sean fila-fila, col-col, fila-col o col-fila.

function pp = puntoprod(v, w)
% Calcula el punto producto de dos vectores que recibe como argumento

% Comprobación del número de argumentos recibidos
if (nargin ~= 2)
   pp = NaN;
   return;
end

% Comprobación de la dimensión de los vectores
if (size(v,1) ~= 1 & size(v,2) ~= 1)
   pp = NaN;
   return;
elseif (size(w,1) ~= 1 & size(w,2) ~= 1)
   pp = NaN;
   return;
elseif (length(v) ~= length(w))
   pp = NaN;
   return;
end

pp = 0;
for index = 1:length(v)
   pp = pp + v(index) * w(index);
end

Ejercicio 6 (Bucles for)¶

Escribir una función que calcule la matriz C producto de dos matrices A y B que recibe como argumentos. La función comprobará, en primer lugar, que el número de argumentos recibidos es el adecuado y que las dimensiones de las matrices A y B son congruentes. En caso de argumentos incorrectos la función devolvera una matriz vacía []. Si los argumentos son correctos, la función calculará la matriz C sin hacer uso de funciones de librería, utilizando bucles y operadores aritméticos sencillos.

function C = prodmat(A, B)
% prodmat.m : Calcula la matriz C, producto de las matrices A y B

% Comprobación del número de argumentos recibidos
if (nargin ~= 2)
   C = [];
   return;
end

% Comprobación de la dimensión de las matrices
% El número de columnas de A debe ser igual al de filas de B
% Si las dimensiones no son las correctas el programa finaliza
if (size(A,2) ~= size(B,1))
   C = [];
   return;
end

% Preparo la matriz solución con las dimensiones adecuadas
% Si son A(m,n) y B(n,p) las dimensiones de C seran C(m,p)
C = zeros(size(A,1),size(B,2));

% Cálculo del producto mediante dos bucles for anidados
for fila = 1: size(A,1)
   for col = 1:size(B,2)
      sum = 0;
      for(index=1:size(A,2))
         sum = sum + A(fila, index) * B(index, col);
         end
         C(fila,col) = sum;
      end
   end

end

Ejercicio 7 (Strings, Funciones)¶

La variable contiene los resultados de los catorce partidos de futbol de la quiniela de una jornada. La variable tiene los resultados apostados en una determinada quiniela. Se pide: (a) Calcule en una sola expresión cuantos resultados hay acertados. (b) Cuales son los numeros de orden de los partidos que se han acertado. (C) Cuáles son los resultados acertados.

numAciertos = sum(result == apuesta)
ans = 8

partidosAcertados = find(result == apuesta)
ans =    1    2    5    6    8    9   12   13

resultadosAcertados = result(result == apuesta)
ans = 1121X11X

Ejercicio 8 (Función rand())¶

Desarrolle una función llamada dado() que devuelva un número aleatorio entre uno y seis, todos ellos con la misma probabilidad de ocurrencia.

function r = dado()
   % dado() : devuelve un numero r al azar entre 1 y 6
   % Fórmulas útiles para números aleatorios
   % Número entre A y A+1 : A + rand();
   % Número entre A y B : A + (B-A) * rand()
   % Número entero entre A y B : round(A + (B-A) * rand())
   n = round(1 + (6-1) * rand());
end

Ejercicio 9 (nargin, función rand())¶

Modifique la función dado() del ejercicio anterior de manera que en caso de recibir un número n como argumento devuelva un vector de n componentes cada una de ellas con el resultado de una tirada de dado. Si no recibe argumentos seguirá devolviendo un único número entre 1 y 6. Si el parámetro n recibido como argumento no es entero, se redondeará al entero más proximo.

function r = dado(n)
   % dado() : devuelve un numero r al azar entre 1 y 6
      if(nargin == 0)
         r = round(1 + (6-1) * rand());
      return;
   end

   if(size(n,1) ~=1 | size(n,2) != 1 | n<=0)
      r = round(1 + (6-1) * rand());
      return;
   end

   n = round(n);
   r = [];
   for ind = 1 : n
      r(ind) = round(1 + 5 * rand());
   end
end

Ejercicio 10 (Algoritmo de máximo)¶

Desarrolle una función que reciba un vector V de números reales y devuelva un vector de dos componentes [max, n], la primera, max, con el valor máximo de las componentes del vector V de entrada y la segunda, n, con la posición que ocupa el máximo dentro del vector V. El problema se resolverá mediante bucles y el algoritmo de máximo.

function [max, pos] = maxx(V)
   max = V(1);
   pos = 1;
   for k = 1 : length(V)
      if V(k) > max
         max = V(k);
         pos = k;
      end
   end
 end

Ejercicio 11 (Diagonal secundaria)¶

Desarrolle una función que reciba una matriz como argumento y devuelva un vector fila con los elementos de la diagonal secundaria.

function dg = diag2(A)
% diag2(A) : Devuelve un vector fila con las componentes de la diagonal secundaria

   numfilas = size(A,1);
   numcolumnas = size(A,2);

   % El primer elemento de la diagonal secundaria es A(1,numcolumnas)
   % El segundo es A(2, numcolumnas-1)
   % El último es A(numfilas, numcolumnas-numfilas) si numcolumnas > numfilas
   % El último es A(numfilas, 1) si numcolumnas <= numfilas
   dg = [];
   k = 0;
   while ((numcolumnas - k) > 0 &  (k+1) <= numfilas )
      dg(k+1) = A(k + 1, numcolumnas - k );
      k = k +1;
   end
end

Otra solución del problema utilizando un bucle for en lugar de un bucle while:

function dg = diag2(A)
% diag2(A) : Devuelve un vector fila con las componentes de la diagonal secundaria

   numfilas = size(A,1);
   numcolumnas = size(A,2);

   dg = [];
   for k = 1 : numfilas
      if (numcolumnas - k + 1) == 0
         break;
      end
      dg(k) = A(k, numcolumnas - k +1);
   end
end

Otra solución, observando que los elementos de la diagonal secundaria siempre suman 1 + numcolumnas:

function dg = diag23(A)
% diag2(A) : Devuelve un vector fila con las componentes de la diagonal secundaria

   numfilas = size(A,1);
   numcolumnas = size(A,2);

   % Los índices de los elementos de la diagonal secundaria
   % siempre suman 1 + numcolumnas

   dg = [];
   contador = 0;
   for fil = 1 : numfilas
      for col = 1 : numcolumnas
         if (fil + col) == (1 + numcolumnas)
            contador = contador + 1;
            dg(contador) = A(fil, col);
         end
      end
   end
end

Ejercicio 12 (Recursividad, Factorial de un número)¶

Podemos calcular el factorial de un número mediante la función factorial(), también a partir de la función prod(1:n). Lo que se pide en este problema es desarrollar una función que calcule el factorial de un número natural mediante recursividad, esto es, mediante una función que se llame a si misma.

function fac = recurfac(n)
% Calcula el factorial de un número n mediante recursividad
   if(n == 0)
      fac = 1;
   else
      fac = n * recurfac(n-1);
   end
end

Gráficos de curvas en coordenadas polares¶

En coordenadas polares, un punto del plano queda definido mediante el valor de la distancia desde el punto al origen de coordenadas y el valor del ángulo que va desde el eje , en sentido antihorario, hasta la línea que une el origen de coordenadas con el punto (ver figura).

Se le llama ecuación polar a la ecuación que define una curva expresada en coordenadas polares. En muchos casos se puede especificar tal ecuación definiendo como una función de . La curva resultante consiste en una serie de puntos en la forma . Cuando una función está expresada en coordenadas polares, al origen de coordenadas se le denomina también polo.

Para transformar las coordenadas polares a coordenadas cartesianas hay que aplicar las siguientes ecuaciones:

Vamos a explicar cómo representar gráficamente con Octave una función expresada en coordenadas polares. Utilizaremos para ello la siguiente función:

Lo primero es discretizar el parámetro . A continuación se obtiene el vector de valores del radio . Una vez que tenemos los dos vectores con los valores de y podemos obtener los vectores con los valores de e , haciendo la transformación a coordenadas cartesianas. Finálmente utilizaremos la función plot() con los valores de e obtenidos. El código Octave es el siguiente:

% Discretizamos ro y theta
theta = linspace(0, 10*pi, 300);
ro = 3 - 5 .* (theta);

% Transformamos a cartesianas
x = ro .* cos(theta);
y = ro .* sin(theta);

% Ploteamos
plot(x,y)
axis equal

Algunos ejemplos de funciones en coordenadas polares son:

Circunferencia: la expresión general de una circunferencia con centro en el punto y radio es:

Cuando la circunferencia tiene centro en el origen de coordenadas, la expresión se simplifica:

Rosa polar: es una expresión matemática con forma de flor. Su expresión es:

donde es cualquier valor constante, incluido el cero. Si es un valor entero, si es impar, la flor tendrá pétalos y si es par, la flor tendrá pétalos. Si es racional pero no entero, la gráfica es similar a una rosa pero con los pétalos solapados. Nótese que estas ecuaciones nunca definen una rosa con 2, 6, 10, 14, etc. pétalos. La variable a representa la longitud de los pétalos de la rosa.

Espiral de Arquímedes: La espiral de Arquímedes es una famosa espiral descubierta por Arquímedes, la cual puede expresarse también como una ecuación polar simple. Se representa con la ecuación:

Un cambio en el parámetro producirá un giro en la espiral, mientras que controla la distancia entre los brazos, la cual es constante para una espiral dada. La espiral de Arquímedes tiene dos brazos, uno para y otro para . Los dos brazos están conectados en el polo. La imagen especular de un brazo sobre el eje vertical produce el otro brazo. Esta curva fue una de las primeras curvas, después de las secciones cónicas, en ser descritas en tratados matemáticos. Además es el principal ejemplo de curva que puede representarse de forma más fácil con una ecuación polar.

Elipse: la expresión polar de una elipse centrada en el origen de coordenadas y de semiejes y es:

Se define la excentricidad de una elipse mediante la siguiente expresión:

La ecuación polar de la elipse en función de su excentricidad queda de la siguiente manera:

La función plot3()¶

• plot3(): La función plot3() dibuja curvas en tres dimensiones, o ternas de valores arbitrarios sin necesidad de que formen una superficie. Por ejemplo, el código siguiente dibuja una espiral en tres dimensiones:

t = 0:0.1:10*pi;
r = linspace (0, 1, length(t));
z = linspace (0, 1, length(t));
plot3 (r.*sin(t), r.*cos(t), z, 'linewidth',2);

• view(az, el): Situa el punto de vista en el azimut az y la elevación el. Haciendo view(2) o view(3) establece los valores por defecto para dos y tres dimensiones.

La función mesh()¶

La funcion mesh() dibuja un gráfico de malla en la forma Z=f(X,Y), tomando la Z como elevación y haciendo el color del mallado proporcional al valor de Z. Si Z es una matriz cuyas dimensiones son [m, n]. El vector de las coordenadas x tendrá que tener m elementos, y el vector de las componentes y tendrá n elementos.

Hay varias formas de llamar a la función mesh():

• mesh(Z): dibuja el gráfico de malla 3D correspondiente a Z=Z(x,y) haciendo x=1:m e y=1:n con [m, n] = size(Z). El color es proporcional al valor de Z.
• mesh(X, Y, Z): dibuja el gráfico de malla 3D con el color determinado por el valor de Z. Si X e Y son vectores las dimensiones tienen que cumplir: m=length(X), n=length(y) y [m, n]=size(Z). En este caso, los puntos de la superfcie serán (X(j), Y(i), Z(i,j)). Sy X, Y, Z son matrices, sus dimensiones serán iguales. En este caso los puntos del mallado serán: (X(i,j), Y(i,j), Z(i,j))
• mesgrid(x,y): Dados dos vectores de coordenadas x e y, la función meshgrid() devuelve una pareja de matrices Mx, My que se corresponden con la red de puntos de intersección de las dos series x e y. Estas matrices se pueden pasar directamente a la función mesh() y también se pueden utilizar para calcular los valores de Z si conocemos su expresión en la forma Z=f(X,Y).

Si mallamos los vectores x=[0 1 2 3] e y=[0 1 2 3] con la función [xx, yy]=meshgrid(x,y), el resultado será el siguiente:

xx =
      0 1 2 3
      0 1 2 3
      0 1 2 3
      0 1 2 3
yy =
      0 0 0 0
      1 1 1 1
      2 2 2 2
      3 3 3 3

Por ejemplo:

tx = linspace (-8, 8, 41)';
ty = linspace (-8, 8, 41)';
[xx, yy] = meshgrid (tx, ty);
r = sqrt (xx .^ 2 + yy .^ 2) + eps;
tz = sin (r) ./ r;
mesh (tx, ty, tz);

Produce el gráfico del conocido sombrero que se ve en la figura:

• meshc(): Es similar a mesh() pero dibuja además las curvas de nivel sobre el plano x-y: El resultado se puede ver en la siguiente figura:

Para representar una superficie expresada en ecuaciones paramétricas se procede de la siguiente forma. Supongamos que las ecuaciones de la superficie expresada en función de dos parámetros son las siguientes:

Generaremos unos rangos para el valor de cada parámetro y utilizaremos a continuación la función meshgrid() para generar las matrices con los mallados de valores de los parámetros. A continuación podemos calcular las matrices que servirán para pasarle a la función mesh().

Veamos un ejemplo. Sea el paraboloide elíptico definido por las siguientes ecuaciones:

Se pide representarlo para

El resultado se puede ver en la siguiente figura:

Comandos utilitarios trabajando con ficheros¶

Cuando estamos trabajando con ficheros de texto hay varios comandos de Octave que conviene conocer.

El comando type permite mostrar el contenido de un fichero de texto en la consola. La forma de utilizar el comando es la siguiente: type filename, siendo filename el nombre de un fichero existente. Si el nombre del fichero incluye su ruta de directorio, se buscará el fichero en la ruta especificada. Si solo se incluye el nombre del fichero, se buscará en el directorio de trabajo.

El directorio de trabajo en el que estamos posicionados se puede consultar con pwd(), que muestra en pantalla la ruta del directorio de trabajo donde estamos posicionados.

El contenido del directorio de trabajo se puede mostrar mediante el comando dir, al estilo windows, o mediante el comando ls, al estilo linux. Los dos comando admiten una mascara que permite especificar el tipo de fichero buscados. Por ejemplo podríamos utilizar dir *.m para mostrar un listado de los ficheros con extensión .m. En la máscara, el arterisco hace la función de comodín. Podríamos leer el comando anterior como: ‘dir lo-que-sea punto m’.

También podemos cambiar el directorio de trabajo desde la consola de Octave. Para ello se utiliza el comando cd, ‘change directory‘, que admite un parámetro a continuación para especificar la ruta a la que se quiere acceder. Hay un atajo que hay que conocer. Si tecleamos cd .., cd punto punto, nos moveremos al directorio padre del directorio en el que estemos posicionados. Para acceder a un determinado directorio utilizaremos: cd dirname, donde dirname es el nombre o la ruta completa del directorio en el que queremos posicionarnos.

Los comandos anteriores admiten ser utilizados en forma de función.

Los comandos load y save¶

Podemos guardar variables en ficheros de texto mediante el comando Octave save. Las variables guardadas las podremos recuperar luego mediante el comando load. El comando save se puede utilizar de distintas maneras:

• save filename: Guarda todas las variables del workspace en el fichero filename, donde filename es el nombre del fichero y opcionalmente su ruta completa.
• save filename v1 v2: Se pueden especificar qué variables concretas queremos guardar en el fichero. En este ejemplo, solo se guardarán el el fichero filename las variables de nombre v1 y v2.
• save options filename v1 v2: donde options son una o más opciones separadas por espacio. Las opciones nos permiten especificar, por ejemplo, que se guarden los ficheros en distintos formatos binarios o que se guarden los datos en formato comprimido zip.

La forma de operar se puede ver en el siguiente ejemplo, donde primero definimos una matriz y a continuación la guardamos en un fichero ‘matriz.txt’. Por último se muestra el contenido del fichero mediante type():

La función fopen()¶

Para poder realizar operaciones de lectura y/o escritura en los ficheros de texto es necesario primero abrir el fichero con fopen(). Una vez abierto el fichero, se procede a realizar las distintas operaciones de lectura y/o de escritura en el mismo. Cuando se termina de operar con el fichero hay que cerrarlo con la funcion fclose().

El procedimiento habitual para leer datos de un fichero se resume en el siguiente esquema:

file = fopen('mifichero.txt', 'r');

% Operaciones de lectura del fichero

fclose(fid);

Ĺa función fopen() tiene la siguiente signatura:

• file = fopen(filename, permission) Abre el fichero de nombre filename con el tipo de acceso especificado por permission. Devuelve file, que es un número, un puntero, que identifica al fichero en ulteriores operaciones de lectura escritura. Si se producen errores al intentar abrir el fichero, la función fopen() devuelve -1. El puntero será un número mayor o igual a 3. Los parametros son los siguientes:

  • filename: cadena de texto con el nombre del fichero que se quiere abrir. Si incluye la ruta del fichero se utilizará, si no se buscará o creara el fichero en el directorio de trabajo.

  • permission: especifica el modo de apertura del fichero. Lo habitual es abrir el fichero para escribir en él o para leer de él, si bien es también posible abrir el fichero en modo lectura y escritura. En los modos de escritura es posible añadir texto a un fichero existente o bien crear uno nuevo para la operación, sobreescribiendo el existente si lo hubiera. permission es una cadena de texto que puede tener los siguientes valores:

    • ‘r’: Abre el fichero para lectura. Es el modo por defecto, si se utiliza la función fopen() sin el argumento permission
    • ‘r+’: Abre el fichero en modo lectura-escritura
    • ‘w’: Abre o crea un nuevo fichero en modo escritura. Si existe se sobreescribe
    • ‘w+’: Abre o crea un fichero para lectura-escritura. Si existe, se sobreescribe
    • ‘a’: Abre o crea un nuevo fichero para escritura. Si existe el fichero, añade al final del mismo.
    • ‘a+’: Abre o crea un nuevo fichero para lectura-escritura. Si existe, se añade al final del mismo.

La función fclose()¶

Sirve para cerrar los ficheros cuando se ha terminado de trabajar con ellos.

• fclose(file): Cierra el fichero file. Devuelve 0 si pudo cerrar el fichero o -1 si hay errores.
• fclose(‘All’): Cierra todos los ficheros que haya abiertos. Devuelve 0 si pudieron cerrar los ficheros o -1 si hay errores.

Lectura de ficheros línea a línea¶

Una vez abierto el fichero para lectura podemos ir leyendo una línea cada vez, que recibiremos como una cadena de texto, con o sin carácter fin de línea, según la función de lectura utilizada:

• line = fgets(file): Lee una línea del fichero file, incluyendo el caracter de fin de línea, y la devuelve en forma de cadena de texto.
• line = fgets(file, nchar): Lee al menos nchar caracteres de la siguiente línea del fichero file. (Si se alcanza el final de línea o el final del fichero devuelve lo que haya leído hasta ahí).
• line = fgetl(file): Lee la siguiente línea del fichero file, sin incluir el caracter fin de línea, y la devuelve como una cadena de caracteres.

Lectura de datos formateados¶

Octave ofrece la posibilidad de leee los datos contenidos en el fichero, uno a uno, y devolverlos en el formato seleccionado, pudiendo de esta manera leer , no solo cadenas de texto, sino también números enteros o doubles.

La función para leer datos formateados es fscanf(). La función fscanf() admite varias signaturas. Una forma habitual de utilizar la función es la siguiente:

fscanf()

Escritura de datos formateados¶

En el ejemplo siguiente creamos una matriz A de tres filas y cuatro columnas. A continuación la escribimos en un fichero con una simple indicación de formato numérico. Cuando hacemos type al fichero recién creado podemos ver la forma en que Octave ejecuta la función fprintf().

Como vemos la sentencia fprintf(file, ‘%d ‘, A) escribe todos los elementos de la matriz A en el fichero, utilizando el formato indicado: numero entero más espacio. Hay que fijarse en el orden en que Octave coge los elementos de A para escribirlos en el fichero. Vemos que coge los elementos de la matriz por columnas, esto es, primero el A(1,1), a continuación el A(2,1), y así hasta agotar la primera columna, leyendo a continuación las siguientes columnas de la misma manera hasta agotar la matriz.

Lo que hace Octave es leer las matrices en el orden el que las tiene guardadas en memoria, y Octave guarda las matrices en memoria por columnas.

Si en lugar del formato especificado en el ejemplo anterior, ejecutamos la sentencia fprintf(file, ‘%d %d %d %d \n’, A), en la que decimos a Octave que cada cuatro números inserte un carácter fin de línea ‘\n’, Octave escribirá tres filas y cuatro columnas, pero no en el orden de filas y columnas en el que tenemos definida la matriz:

Una vez más Octave rellena el formato indicado cogiendo los elementos de la matriz por columnas.

Utilización de paquetes adicionales¶

Octave-Forge [1] es el repositorio donde se centraliza el desarrollo colaborativo de paquetes de funcionalidades adicionales para Octave.

Los paquetes de Octave-Forge amplian la funcionalidades estándar de Octave, añadiendo utilidades y funciones específicas para distintos ambitos de aplicación. Todos los paquetes de Octave-Forge se pueden instalar mediante el mecanismo de instalación de paquetes incorporado en el propio Octave. Para instalar un paquete, desde la consola de Octave, teclear:

pkg install -forge package_name

donde package-name es el nombre del paquete a instalar. Octave descargará e instalará el paquete, quedando las funciones contenidas en el mismo disponibles para su utilización.

Hay multitud de paquetes con funcionalidades adicionales que se pueden incorporar a la instalación de Octave. La lista de paquetes disponibles se puede consultar en [2]. Algunos ejemplos de paquetes adicionales que podemos instalar:

• image: procesamiento de imágenes
• io: funciones de entrada-salida en distintos formatos: csv, xml, xls, json y otros.
• java: interface para utilizar funciones java
• mechanics: incluye, entre otras, funciones útilitarias para calculo de estructuras
• octproj: transformaciones entre proyecciones cartográficas
• splines: funciones para manipulación de splines

Podemos comprobar la lista de paquetes que tenemos instalados haciendo:

pkg list

Además de la lista de paquetes anterior, hay una lista adicional de paquetes no incluidos en el sistema de paquetes pkg, pero cuyo código fuente se puede consultar y descargar en [3].

[1] Octave-Forge: http://octave.sourceforge.net/index.html

[2] Octave-Forge: Lista de paquetes: http://octave.sourceforge.net/packages.php

[3] Paquetes adicionales: http://sourceforge.net/p/octave/_list/hg

Utilidades para interactuar con el Sistema Operativo¶

Podemos ejecutar comandos del sistema operativo mediante la función system(cmd, flag). El primer argumento será una cadena de caracteres con el comando a ejecutar. El segundo parámetro es opcional, y si está presente, la salida del comando se devolverá como una cadena de texto. Si este segundo parámetro no está presente la salida del comando se enviará a la salida estándar. El valor de este segundo parámetro es cualquiera, por ejemplo un ‘1’.

Podemos recuperar el valor de una variable del sistema mediante la función getenv(‘var_name’), por ejemplo, en Linux, para recuperar el path:

Podemos conocer el directorio actual mediante pwd(), cambiar de directorio mediante cd <dir_name>, listar el contenido de un directorio mediante dir o ls y crear un nuevo directorio mediante mkdir <di_name>.

La función computer() nos informará del ordenador en el que se está ejecutando Octave.

El copo de nieve de Koch¶

La curva fractal ‘Koch snowflake‘, también conocida como la ‘isla de Koch’, fue descrita por primera vez por Helge von Koch en 1904. Se construye partiendo de un triángulo equilátero. En cada lado del triángulo se sustituye la tercera parte central por un nuevo triángulo equilátero de lado un tercio del primero. Este proceso se repite indefinidamente.

El resultado es una curva cuya longitud tiende a infinito, pero encerrando un área que tiende a una cantidad finita. Podemos implementar una función que dibuje una curva de Koch mediante la utilización de la tortuga creada en el apartado anterior. Para ello crearemos en primer lugar la función drawkoch() que dibuja, de forma recursiva, el segmento básico de la curva de Koch para un determinado nivel de profundidad:

function drawkoch(length, depth)
   if depth == 0
      turtle('AVANZA', length);
   else
      drawkoch(length/3, depth-1);
      turtle('GIRA', 60);
      drawkoch(length/3, depth-1);
      turtle('GIRA', -120);
      drawkoch(length/3, depth-1);
      turtle('GIRA', 60);
      drawkoch(length/3, depth-1);
   end
end

Podemos probar la función recién creada enlazando los tres primeros niveles de segmentos de Koch. La secuencia de instrucciones y el resultado se puede apreciar en la siguiente figura:

Para dibujar el copo de nieve completo, con un determinado nivel de profundidad tendremos que dibujar tres tramos de Koch, cada uno girado -120 grados respecto del anterior. Podemos implementar una función que denominaremos koch(), que utiliza la función anterior para dibujar el copo de nieve. El código de la función koch() será el siguiente:

function koch(length, level)
% Dibuja la curva de Koch para un nivel determinado

   hold off;
   close;
   hold on;
   xlim([-10, 2*length*sind(60)+10]);
   ylim([-10-length*sind(60), 10+length*sind(60)]);

   global x y angle;
   turtle('INIT');

   drawkoch(length, level);
   turtle('GIRA', -120);
   drawkoch(length, level);
   turtle('GIRA', -120);
   drawkoch(length, level);
end

El resultado de la función para los primeros niveles se puede ver en la siguiente figura:

La curva C de Lévy¶

El último ejemplo que vamos a mostrar es el de la curva C de Lévy, descrita por primera vez por Ernesto Cesàro en 1906 y Georg Faber en 1910, pero que lleva el nombre del matemático francés Paul Lévy, que fue el primero que proporcionó un método geométrico para construirla.

La curva C comienza con un segmento recto de cierta longitud. Utilizando este segmento como hipotenusa, se construye un triángulo isósceles con ángulos 45º, 90º y 45º. Sustituimos el segmento original por los dos catetos del triángulo. El proceso se repite indefinidamente con cada uno de los segmentos que se generan. En la figura siguiente se muestran los primeros niveles de la curva C:

El código que se muestra a continuación dibuja la curva C para un determinado nivel de profundidad:

function ccurve(length,level)
   if level == 0
      turtle('AVANZA',length);
   else
      turtle('GIRA',45);
      ccurve(length/sqrt(2), level-1);
      turtle('GIRA',-90);
      ccurve(length/sqrt(2), level-1);
      turtle('GIRA',45);
   end
end

El resultado para el nivel de profundidad 12 se muestra en la siguiente figura:

Introducción¶

Actualmente los Modelos Digitales de Elevaciones, MDE, sustituyen a los antiguos mapas de curvas de nivel. Consiste en dividir el terreno mediante una cuadrícula rectángular y asignar una cota a cada rectángulo resultante de la división. Esto permite representar el MDE mediante una matriz. Un formato habitual para este tipo de mapas es el formato ESRI ASC. Consiste en un fichero de texto en el que se guardan los valores de las alturas en formato coma flotante junto con unas lineas al principio que indican el número de filas y columnas del MDE, las coordenadas de la esquina inferior izquierda del mapa, (Lower Left Corner), el tamaño en metros del lado de los rectángulos de la rejilla, (paso de malla), y a veces algún dato más.

La figura siguiente muestra la cabecera de uno de estos ficheros .asc según se ve al abrirlo con un editor de texto y al abrirlo con el programa gvSIG:

Actualmente es posible descargar modelos digitales de elevaciones a nivel mundial con paso de malla 30 metros, desde el portal de la NASA http://asterweb.jpl.nasa.gov/gdem.asp. En España está disponible el MDE con paso de malla 5 metros. Se puede descargar desde el portal del Instituto Geográfico Nacional http://centrodedescargas.cnig.es/CentroDescargas.

Problema¶

Se pide desarrollar una rutina en Octave-Matlab que lea el fichero mdt25.asc y realice una representación tridimensional del mismo.

Solución¶

Si abrimos el fichero mdt25.asc con un editor de texto, veremos que tiene seis líneas de metadatos al principio, y luego vienen las verdaderas líneas de datos. Los metadatos del fichero mdt25.asc son los siguientes:

• NCOLS, NROWS: Número de filaas y columnas de la rejilla rectángular en que se ha dividido el terreno, en este caso 601x410
• XLLCORNER, YLLCORNER: Coordenadas UTM de la esquina inferior izquierda del MDE
• CELLSIZE: Tamaño del lado de la celda, en este caso cuadrada de lado 25 metros
• NODATA_VALUE: Las celdas que tengan este valor quire decir que no hay dato de altura de la celda.

En la solución que vamos a desarrollar aquí no haremos uso de los metadatos, salvo para saber el número de filas y columnas de la matriz, dato que miraremos abriendo el fichero con un editor de texto. La solución propuesta abre el fichero mdt25.asc, se salta las primeras seis líneas del fichero utilizando fgets() y a continuación lee los valores de las alturas en una matriz. Dada la forma de trabajar de la función fscanf(), que va leyendo valores y rellenando por columnas, es necesario leer la matriz con tamaño [columnasxfilas] y luego trasponerla. A continuación se muestra el código de la rutina Octave, y de la figura resultante:

function readmdt(filename)
% Lee un fichero raster en formato ASC

   file = fopen(filename,'r');
   if file <3
      fprintf('Error al abrir fichero %s', filename)
      return;
   end

   for i=1:6
      metadato = fgets(file);
      fprintf("%s", metadato)
   end

   filas = 410;
   columnas = 601;
   A = fscanf(file, '%f', [columnas, filas]);
   A = A';

   hold off;
   close;
   hold on;
   xlim([0,610]);
   ylim([0,410]);
   xlabel('x');
   ylabel('y');
   view(240,60);

   meshz(A);

   fclose(file);

end