Matemáticas con Mathematica

Práctica 5. Elementos de programación con Mathematica

Francisco Javier Pérez González
Departamento de Análisis Matemático
Universidad de Granada

Introducción

Las expresiones que escribes para que un ordenador realice una determinada tarea consituyen un programa y reciben el nombre de código fuente. Lo primero que necesitas para escribir un programa es un lenguaje de programación que te permita comunicarte con el ordenador. Hay muchos lenguajes de programación con diferencias notables entre ellos. Mathematica tiene su propio lenguaje de programación: es el que usas cuando escribes comandos y expresiones que pueden ser ejecutados en una celda de input. Los ordenadores no entienden el código fuente escrito en ninguno de los lenguajes de programación usuales, dicho código debe ser traducido a código máquina que es el único que entiende un microprocesador como el que tiene el ordenador que estás usando. 
Los lenguajes de programación se dividen en dos grandes categorías dependiendo de la forma en que se haga la conversión de código fuente a código máquina. Cuando ésta debe hacerse (una sola vez), antes de ejecutar el programa en el ordenador, por medio de un compilador que realiza la conversión y produce un archivo que contiene el código máquina el cual podrá ser ejecutado posteriormente sin necesidad de nueva conversión, se dice que se trata de un lenguaje compilado. Lenguajes como C, Fortran o Pascal son lenguajes compilados. Cuando la conversión debe hacerse sobre la marcha al mismo tiempo que se ejecuta el programa (y cada vez que se ejecuta el programa), se dice que se trata de un lenguaje interpretado. Puesto que la conversión requiere un tiempo, los lenguajes compilados son mucho más rápidos que los interpretados. Mathematica es un lenguaje interpretado aunque dispone de recursos que permiten "compilar expresiones" y traducirlas internamente a un lenguaje próximo al código máquina. Si lo piensas un poco te convencerás de que un programa como Mathematica tiene que ser así: ¿te imaginas que tuvieras que compilar el código de cada celda de input antes de ejecutarla?
Un mismo programa puede implementarse de muy distintas formas. En este sentido suelen distinguirse varios estilos de programación (paradigmas de programación) que se diferencian por la forma de organizar el código fuente y en la preferencia por cierto tipo de acciones o de datos. Una virtud de Mathematica es que permite programar según diferentes principios y paradigmas. 
En esta práctica vamos a considerar algunos aspectos propios de un estilo de programación, llamado programación procedimental, que usa con frecuencia estructuras condicionales ("If[]", "Which[]", "While[]") que permiten evaluar ciertas expresiones dependiendo del resultado de un test, así como estructuras recursivas y bucles ("Do[]", "For[]"). Pascal, Fortran, C son lenguajes de programación procedimental.

Reasignación de variables

En las prácticas anteriores has aprendido a usar el comando de asignación "=" para declarar e inicializar variables. Es frecuente que queramos cambiar el valor de una variable por medio de una operación aritmética y asignarle el nuevo valor a dicha variable. Podemos hacer esto usando el comando "=" de la misma forma que antes lo usamos para declarar la variable pero no debes olvidar que el lado de la derecha se evalúa antes que la nueva asignación sea hecha. Considera la siguiente celda:

var = 10;
var = var + 10

(Observa que el ; suprime la salida de la primera línea). Lo que ocurre en la segunda línea es que "var+10" es evalúado al entero 20 el cual es asignado a "var".  Desde un punto de vista matemático puede parecer extraño el uso del igual en la expresión "var=var+10" pero no olvides que "=" en Mathematica no significa "igual" sino que es un comando de asignación por el cual se atribuye un valor a una expresión (recuerda que en Mathematica el doble signo "==" es el test de igualdad).

Aquí tienes otros ejemplos:

var = var * 2

var = var / 5

En los bucles de programación suelen usarse iteradores que no son sino variables que se inicializan en un valor y se van incrementando o decrementando hasta que alcanzan un cierto valor en cuyo momento el bucle se detiene.  Una expresión cómoda para incrementar una variable en una unidad es "++". Por ejemplo "i++" indica que se aumente el valor actual de i en una unidad. Análogo significado tiene "i--". Debes tener en cuenta que el incremento (o decremento) tendrá lugar después de que el valor actual de i se haya usado en la expresión correspondiente.

var = 9var++

Aunque obtenemos 9 como valor de var su valor actual es 10.

var

Lo que ha sucedido es que al ejecutar "var++" el valor de esta expresión era el de var, que era 9, y después que se hizo la evaluación (que devolvió el valor 9), var fue incrementada.

Variables locales y variables globales

Cada vez que declaras una variable o una función Mathematica la guarda en algún contexto. Por defecto el contexto en el que se guardan todas las variables y funciones definidas por el usuario es el contexto Global. Esto es así para que dichas variables y funciones puedan ser usadas en cualquier celda en cualquier momento.

?var

Mathematica te dice que var es una variable almacenada en el contexto Global cuyo valor actual es 10. Muchas veces es conveniente usar alguna variable o alguna función para realizar una tarea muy específica y no nos interesa que dicha variable o función pueda ser accesible desde otras instancias del programa. Para ello Mathematica dispone del comando Module[{x,y,... },expr]que especifica que los símbolos x, y,... en expr deben ser tratados como variables locales.

Module[{var}, var = 0; var](* Mathematica nos da el valor de la variable local "var" *)

var

Aquí tienes otro ejemplo para que entiendas la diferencia entre variable global y local.

\(t = 5; 
    \  (*\ declaramos\ "t"\ como\ variable\ global\ con\ valor\ inicial\ 5
        \ *) 
f[x_] := Module[{t}, t = \((1 + x)\)^2; Expand[t]] (*\ 
      declaramos "t" como\ una\ variable\ local\ *) \)f[a + b]

?t

Clear["@"]

Evaluación condicional

If

El uso más frecuente del comando "If[]" es con tres argumentos.

If[test,expr1,expr2]

El primer argumento debe ser una expresión boleana. Si test=True,  expr1 es ejecutada, mientras que si test=False, expr2 es ejecutada. Las expresiones expr1 y expr2 pueden estar formadas por varias instrucciones separadas por punto y coma.

Por ejemplo, podemos usar "If[]"  para definir una función que devuelva 1 si su argumento es positivo y devuelva -1 si es negativo.

Clear[f]
f[x_] := If[x>0,1,-1] (* observa que es obligado el uso de := si no lo usamos Mathematica trataría de evaluar inmediatamente la condición x>0 lo que daría un error *)
f[-.2]
f[.01]

Podemos modificar la función para que devuelva 0 si su argumento es igual a 0. Para ello podemos anidar dos "If[]" como sigue.

Clear[f]
f[x_]:= If[x>0,1,If[x==0,0,-1]]

Si el argumento es mayor que 0, "f[]" devuelve un 1; si no, si es igual a cero, devuelve un 0; en otro caso devuelve un -1.

Which

Este comando se usa cuando queremos elegir entre varias acciones posibles dependiendo del resultado de varios tests. En la función antes definida hemos necesitado dos tests: uno para comprobar si el argumento es mayor que 0 y otro para comprobar si es igual a 0. En estos casos el comando "Which[]" es apropiado porque permite hacer todo el trabajo en una estructura sintáctica sencilla. El comando "Which[]" considera sus argumentos por pares formados por un test y una expresión. La sintaxis es

Which[test1,expr1,test2,expr2,...]

y devuelve el valor de la primera expresión para la que el correspondiente test da True. Las expresiones expr1,expr2,... pueden estar formadas por varias instrucciones separadas por punto y coma. Podríamos haber definido nuestra función anterior como sigue:

Clear[f]
f[x_]:= Which[x>0,1,x==0,0,x<0,-1]

Podemos usar un "Which[]" para definir una función que calcule el factorial de un número.

Clear[fac]
fac[n_]:= Which[n>1,n*fac[n-1],n==1,1]

For

El comando For se usa para realizar bucles. Un “bucle” es un proceso repetitivo que se realiza un cierto número de veces. La sintaxis de la orden es

For[expresióninicial,test,incr,expr]

• expresióninicial nos sitúa en las condiciones de comienzo del bucle.
• test dirá a Mathematica el momento de detener el proceso.
• incr expresará la forma de modificar la condición inicial.
• expr dirá a Mathematica lo que tiene que realizar en cada paso; la expresión puede estar compuesta de varias instrucciones separadas mediante punto y coma.
El orden de evaluación es test->expr->incr.El proceso termina cuando test da False. Un aspecto importante de For es que su resultado es Null es decir, que no proporciona resultado en pantalla salvo que tú se lo indiques de alguna forma explícita. Vamos a ver algunos ejemplos del funcionamiento de este comando. Ejecuta la siguiente celda.

Clear[i]For[i = 1, i < 500, i++, i = i^2]

Aparentemente no ha pasado nada pero eso no es así. Comprueba el valor actual de i (trata de averiguarlo antes de preguntárselo a Mathematica).

i

¿No te lo esperabas? Bien, siempre que quieras saber con total claridad todas las operaciones que ha realizado Mathematica para ejecutar un comando tienes a tu disposición la orden Trace (un impagable regalo de Mathematica). Observa su funcionamiento.

Trace[For[i = 1, i < 500, i++, i = i^2]]

Demasiado largo para explicarlo con palabras pero no sigas adelante antes de entender claramente el significado de todos y cada uno de los pasos dados por Mathematica. Fíjate especialmente en la forma en que se va incrementando el valor de i después de realizar la operación indicada y por qué el valor actual de i es 677. Una forma de lograr que For proporcione una salida en pantalla  es usando el comando Print[expr]que, precisamente, sirve para imprimir en pantalla.

For[i = 1, i < 500, i++, i = i^2; Print[i]]

El último valor de i que Mathematica ha impreso en pantalla es 676 = (((25 + 1))^2) porque ese fue el último valor calculado de i^2. Sin embargo el valor actual de i es 677 porque después de hacer i=i^2=676 ha incrementado el valor en una unidad al ejecutar la orden i++. La orden Print debes usarla con prudencia.

El comando For también puede usarse en la forma For[expresióninicial,test,incr]. Aquí tienes un ejemplo de ese uso de For con el que podemos comprobar que la sucesión (((1 + 1/n))^n) converge a E aunque muy lentamente. Lo que hacemos es calcular el primer número n tal que  | (((1 + 1/n))^n) - E |  < (10^(-4)).

For[n = 1, Abs[(1.  + 1/n)^n - E] ≥ 10^-4, 
      n = n + 1]  nAbs[(1.  + 1/n)^n - E]

La sucesión (∑_(k = 0)^n) (1/(k!)) converge a E mucho más rápidamente.

For[n = 1, 
      Abs[(∑_(k = 0)^n) (1. /(k!)) - E] ≥ 10^-4, 
      n = n + 1]nAbs[(∑_(k = 0)^n) (1. /(k!)) - E]

While

While[test,expr] permite realizar repetidamente las instrucciones indicadas en expr mientras test da True y se detiene cuando Test da False. Aquí tienes un ejemplo del uso de While para aproximar, usando el método de Newton-Raphson, una solución de la ecuación (x^3) + 1 = 0 partiendo de un valor inicial b.

\(f[b_] := 
      Module[{x, k = 0}, 
        x = b;  (*\ 
          Usando\ Module\ las\ variables\ quedan\ localizadas . \ Las\ 
              variables\ locales\ son\ x\ y\ k . \ La\ variable\ k\ se\ ha\ 
              inicializado\ con\ k = 0\ *) 
		
        While[Abs[x\^3 + 1] > 0.001\  && \ k ≤ 50\  && \ 
            x ≠ 0, x = \(\(-1\) + 2  x\^3\)\/\(3\ x\^2\); 
          k = k + 1]; 

        		{k, x}\  (*\ 
          pedimos\ que\ nos\ devuelva\ el\ útimo\ valor\ de\ k\ que\ 
              es\ el\ número\ de\ iteraciones\ realizadas\  - 
            observa\ que\ k\ empieza\ en\ 0 - 
            y\ de\ x\ que\ es\ el\ valor\ aproximado\ de\ la\ raíz\ 
              obtenido\ *) ]\)

Prueba distintos valores para b no demasiado lejanos de la raíz x = -1.

f[- .1]f[ .15]f[ .55]

También podemos hacer lo mismo con un For.

Clear[f]\(f[b_] := 
      Module[{x, k}, 
		
        For[x = b; k = 0, 
          Abs[x\^3 + 1] > 0.001\  && \ k ≤ 50\  && \ 
            x ≠ 0, \(k++\), x = \(\(-1\) + 2  x\^3\)\/\(3\ x\^2\)]; 
        
		{k, x}]\)

f[- .1]f[ .15]f[ .55]

Break

Este comando, cuya sintaxis es simplemente Break[], tiene como única finalidad interrumpir la ejecución de un bucle. Esto puede ser interesante si, por ejemplo, creemos que el resultado que buscamos puede alcanzarse antes de que el test del bucle sea False. En el siguiente ejemplo se puede comprender el uso de esta orden. Recuerda que en la práctica primera comprobamos que es posible sumar una cantidad muy pequeña al número 1. de forma que el resultado de la suma, debido a los errores de redondeo de la aritmética de punto flotante, seguía siendo igual a 1. Podemos calcular el orden de magnitud de esa cantidad tan pequeña con el siguiente programa.

\(For[i = 1, i ≤ 100, \(i++\), 
      
If[1.  + 10\^\(-i\) == 1. , 
        
\ \ 
        Print[\*""El exponente i con el que se logra la igualdad 1.+ (10^(-i))== 1. es i="", i]; 
        
Break[],  (*\ 
          para\ interrumpir\ la\ ejecución\ del\ bucle\ cuando\ se\ 
            dé\ la\ igualdad\ *) 
Print[
          "Para ", 10. \^\(-i\), " no se da la igualdad porque ", 
          N[1.  + 10\^\(-i\), 14], " es distinto de ", \ \ 1. ]
          
]
]\)

Aquí tienes un uso de While para definir una función que detecta el primer elemento negativo de una lista.

primerNeg[x_] := Module[{i}, 
(* Usamos Module para localizar la variable i *)
   i = 1;
   While[ x[[i]]>=0, i=i+1 ];
   x[[i]]
]

primerNeg[{3,5,2,-7,4,3,11,-2,0,-1,8}]

La función anterior no está del todo bien definida. Compruébalo.

primerNeg[{3,5,2,7,4,3,11,2,0,1,8}]

Lo que ha ocurrido es que, al no haber ningún elemento negativo en la lista, el iterador ha alcanzado el valor 12 y la función ha tratado de comprobar si el elemento que ocupa dicho lugar en la lista es ≥0 pero la lista solamente tiene 11 elementos.

Length[{3, 5, 2, 7, 4, 3, 11, 2, 0, 1, 8}]

Uno de los ejercicios que te propongo es que modifiques esta función para evitar ese error.

Principios generales para trabajar con bucles iterativos

1. Debemos ininializar una o más variables (iteradores) antes de empezar el bucle.  Es muy conveniente usar Module para localizar dichas variables.

2. Modificamos los iteradores cada vez que se realiza el bucle. Lo más frecuente es incrementar o decrementar las iteradores.

3. Después de que el bucle haya terminado debemos ejecutar algún comando para pedir una respuesta. Este comando puede consistir simplemente en pedir los últimos valores de los iteradores o de las operaciones realizadas. Por ellos mismos, los bucles no devuelven ningún valor (salida Null). En el ejemplo anterior hemos puesto "x[[i]]" para que el bucle nos devuelva el primer elemento negativo de la lista.

4. Debes asegurarte de que el bucle terminará en un número razonable de pasos y evitar entrar en un bucle infinito.

Un ejemplo más

La función "PrimeQ[]" es un test boleano que comprueba si un número es primo.

?PrimeQ

Vamos a usar dicha función para definir otra que nos devuelva el primer número primo mayor que un número dado.

primerPrimoDespues[n_] := 
   Module[{i},
      i = n+1;
      While[ Not[PrimeQ[i] ], i++];
      i
   ]

primerPrimoDespues[2004]

Dos errores que debes evitar.

Un bucle no devuelve ningún valor si no lo pides

Clear[primerPrimoDespues]
primerPrimoDespues[n_] := 
   Module[{i},
      i = n+1;
      While[ Not[PrimeQ[i] ], i++]
   ] (* esta función hace correctamente su trabajo pero no devuelve nada *)

La forma de arreglar esto es pedirle a Mathematica que te proporcione el resultado deseado.

Clear[primerPrimoDespues]
primerPrimoDespues[n_] := 
   Module[{i},
      i = n+1;
      While[ Not[PrimeQ[i] ], i++]; i (* esta es la pequeña diferencia *)
   ]

Mathematica no permite que en un bucle se usen parámetros formales

Clear[primerPrimoDespues]
primerPrimoDespues[n_] := 
   Module[{},
      n = n+1;
      While[ Not[PrimeQ[n] ], n++];
      n
   ] (* esta función da un error si tratas de usarla *)

La forma de arreglar esto es fácil: basta localizar las variables dentro de "Module[ ]".

Clear[primerPrimoDespues]
primerPrimoDespues[m_] := 
   Module[{n=m},(* esta es la pequeña diferencia *)
      n = n+1;
      While[ Not[PrimeQ[n] ], n++];
      n
   ]

Do

Se trata de un comando que permite hacer repetidamente una operación. Su sintaxis es muy sencilla Do[expr,{iter}]donde expr es una o varias instrucciones separadas por punto y coma y {iter} puede ser de varias formas:
{iter}={número} Repite las operaciones indicadas tantas veces como induca "número".
{iter}={i,número} Repite las operaciones indicadas desde i=1 hasta i=número incrementando i de uno en uno.
{iter}={i,número1,numero2} Repite las operaciones indicadas desde i=numero1 hasta i=número2 incrementando i de uno en uno.
{iter}={i,número1,numero2,incr} Repite las operaciones indicadas desde i=numero1 hasta i=número2 incrementando i de incr en incr.

Aquí tienes un ejemplo del uso de Do para sumar los cuadrados de los números pares desde 2000 a 3000. Observa que debemos inicializar una variable "sum=0" en la que iremos almacenando los cálculos y al final le pedimos que nos devuelva el valor de dicha variable.

Module[ {sum,i},
   sum = 0;
   Do[sum = sum + i^2, {i,2000,3000,2}];
   sum
]

Programación basada en listas

Mathematica proporciona la posibilidad de programar usando listas. Con frecuencia es posible sustituir procesos iterativos y bucles por el uso apropiado de listas. Supongamos, por ejemplo, que queremos programar el cálculo de los puntos medios de una lista de números; es decir, se trata de calcular los puntos medios de cada dos números consecutivos de la lista dada. Hagamos esto de dos formas disferentes. En la primera usamos programación procedimental y en la segunda programación basada en listas.

\(puntosmedios[lista_] := 
      Module[{x},  (*\ localizamos\ la\ variable\ "x"\ *) \ 
			
        x = Table[0, {Length[lista] - 1}]\  (*\ valor\ inicial\ de\ x\ *) ; 

        		For[k = 1, k ≤ \ Length[lista] - 1, 
          x[\([k]\)] = \((lista[\([k]\)] + lista[\([k + 1]\)])\)*0.5; 
          \(k++\)] (*\ calculamos\ los\ puntos\ medios\ *) \ ; 
		x
        \  (*\ pedimos\ a\ Mathematica\ que\ devuelva\ el\ valor\ actual\ de\ 
            x\ *) ]\)

puntosmedios[Range[0, 100]]

nuevospuntosmedios[lista_] := 
      (Drop[lista, 1] + Drop[lista, -1])*0.5

nuevospuntosmedios[Range[0, 100]]

Pidamos a Mathematica que nos diga el tiempo que tarda en calcular con los comandos que acabamos de definir los puntos medios de una lista con los primeros cien mil números.

Timing[puntosmedios[Range[0, 100000]]][[1]]

Timing[nuevospuntosmedios[Range[0, 100000]]][[1]]

Como puedes comprobar la segunda forma es mucho más eficaz que la primera.

Ejercicios

1 Define una función que devuelva el primer elemento negativo de una lista y en caso de que la lista no tenga elementos negativos devuelva la expresión "la lista no tiene elementos negativos".

2 Dado un número x>0, se puede consguir que la suma  1/1 + 1/2 + 1/3 + 1/4 + ...+1/n  sea mayor que x tomando un número n suficientemente grande. Define una función que calcule el primer número n que satisface dicha condición. Calcula n para x=10,11,13. 
Sugerencia. Usa números reales en tu función para que Mathematica no calcule simbólicamente.

3 Define una función que calcule el menor número natural n tal que (2^n) es mayor que un número dado x.

4 Los números de Fibonacci son los elementos de la sucesión 1,1,2,3,5,8,13,21,34,... cuyos primeros dos elementos son iguales a 1 y cada elemento posterior es la suma de los dos que le preceden. Define una función que calcule el n-ésimo número de Fibonacci F[n]. Utiliza ListPlot[ ] (consulta la ayuda de Mathematica) para representar los números de Fibonacci desde F[200] hasta F[250]. Represéntalos en rojo con un tamaño apropiado.