Source: http://docplayer.es/12709-Primeros-pasos-en-maxima.html
Timestamp: 2017-09-25 22:42:21+00:00

Document:
Primeros pasos en Maxima - PDF
Download "Primeros pasos en Maxima"
Gloria Saavedra Roldán
1 Primeros pasos en Maxima Mario Rodríguez Riotorto 13 de febrero de 015
2 Copyright c Mario Rodríguez Riotorto Este documento es libre; se puede redistribuir y/o modificar bajo los términos de la GNU General Public License tal como lo publica la Free Software Foundation. Para más detalles véase la GNU General Public License en This document is free; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation. See the GNU General Public License for more details at
4 Índice general 1. Introducción 3. Despegando con Maxima 7.1. Instalación Una primera sesión Números y sus operaciones Enteros y decimales Números complejos Combinatoria Estructuras de datos Listas Arrays Cadenas Conjuntos Grafos Álgebra Transformaciones simbólicas Ecuaciones Inecuaciones Matrices Patrones y reglas Cálculo Funciones matemáticas Límites Derivadas Integrales Ecuaciones diferenciales Vectores y campos vectoriales
5 ÍNDICE GENERAL 7. Gráficos El módulo plot El módulo draw Probabilidades y estadística Probabilidad Estadística descriptiva Estadística inferencial Otros métodos numéricos Interpolación Ajuste por mínimos cuadrados Optimización con restricciones Programación lineal Programando en Maxima Programación a nivel de Maxima Programación a nivel de Lisp
6 Capítulo 1 Introducción Este es un manual introductorio al entorno de cálculo simbólico Maxima, directo sucesor del legendario Macsyma. El objetivo del manual es facilitar el acceso a este programa a todas aquellas personas que por vez primera se acercan a él. Maxima es un programa cuyo objeto es la realización de cálculos matemáticos, tanto simbólicos como numéricos; es capaz de manipular expresiones algebraicas y matriciales, derivar e integrar funciones, realizar diversos tipos de gráficos, etc. Su nombre original fue Macsyma (MAC s SYmbolic MAnipulation System, donde MAC, Machine Aided Cognition, era el nombre del Laboratory for Computer Science del MIT durante la fase inicial del proyecto Macsyma). Se desarrolló en estos laboratorios a partir del año 1969 con fondos aportados por varias agencias gubernamentales norteamericanas (National Aeronautics and Space Administration, Office of Naval Research, U.S. Department of Energy y U.S. Air Force). El concepto y la organización interna del programa están basados en la tesis doctoral que Joel Moses elaboró en el MIT sobre integración simbólica. Según Marvin Minsky, director de esta tesis, Macsyma pretendía automatizar las manipulaciones simbólicas que realizaban los matemáticos, a fin de entender la capacidad de los ordenadores para actuar de forma inteligente. El año 198 es clave. El MIT transfiere una copia de Macsyma a la empresa Symbolics Inc. para su explotación económica, haciendo el código propietario, y otra al Departamento de Energía, copia ésta que será conocida con el nombre de DOE-Macsyma. En 199 la versión comercial de Macsyma sería adquirida por una empresa que se llamaría precisamente Macsyma Inc, y el programa iría perdiendo fuelle progresivamente ante la presencia en el mercado de otros programas similares como Maple o Mathematica, ambos los dos inspirados en sus orígenes en el propio Macsyma. Pero ocurrieron dos historias paralelas. Desde el año 198, y hasta su fallecimiento en el 001, William Schelter de la Universidad de Texas man- 3
7 4 CAPÍTULO 1. INTRODUCCIÓN tuvo una versión de este programa adaptada al estándar Common Lisp en base a DOE-Macsyma, la cual ya se conoce con el nombre de Maxima para diferenciarla de la versión comercial. En el año 1998 Schelter había conseguido del DOE permiso para distribuir Maxima bajo la licencia GNU-GPL (www.gnu.org); con este paso, muchas más personas empezaron a dirigir su mirada hacia Maxima, justo en el momento en el que la versión comercial estaba declinando. Actualmente, el proyecto está siendo mantenido por un grupo de desarrolladores originarios de varios países, asistidos y ayudados por otras muchas personas interesadas en Maxima y que mantienen un cauce de comunicación a través de la lista de correo Además, desde el año 006, se ha activado una lista de correos para usuarios de habla hispana en Puesto que Maxima se distribuye bajo la licencia GNU-GPL, tanto el código fuente como los manuales son de libre acceso a través de la página web del proyecto, Uno de los aspectos más relevantes de este programa es su naturaleza libre; la licencia GPL en la que se distribuye brinda al usuario ciertas libertades: libertad para utilizarlo, libertad para modificarlo y adaptarlo a sus propias necesidades, libertad para distribuirlo, libertad para estudiarlo y aprender su funcionamiento. La gratuidad del programa, junto con las libertades recién mencionadas, hacen de Maxima una formidable herramienta pedagógica, de investigación y de cálculo técnico, accesible a todos los presupuestos, tanto institucionales como individuales. No quiero terminar esta breve introducción sin antes recordar y agradecer a todo el equipo humano que ha estado y está detrás de este proyecto, desde el comienzo de su andadura allá por los años 70 hasta el momento actual; incluyo también a todas las personas que, aún no siendo desarrolladoras del programa, colaboran haciendo uso de él, presentando y aportando continuamente propuestas e ideas a través de las listas de correos. Un agradecimiento especial a Robert Dodier y James Amundson, quienes me brindaron la oportunidad de formar parte del equipo de desarrollo de Maxima y me abrieron las puertas a una experiencia verdaderamente enriquecedora.
8 5 Finalmente, un recuerdo a William Schelter. Gracias a él disponemos hoy de un Maxima que, paradójicamente, por ser libre, no pertenece a nadie pero que es de todos. Ferrol-A Coruña
9 6 CAPÍTULO 1. INTRODUCCIÓN
10 Capítulo Despegando con Maxima.1. Instalación Maxima funciona en Windows, Mac y Linux. En Windows, la instalación consiste en descargar el binario exe desde el enlace correspondiente en la página del proyecto y ejecutarlo. En Mac, es necesario descargar el fichero con extensión dmg, hacer doble clic sobre él y arrastrar los iconos de Maxima, wxmaxima y Gnuplot a la carpeta Aplicaciones. El proceso detallado se explica aquí: En Linux, la mayoría de las distribuciones tienen ficheros precompilados en sus respectivos repositorios con las extensiones rpm o deb, según el caso. Sin embargo, en algunas distribuciones las versiones oficiales no suelen estar muy actualizadas, por lo que podría ser recomendable proceder con la compilación de las fuentes, tal como se expone a continuación. Las siguientes indicaciones sobre la compilación hacen referencia al sistema operativo Linux. Si se quiere instalar Maxima en su estado actual de desarrollo, será necesario descargar los ficheros fuente completos desde el repositorio de Sourceforge y proceder posteriormente a su compilación. Se deberá tener operativo un entorno Common Lisp en la máquina (clisp, cmucl, sbcl o gcl son todas ellas alternativas válidas), así como todos los programas que permitan ejecutar las instrucciones que se indican a continuación, incluidas las dependencias tcl-tk y gnuplot. Grosso modo, estos son los pasos a seguir para tener Maxima operativo en el sistema 1 : 1. Descargar las fuentes a un directorio local siguiendo las instrucciones que se indican en el enlace al repositorio de la página del proyecto. 1 Se supone que se trabaja con clisp, con codificación unicode y que queremos el sistema de ayuda en inglés y español. 7
11 8 CAPÍTULO. DESPEGANDO CON MAXIMA. Acceder a la carpeta local de nombre maxima y ejecutar las instrucciones./bootstrap./configure --enable-clisp --enable-lang-es-utf8 make make check sudo make install make clean 3. Para ejecutar Maxima desde la línea de comandos (Figura.1) se deberá escribir maxima si se quiere trabajar con la interfaz gráfica (Figura.), teclear xmaxima Otro entorno gráfico es Wxmaxima (Figura.3), que se distribuye conjuntamente con los ejecutables para Windows y Mac. En Linux, una vez instalado Maxima, se podrá instalar este entorno separadamente. Una cuarta alternativa, que resulta muy vistosa por hacer uso de TEX, es la de ejecutar Maxima desde el editor Texmacs 3 (Figura.4). Emacs es un editor clásico en sistemas Unix. Aunque hay varias formas de ejecutar Maxima desde Emacs, con el modo Imaxima 4 (Figura.5) se consigue que los resultados estén formateados en TEX. También se han desarrollado en los últimos años entornos web para Maxima 5... Una primera sesión En esta sesión se intenta realizar un primer acercamiento al programa a fin de familiarizarse con el estilo operativo de Maxima, dejando sus habilidades matemáticas para más adelante. Una vez accedemos a Maxima, lo primero que vamos a ver será la siguiente información: Maxima using Lisp CLISP.48 ( )
12 .. UNA PRIMERA SESIÓN 9 Figura.1: Maxima desde la línea de comandos. Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) En primer lugar nos informa sobre la versión de Maxima que se está ejecutando y la URL del proyecto, a continuación indica qué entorno Lisp está dándole soporte al programa y la nota legal sobre la licencia, por supuesto GNU-GPL. Finaliza con un recordatorio a Schelter y cómo debemos proceder si encontramos un fallo en el programa, ayudando así a su mejora. Tras esta cabecera, muy importante, aparece el indicador (%i1) esperando nuestra primera pregunta (i de input). Si queremos calcular una simple suma tecleamos la operación deseada seguida de un punto y coma (;) y una
13 10 CAPÍTULO. DESPEGANDO CON MAXIMA Figura.: Xmaxima, el entorno gráfico basado en tcl-tk. pulsación de la tecla retorno (%i1) ; ( %o1) 68 (%i) Vemos que el resultado de la suma está etiquetado con el símbolo (%o1) indicando que es la salida correspondiente a la primera entrada (o de output). A continuación (%i) indica que Maxima espera nuestra segunda instrucción. El punto y coma actúa también como un separador cuando escribimos varias instrucciones en un mismo renglón. Nuestra siguiente operación consistirá en asignar el valor a la variable x y el a la y, solicitando luego su producto,
14 .. UNA PRIMERA SESIÓN 11 Figura.3: Wxmaxima, el entorno gráfico basado en wxwidgets. (%i) x:34578; y:984003; x*y; ( %o) ( %o3) ( %o4) conviene prestar atención al hecho de que la asignación de un valor a una variable se hace con los dos puntos (:), no con el signo de igualdad, que se reserva para las ecuaciones. Es posible que no deseemos los resultados intermedios que Maxima va calculando o, como en este caso, las asignaciones a las variables que va
15 1 CAPÍTULO. DESPEGANDO CON MAXIMA Figura.4: Maxima desde Texmacs. haciendo; en tales situaciones debe hacerse uso del delimitador $, que no devuelve los resultados que va calculando. Repitiendo el cálculo de la forma (%i5) x:34578$ y:984003$ x*y; ( %o7) podemos conseguir una salida más limpia. Las asignaciones a variables se mantienen activas mientras dure la sesión con Maxima, por lo que podemos restar las variables x e y anteriores (%i8) x-y; ( %o13) Esto tiene un riesgo; si queremos resolver la ecuación x 3x + 1 = 0,
16 .. UNA PRIMERA SESIÓN 13 Figura.5: Maxima en Emacs. (%i9) solve(x^-3*x+1=0,x); A number was found where a variable was expected - solve -- an error. Quitting. To debug this try debugmode(true); nos devuelve un mensaje de error, ya que donde se supone que hay una incógnita, lo que realmente encuentra es un número, en este caso Hay dos formas de resolver esta situación; la primera consiste en utilizar el operador comilla simple ( ), que evita la evaluación de las variables: (%i9) solve( x^-3* x+1=0, x); ( %o9) [ 5 3 x =, x = ] 5 + 3
17 14 CAPÍTULO. DESPEGANDO CON MAXIMA (%i10) x; ( %o10) como se ve, x mantiene su valor original. La segunda forma de evitar este problema consiste en vaciar el contenido de la variable x mediante la función kill, (%i11) kill(x)$ (%i1) solve(x^-3*x+1=0,x); ( %o1) (%i13) x; [ 5 3 x =, x = ] ( %o13) x en este caso, x ha perdido su valor. Las etiquetas (%in) y (%on), siendo n un número entero, pueden ser utilizadas a lo largo de una sesión a fin de evitar tener que volver a escribir xy expresiones; así, si queremos calcular el cociente x y, con x = y y = , podemos aprovechar los resultados (%o7) y (%o8) y hacer (%i14) %o7/%o8; ( %o14) Las etiquetas que indican las entradas y salidas pueden ser modificadas a voluntad haciendo uso de las variables inchar y outchar, (%i15) inchar; ( %o15) %i (%i16) outchar; ( %o16) %o (%i17) inchar: menda; ( %o17) menda (menda18) outchar: lerenda;
18 .. UNA PRIMERA SESIÓN 15 (lerenda18) lerenda (menda19) +6; (lerenda19) 8 (menda0) inchar: %i$ outchar: %o$ (%i) Al final se han restaurado los valores por defecto. En caso de que se pretenda realizar nuevamente un cálculo ya indicado deberán escribirse dos comillas simples ( ), no comilla doble, junto con la etiqueta de la instrucción deseada (%i3) %i8; ( %o3) x Obsérvese que al haber dejado a x sin valor numérico, ahora el resultado es otro. Cuando se quiera hacer referencia al último resultado calculado por Maxima puede ser más sencillo hacer uso del símbolo %, de forma que si queremos restarle la x a la última expresión podemos hacer (%i4) %-x; ( %o4) Ciertas constantes matemáticas de uso frecuente tienen un símbolo especial en Maxima: la base de los logaritmos naturales (e), el cociente entre la longitud de una circunferencia y su diámetro (π), la unidad imaginaria (i = 1), la constante de Euler-Mascheroni (γ = 0 e x ln xdx) y la razón aurea (φ = 1+ 5 ), se representarán por %e, %pi, %i, %gamma y %phi, respectivamente. Es muy sencillo solicitar a Maxima que nos informe sobre alguna función de su lenguaje; la técnica nos servirá para ampliar información sobre cualquier instrucción a la que se haga referencia en este manual, por ejemplo, (%i5) describe(sqrt); -- Función: sqrt (<x>) Raíz cuadrada de <x>. Se representa internamente por <x>^(1/). Véase también rootscontract. Si la variable radexpand vale true hará que las raíces
19 16 CAPÍTULO. DESPEGANDO CON MAXIMA n -ésimas de los factores de un producto que sean potencias de n sean extraídas del radical; por ejemplo, sqrt(16*x^) se convertirá en 4*x sólo si radexpand vale true. There are also some inexact matches for sqrt. Try?? sqrt to see them. (%o5) true nos instruye sobre la utilización de la función sqrt. Alternativamente, se puede utilizar para el mismo fin el operador interrogación (?), (%i6)? sqrt -- Función: sqrt (<x>) Raíz cuadrada de <x>. Se representa internamente por <x>^(1/). Véase también rootscontract. Si la variable radexpand vale true hará que las raíces n -ésimas de los factores de un producto que sean potencias de n sean extraídas del radical; por ejemplo, sqrt(16*x^) se convertirá en 4*x sólo si radexpand vale true. There are also some inexact matches for sqrt. Try?? sqrt to see them. (%o6) true Cuando no conocemos el nombre completo de la función sobre la que necesitamos información podremos utilizar el operador de doble interrogación (??) que nos hace la búsqueda de funciones que contienen en su nombre cierta subcadena, (%i7)?? sqrt 0: isqrt (Operadores generales) 1: sqrt (Operadores generales) : sqrtdenest (Funciones y variables para simplification) 3: sqrtdispflag (Operadores generales) Enter space-separated numbers, all or none : 3 -- Variable opcional: sqrtdispflag Valor por defecto: true
20 .. UNA PRIMERA SESIÓN 17 Si sqrtdispflag vale false, hará que sqrt se muestre con el exponente 1/. (%o7) true Inicialmente muestra un sencillo menú con las funciones que contienen la subcadena sqrt y en el que debemos seleccionar qué información concreta deseamos; en esta ocasión se optó por 3 que es el número asociado a la variable opcional sqrtdispflag. A fin de no perder los resultados de una sesión con Maxima, quizás con el objeto de continuar el trabajo en próximas jornadas, podemos guardar en un archivo aquellos valores que nos puedan interesar; supongamos que necesitamos almacenar el contenido de la variable y, tal como lo definimos en la entrada (%i5), así como el resultado de la salida (%o9), junto con la instrucción que la generó, la entrada (%i9), (%i8) save("mi.sesion",y,resultado=%o9,ecuacion=%i9)$ Véase que el primer argumento de save es el nombre del archivo, junto con su ruta si se considera necesario, que la variable y se escribe tal cual, pero que las referencias de las entradas y salidas deben ir acompañadas de un nombre que las haga posteriormente reconocibles. Por último, la forma correcta de abandonar la sesión de Maxima es mediante (%i9) quit(); Abramos ahora una nueva sesión con Maxima y carguemos en memoria los resultados de la anterior sesión, (%i1) load("mi.sesion")$ (%i) y; ( %o) (%i3) ecuacion; ( %o3) solve ( x 3 x + 1 = 0, x ) (%i4) resultado; ( %o4) [ ] x =, x = Si ahora quisiéramos recalcular las soluciones de la ecuación, modificándola previamente de manera que el coeficiente de primer grado se sustituyese por otro número, simplemente haríamos
21 18 CAPÍTULO. DESPEGANDO CON MAXIMA (%i5) subst(5,3,ecuacion); ( %o5) solve ( x 5 x + 1 = 0, x ) (%i6) ev(%, nouns); ( %o6) [ 1 5 x =, x = ] La función save guarda las expresiones matemáticas en formato Lisp, pero a nivel de usuario quizás sea más claro almacenarlas en el propio lenguaje de Maxima; si tal es el caso, mejor será el uso de stringout, cuyo resultado será un fichero que también se podrá leer mediante load. Continuemos experimentando con esta segunda sesión de Maxima que acabamos de iniciar. Maxima está escrito en Lisp, por lo que ya se puede intuir su potencia a la hora de trabajar con listas; en el siguiente diálogo, el texto escrito entre las marcas /* y */ son comentarios que no afectan a los cálculos, (%i7) /* Se le asigna a la variable x una lista */ x: [cos(%pi), 4/16, [a, b], (-1)^3, integrate(u^,u)]; ( %o7) [ 1, ] 14 u3, [a, b], 1, 3 (%i8) /* Se calcula el número de elementos del último resultado */ length(%); ( %o8) 5 (%i9) /* Se transforma una lista en conjunto, */ /* eliminando redundancias. Los */ /* corchetes se transforman en llaves */ setify(x); ( %o9) { 1, } 14 u3, [a, b], 3 Vemos a continuación cómo se procede a la hora de hacer sustituciones en una expresión, (%i10) /* Sustituciones a hacer en x */ x, u=, a=c;
22 .. UNA PRIMERA SESIÓN 19 ( %o10) [ 1, 1 4, [c, b], 1, 8 ] 3 (%i11) /* Forma alternativa para hacer lo mismo */ subst([u=, a=c],x); ( %o11) [ 1, 1 4, [c, b], 1, 8 ] 3 (%i1) %o7[5], u=[1,,3]; ( %o1) [ 1 3, 8 ] 3, 9 En esta última entrada, %o7 hace referencia al séptimo resultado devuelto por Maxima, que es la lista guardada en la variable x, de modo que %o7[5] es el quinto elemento de esta lista, por lo que esta expresión es equivalente a a x[5]; después, igualando u a la lista [1,,3] este quinto elemento de la lista es sustituído sucesivamente por las tres cantidades. Tal como ya se ha comentado, Maxima utiliza los dos puntos (:) para hacer asignaciones a variables, mientras que el símbolo de igualdad se reserva para la construcción de ecuaciones. Un aspecto a tener en cuenta es que el comportamiento de Maxima está controlado por los valores que se le asignen a ciertas variables globales del sistema. Una de ellas es la variable numer, que por defecto toma el valor lógico false, lo que indica que Maxima evitará dar resultados en formato decimal, prefiriendo expresiones simbólicas y fraccionarias, (%i13) numer; ( %o13) false (%i14) sqrt(8)/1; ( %o14) (%i15) numer:true$ (%i16) sqrt(8)/1; 1 3 ( %o16) (%i17) numer:false$ /*se reinstaura valor por defecto */ Sin embargo, en general será preferible el uso de la función float,
23 0 CAPÍTULO. DESPEGANDO CON MAXIMA (%i18) float(sqrt(8)/1); ( %o18) Las matrices también tienen su lugar en Maxima; la manera más inmediata de construirlas es con la instrucción matrix, (%i19) A: matrix([sin(x),,%pi],[0,y^,5/8]); ( %o19) ( ) sin x π 0 y 5 8 (%i0) B: matrix([1,],[0,],[-5,integrate(x^,x,0,1)]); ( %o0) (%i1) A.B; /* producto matricial */ ( %o1) ( sin x 5 π sin x + π ) 5 8 y En cuanto a las capacidades gráficas de Maxima, se remite al lector a la sección correspondiente. A fin de preparar publicaciones científicas, si se quiere redactar documentos en formato L A TEX, la función tex de Maxima será una útil compañera; en el siguiente ejemplo se calcula una integral y a continuación se pide la expresión resultante en formato TEX: (%i) integrate(sqrt(a+x)/x^5,x,1,) = integrate(sqrt(+x)/x^5,x,1,)$ (%i3) tex(%); $$\int_{1}^{}{{{\sqrt{x+4}}\over{x^5}}\;dx}={{5\,\sqrt{}\,\log \left(5-\,\sqrt{}\,\sqrt{3}\right)+548\,\sqrt{3}}\over{048}}-{{15 \,\sqrt{}\,\log \left(3-\,\sqrt{}\right)+44}\over{6144}}$$ Al pegar y copiar el código encerrado entre los dobles símbolos de dólar a un documento L A TEX, el resultado es el siguiente: 1 x + a x 5 dx = 5 ( ) log log ( 3 )
24 .. UNA PRIMERA SESIÓN 1 Un comentario sobre la entrada %i. Se observa en ella que se ha escrito dos veces el mismo código a ambos lados del signo de igualdad. La única diferencia entre ambos es el apóstrofo que antecede al primer integrate; éste es el operador de comilla simple, ya anteriormente aparecido, que Maxima utiliza para evitar la ejecución de una instrucción, de forma que en el resultado la primera integral no se calcula, pero sí la segunda, dando lugar a la igualdad que se obtiene como resultado final. Las expresiones que en Maxima se marcan con la comilla simple para no ser evaluadas reciben el nombre de expresiones nominales, en contraposición a las otras, conocidas como expresiones verbales. Ya se ha comentado más arriba cómo solicitar ayuda de Maxima. Algunos usuarios encuentran engorroso que la ayuda se interfiera con los cálculos; esto se puede arreglar con un mínimo de bricolage informático. Al tiempo que se instala Maxima, se almacena también el sistema de ayuda en formato html en una determinada carpeta; el usuario tan sólo tendrá que buscar esta carpeta, en la que se encontrará la página principal maxima.html, y hacer uso de la función system de Maxima que ejecuta comandos del sistema, la cual admite como argumento una cadena alfanumérica con el nombre del navegador (en el ejemplo, mozilla) seguido de la ruta hacia el documento maxima.html, (%i4) system("(mozilla /usr/local/share/maxima/5.5.0/doc/html/es.utf8/maxima.html)&")$ Este proceso se puede automatizar si la instrucción anterior se guarda en un fichero de nombre maxima-init.mac y se guarda en la carpeta hacia la que apunta la variable maxima_userdir; por ejemplo, en Linux, (%i5) maxima_userdir; ( %o5) /home/xefe/.maxima El fichero maxima-init.mac se lee cada vez que arranca Maxima y se ejecutan las instrucciones contenidas en él; se puede utilizar para cargar automáticamente módulos adicionales que vayamos a utilizar con frecuencia, funciones especiales definidas por nosotros o, como en este caso, para que se abra el navegador con las páginas del manual.
25 CAPÍTULO. DESPEGANDO CON MAXIMA
26 Capítulo 3 Números y sus operaciones 3.1. Enteros y decimales Maxima puede trabajar con números enteros tan grandes como sea necesario, (%i1) 00!; ( %o1) Los operadores aritméticos más comunes son: + suma - resta * producto / división ^ o ** potencia Maxima devuelve los resultados de forma exacta, sin aproximaciones decimales; así tenemos que para realizar el cálculo de la expresión haremos ( ) (%i) ((3^7/(4+3/5))^(-1)+7/9)^3; 3 3
27 4 CAPÍTULO 3. NÚMEROS Y SUS OPERACIONES ( %o) obteniendo el resultado en forma de fracción simplificada. De todos modos podemos solicitar el resultado en forma decimal; por ejemplo, si queremos la expresión decimal del último resultado, (%i3) float(%); ( %o3) Maxima puede trabajar con precisión arbitraria. Para calcular el valor del cociente e π con 100 cifras decimales, deberemos especificar primero la precisión requerida asignándole a la variable fpprec el valor 100 y a continuación realizar el cálculo, solicitando la expresión decimal con una llamada a la función bfloat: (%i4) fpprec:100$ bfloat(%e / %pi); ( %o5) B 10 1 Nótese que cuando un resultado es devuelto en el formato bfloat se escribe la letra B en lugar de la clásica E para indicar el exponente. En la instrucción %i4 se establece la precisión deseada y luego se transforma el valor simbólico %e / %pi mediante la función bfloat. Recuérdese que el símbolo $ sirve como delimitador entre instrucciones. Veamos cómo factorizar un números tal como 00! en sus factores primos, (%i6) factor(%o1); ( %o6) En la siguiente factorización hacemos uso de la variable global showtime para saber el tiempo que le lleva ejecutar el cálculo, (%i7) showtime:true$ Evaluation took 0.00 seconds (0.00 elapsed) using 80 bytes. (%i8) factor(^300-1); Evaluation took 0.86 seconds (0.85 elapsed) using MB. ( %o8)
28 3.1. ENTEROS Y DECIMALES 5 (%i9) /* desactivamos la información de los tiempos */ showtime: false$ El texto que se escribe entre /* y */ son comentarios que Maxima ignora. Además de la variable opcional showtime, otro mecanismo de control de los tiempos de ejecución es la función time, (%i10) time(%o8); ( %o10) [ ] En relación con los números primos, para saber si un número lo es o no, (%i11) primep(3^56-1); ( %o11) false Y para solicitarle a Maxima que compruebe si un número es par o impar necesitamos echar mano de las funciones evenp or oddp, respectivamente, (%i1) evenp(4); ( %o1) true (%i13) oddp(31); ( %o13) true Le solicitamos a Maxima que nos haga una lista con todos los divisores de un número, (%i14) divisors(100); ( %o14) {1,, 4, 5, 10, 0, 5, 50, 100} En lo que concierne a la división, puede ser necesario conocer el cociente entero y el resto correspondiente; para ello se dispone de las funciones quotient y remainder, (%i15) quotient(5689,14); ( %o15) 6 (%i16) remainder(5689,14); ( %o16) 15
29 6 CAPÍTULO 3. NÚMEROS Y SUS OPERACIONES 3.. Números complejos Como ya se comentó más arriba, la unidad imaginaria 1 se representa en Maxima mediante el símbolo %i, (%i1) %i^; ( %o1) 1 Se soportan las operaciones básicas con números complejos, (%i) z1:3+5*%i$ z:1/-4*%i$ (%i4) z1 + z; /*suma*/ /*z1 y z*/ ( %o4) i + 7 (%i5) z1 - z; /*resta*/ ( %o5) 9 i + 5 (%i6) z1 * z; ( %o6) /*multiplicación*/ ( ) 1 4 i (5 i + 3) (%i7) z1 / z; /* división */ ( %o7) 5 i i Es posible que estos dos últimos resultados nos parezcan frustrantes y los deseemos en otro formato; a tal efecto podemos pedir a Maxima que nos devuelva (%o6) y (%o7) en forma cartesiana, (%i8) rectform(%o6); ( %o8) (%i9) rectform(%o7); i ( %o9) 58 i Las funciones realpart e imagpart extraen del número complejo sus partes real e imaginaria, respectivamente,
30 3.. NÚMEROS COMPLEJOS 7 (%i10) realpart(%o7); ( %o10) (%i11) imagpart(%o7); ( %o11) Antes hemos optado por los resultados en formato cartesiano; pero también es posible solicitarlos en su forma polar, (%i1) polarform(%o7); ( %o1) 9 i 34 e (π arctan( 37)) 65 (%i13) rectform(%); /* forma rectangular */ ( %o13) i (%i14) float(%); /* formato decimal */ ( %o14) i El módulo de un número complejo se calcula con la función abs y se le puede pasar el número tanto en su forma cartesiana como polar, (%i15) abs(%o9); ( %o15) (%i16) abs(%o1); ( %o16) El argumento de un número complejo se calcula con la función carg, que devuelve un ángulo en el intervalo ( π, π], (%i17) carg(%o9);
31 8 ( %o17) π arctan CAPÍTULO 3. NÚMEROS Y SUS OPERACIONES ( ) 9 37 Por último, el conjugado de un número complejo se calcula con la función conjugate. Como se ve en este ejemplo, el resultado dependerá del formato del argumento, cartesiano o polar, (%i18) conjugate(%o9); /*forma cartesiana*/ ( %o18) 58 i (%i19) conjugate(%o1); /*forma polar*/ ( %o19) 9 i 34 e arctan( 37) Combinatoria Ya hemos visto el cálculo del factorial de un número natural, bien con el operador!, bien con la función factorial; así podemos calcular el número de permutacions de n elementos, P n = n!, (%i1) /* permutaciones de 3 elementos */ factorial(3); ( %o1) Esto nos da el número de ordenaciones diferentes que se pueden hacer con n objetos; si queremos generarlas podemos hacer uso de la función permutations, (%i) permutations([1,,3]); ( %o) {[1,, 3], [1, 3, ], [, 1, 3], [, 3, 1], [3, 1, ], [3,, 1]} Cuando se trabaja con factoriales, la función minfactorial puede ayudarnos a simplificar algunas expresiones, (%i3) n!/(n+)!; ( %o3) n! (n + )!
32 3.3. COMBINATORIA 9 (%i4) minfactorial (%); ( %o4) 1 (n + 1) (n + ) El cálculo del número de variaciones con repetición de m elementos tomados de n en n, V R n m = m n, no necesita de ninguna función especial, (%i5) /* VR de 17 tomados de 4 en 4 */ 17^4; ( %o5) 8351 Sin embargo, para calcular el número de variaciones sin repetición de m elementos tomados de n en n, V n m = m!/(m n)!, debemos hacer uso de la función permutation definida en el módulo functs, que tendremos que cargar en memoria previamente, (%i6) load("functs")$ (%i7) /* V de 17 tomados de 4 en 4 */ permutation(17,4); ( %o7) 5710 Nótese que podemos calcular P 3 con esta misma función, sin llamar al factorial, (%i8) /* permutaciones de 3 elementos, otra vez */ permutation(3,3); ( %o8) El número de combinaciones de m elementos tomados de n en n, C n m = V n m/p n, se puede calcular con la función binomial, (%i9) /* combinaciones de 30 elementos, tomados de 5 en 5 */ binomial(30,5); ( %o9) El número de permutaciones de n elementos, con repeticiones a 1,..., a k, P a 1,...,a r n = (a a r )!, a 1!... a r! siendo a a r = n, se calcula con la función multinomial_coeff,
33 30 CAPÍTULO 3. NÚMEROS Y SUS OPERACIONES (%i10) /* permutaciones con repetición de 5 elementos, con repeticiones, y 1 */ multinomial_coeff(,,1); ( %o10) 30 La función stirling1 calcula el número de Stirling de primera especie, s(n, k), que se define de la forma ( 1) n k s(n, k) = Card{σ σ es permutación de {1,..., n} σ tiene k ciclos} (%i11) stirling1(5,); ( %o11) 50 Por otro lado, la función stirling calcula el número de Stirling de segunda especie, S(n, k), que se define de la forma S(n, k) = Card{π π es partición de {1,..., n} Card(π) = k} (%i1) stirling(3,); ( %o1) 3 Esto es, hay 3 particiones de {1,, 3} formadas por dos subconjuntos, las cuales podemos obtener invocando la función set_partitions, (%i13) set_partitions ({1,,3}, ); ( %o13) {{{1}, {, 3}}, {{1, }, {3}}, {{1, 3}, {}}} El número de Bell, que calculamos con la función belln es el número total de particiones que admite un conjunto finito, B(n) = n S(n, k). k=1 (%i14) /* número de particiones en un conjunto de cardinalidad 3 */ belln(3); ( %o14) 5 En relación con la combinatoria, véanse también las secciones 4.4 y 4.5, sobre conjuntos y grafos, respectivamente.
34 Capítulo 4 Estructuras de datos 4.1. Listas Las listas son objetos muy potentes a la hora de representar estructuras de datos; de hecho, toda expresión de Maxima se representa internamente como una lista, lo que no es de extrañar habida cuenta de que Maxima está programado en Lisp (List Processing). Veamos cómo podemos ver la representación interna, esto es en Lisp, de una sencilla expresión tal como 1 + 3a, (%i1)?print(1+3*a)$ ((MPLUS SIMP) 1 ((MTIMES SIMP) 3 $A)) Pero al nivel del usuario que no está interesado en las interioridades de Maxima, también se puede trabajar con listas como las definidas a continuación, siempre encerradas entre corchetes, (%i) r:[1,[a,3],sqrt(3)/,"don Quijote"]; ( %o) [ ] 3 1, [a, 3],, Don Quijote Vemos que los elementos de una lista pueden a su vez ser también listas, expresiones matemáticas o cadenas de caracteres incluidas entre comillas dobles, lo que puede ser aprovechado para la construcción y manipulación de estructuras más o menos complejas. Extraigamos a continuación alguna información de las listas anteriores, (%i3) listp(r); /* es r una lista? */ ( %o3) true (%i4) first(r); /* primer elemento */ 31
35 3 CAPÍTULO 4. ESTRUCTURAS DE DATOS ( %o4) 1 (%i5) second(r); /* segundo elemento */ ( %o5) [a, 3] (%i6) third(r); /*...hasta tenth */ 3 ( %o6) (%i7) last(r); /* el último de la lista */ ( %o7) Don Quijote (%i8) rest(r); /* todos menos el primero */ ( %o8) [ ] 3 [a, 3],, Don Quijote (%i9) part(r,3); /* pido el que quiero */ 3 ( %o9) (%i10) length(r); /* cuántos hay? */ ( %o10) 4 (%i11) reverse(r); /* le damos la vuelta */ ( %o11) [ ] 3 Don Quijote,, [a, 3], 1 (%i1) member(a,r); /* es a un elemento?*/ ( %o1) false (%i13) member([a,3],r); /* lo es [a,3]? */ ( %o13) true
36 4.1. LISTAS 33 (%i14) /* asigno valor a q y ordeno */ sort(q: [7,a,1,d,5,3,b]); ( %o14) [1, 3, 5, 7, a, b, d] (%i15) delete([a,3],r); /* borro elemento */ ( %o15) [ ] 3 1,, Don Quijote Nótese que en todo este tiempo la lista r no se ha visto alterada, (%i16) r; ( %o16) [ ] 3 1, [a, 3],, Don Quijote Algunas funciones de Maxima permiten añadir nuevos elementos a una lista, tanto al principio como al final, (%i17) cons(1+,q); ( %o17) [3, 7, a, 1, d, 5, 3, b] (%i18) endcons(x,q); ( %o18) [7, a, 1, d, 5, 3, b, x] (%i19) q; ( %o19) [7, a, 1, d, 5, 3, b] En este ejemplo hemos observado también que la lista q no ha cambiado. Si lo que queremos es actualizar su contenido, (%i0) q: endcons(x,cons(1+,q))$ (%i1) q; ( %o1) [3, 7, a, 1, d, 5, 3, b, x] Es posible unir dos listas, (%i) append(r,q);
37 34 ( %o) CAPÍTULO 4. ESTRUCTURAS DE DATOS [ ] 3 1, [a, 3],, Don Quijote, 3, 7, a, 1, d, 5, 3, b, x Cuando los elementos de una lista van a obedecer un cierto criterio de construcción, podemos utilizar la función makelist para construirla, (%i3) s:makelist(+k*,k,0,10); ( %o3) [, 4, 6, 8, 10, 1, 14, 16, 18, 0, ] donde le hemos indicado a Maxima que nos construya una lista con elementos de la forma +*k, de modo que k tome valores enteros de 0 a 10. La función apply permite suministrar a otra función todos los elementos de una lista como argumentos, así podemos sumar o multiplicar todos los elementos de la lista s recién creada, (%i4) apply("+",s); ( %o4) 13 (%i5) apply("*",s); ( %o5) aunque estas dos operaciones hubiese sido quizás más apropiado haberlas realizado con las funciones sum y product. A veces interesará aplicar una misma función a varios elementos de una lista de forma independiente, para lo que haremos uso de map. A continuación un ejemplo de cálculo de factoriales y otro trigonométrico; en este último, no es necesario aplicar map, ya que algunas funciones se distribuyen automáticamente sobre los elementos de una lista, (%i6) map("!",s); ( %o6) [, 4, 70, 4030, , , , , , , ] (%i7) sin(s); ( %o7) [sin, sin 4, sin 6, sin 8, sin 10, sin 1, sin 14, sin 16, sin 18, sin 0, sin ] A veces se quiere utilizar en map una función que no está definida y que no va a ser utilizada más que en esta llamada a map. En tales casos se puede hacer uso de las funciones lambda; supóngase que a cada elemento x de una lista se le quiere aplicar la función f(x) = sin(x) + 1, sin necesidad de haberla definido previamente,
38 4.1. LISTAS 35 (%i8) map(lambda([x], sin(x) + 1/),[7,3,4,9]); ( %o8) [ sin 7 + 1, sin 3 + 1, sin 4 + 1, sin ] Como se ve, la función lambda admite dos argumentos; el primero es una lista (no se puede evitar la redundancia) de argumentos, siendo el segundo la expresión que indica qué se hace con los elementos de la lista anterior. Por último, las listas también pueden ser utilizadas en operaciones aritméticas. Algunas de las operaciones que se aplican a continuación no están definidas matemáticamente, por lo que Maxima se toma la libertad de definirlas por su cuenta, lo que redunda en un ahorro de código por parte del usuario y en un estilo de programación más compacto, (%i9) [1,,3]+[a,b,c]; ( %o9) [a + 1, b +, c + 3] (%i30) [1,,3]*[a,b,c]; ( %o30) [a, b, 3 c] (%i31) [1,,3]/[a,b,c]; ( %o31) [ 1 a, b, 3 ] c (%i3) [1,,3]-[a,b,c]; ( %o3) [1 a, b, 3 c] (%i33) [1,,3].[a,b,c]; /* producto escalar */ ( %o33) 3 c + b + a (%i34) [a,b,c]^3; ( %o34) [ a 3, b 3, c 3] (%i35) 3^[a,b,c];
39 36 ( %o35) CAPÍTULO 4. ESTRUCTURAS DE DATOS [3 a, 3 b, 3 c] En Maxima, las listas se pueden interpretar como vectores; para más información, acúdase a la Sección 6.6. Para que estas operaciones puedan realizarse sin problemas, la variable global listarith debe tomar el valor true, en caso contrario el resultado será bien distinto, (%i36) listarith:false$ (%i37) [1,,3]+[4,5,6]; ( %o37) [4, 5, 6] + [1,, 3] (%i38) listarith:true$ Como ya se vió al comienzo de esta sección, una lista puede ser elemento de otra lista, si queremos deshacer todas las listas interiores para que sus elementos pasen a formar parte de la exterior, (%i39) flatten([1,[a,b],,3,[c,[d,e]]]); ( %o39) [1, a, b,, 3, c, d, e] El cálculo del módulo de un vector se puede hacer mediante la definición previa de una función al efecto: (%i40) modulo(v):= if listp(v) then sqrt(apply("+",v^)) else error("mucho ojito: ", v, " no es un vector!!!!")$ (%i41) xx:[a,b,c,d,e]$ (%i4) yy:[3,4,-6,0,4/5]$ (%i43) modulo(xx-yy); ( %o43) ( e 4 ) + d 5 + (c + 6) + (b 4) + (a 3) 4.. Arrays Los arrays son estructuras que almacenan datos secuencial o matricialmente. Se diferencian de las listas en que éstas pueden alterar su estructura interna y los arrays no; pero la ventaja frente a las listas estriba en que el acceso a un dato cualquiera es directo, mientras que en las listas, al ser estructuras dinámicas, para acceder al n-ésimo dato, un puntero debe recorrer
40 4.. ARRAYS 37 todos los elementos de la lista desde el que ocupa el primer lugar. Así pues, se recomienda el uso de arrays cuando se vaya a utilizar una secuencia larga de expresiones cuya estructura se va a mantener inalterable; en general, se mejorarán los tiempos de procesamiento. Si queremos construir un array, lo debemos declarar previamente mediante la función array, indicando su nombre y dimensiones; a continuación le asignaríamos valores a sus elementos, (%i1) array(arreglillo,,); ( %o1) arreglillo (%i) arreglillo[0,0]: 34$ (%i3) arreglillo[1,]:x+y$ (%i4) arreglillo[,]:float(%pi)$ (%i5) listarray(arreglillo); ( %o5) [34, #####, #####, #####, #####, y + x, #####, #####, ] En el ejemplo, hemos definido un array de dimensiones 3 por 3, ya que la indexación de los arrays comienza por el cero, no por el uno como en el caso de las listas; a continuación hemos asignado valores a tres elementos del array: un entero, una expresión simbólica y un número decimal en coma flotante. Finalmente, la función listarray nos permite ver el contenido del array, utilizando el símbolo ##### para aquellos elementos a los que aún no se les ha asignado valor alguno. Si ya se sabe de antemano que el array va a contener únicamente números enteros o números decimales, podemos declararlo al definirlo para que Maxima optimice el procesamiento de los datos. A continuación se declara un array de números enteros haciendo uso del símbolo fixnum (en el caso de números decimales en coma flotante se utilizaría flonum) y a continuación la función fillarray rellenará el array con los elementos previamente almacenados en una lista, (%i6) lis: [6,0,-4,456,56,-99]; ( %o6) [6, 0, 4, 456, 56, 99] (%i7) array(arr,fixnum,5); ( %o7) arr (%i8) fillarray(arr,lis);
41 38 CAPÍTULO 4. ESTRUCTURAS DE DATOS ( %o8) arr (%i9) arr[4]; ( %o9) 56 (%i10) listarray(arr); ( %o10) [6, 0, 4, 456, 56, 99] Hasta ahora hemos hecho uso de los llamados array declarados. Si no hacemos la declaración pertinente mediante la función array, y en su lugar realizamos asignaciones directas a los elementos de la estructura, creamos lo que se conoce como array de claves (hash array) o, en la terminología de Maxima, array no declarado; estos arrays son mucho más generales que los declarados, puesto que pueden ir creciendo de forma dinámica y los índices no tienen que ser necesariamente números enteros, (%i11) tab[9]:76; ( %o11) 76 (%i1) tab[-58]:/9; ( %o1) (%i13) tab[juan]:sqrt(8); 9 ( %o13) 3 (%i14) tab["cadena"]: a-b; ( %o14) a b (%i15) tab[x+y]:5; ( %o15) 5 (%i16) listarray(tab); ( %o16) [ ] 9, 76, a b, 3, 5
42 4.3. CADENAS 39 (%i17) tab["cadena"]; ( %o17) a b Esta generalidad en la estructura de los arrays no declarados tiene la contrapartida de restarle eficiencia en el procesamiento de los datos almacenados, siendo preferible, siempre que sea posible, la utilización de los declarados. Ya por último, junto con los arrays declarados y no declarados mencionados, es posible trabajar a nivel de Maxima con arrays de Lisp. En este caso su construcción se materializa con la función make_array, (%i18) a: make_array(flonum,5); ( %o18) {Array : # ( )} (%i19) a[0]:3.4; ( %o19) 3.4 (%i0) a[4]:6.89; ( %o0) 6.89 (%i1) a; ( %o1) {Array : # ( )} Nótense aquí un par de diferencias respecto de la función array: el número 5 indica el número de elementos (indexados de 0 a 4) y no el índice máximo admisible; por otro lado, la función no devuelve el nombre del objeto array, sino el propio objeto Cadenas Una cadena se construye escribiendo un texto o cualquier otra expresión entre comillas dobles. Vamos a hacer algunos ejemplos básicos, (%i1) xx: "Los "; ( %o1) Los (%i) yy: " cerditos";
43 40 CAPÍTULO 4. ESTRUCTURAS DE DATOS ( %o) cerditos (%i3) zz: concat(xx,sqrt(9),yy); ( %o3) Los 3 cerditos la función concat admite varios parámetros, los cuales evalúa y luego concatena convirtiéndolos en una cadena. La función stringp nos indica si su argumento es o no una cadena, (%i4) stringp(zz); ( %o4) true Para saber el número de caracteres de una cadena haremos uso de slength (%i5) slength(zz); ( %o5) 14 Una vez transformada en cadena una expresión de Maxima, deja de ser evaluable, (%i6) concat(sqrt(5)); ( %o6) 5 (%i7) % + ; ( %o7) 5 + (%i8) stringp(%); ( %o8) false el resultado no se simplifica porque uno de los sumandos es una cadena y no se reconoce como cantidad numérica. La función charat devuelve el carácter que ocupa una posición determinada, (%i9) charat(zz,7); ( %o9) c También podemos solicitar que Maxima nos devuelva una lista con todos los caracteres o palabras que forman una cadena,
44 4.3. CADENAS 41 (%i10) charlist(zz); ( %o10) [L, o, s,, 3,, c, e, r, d, i, t, o, s] (%i11) tokens(zz); ( %o11) [Los, 3, cerditos] O podemos transformar letras minúsculas en mayúsculas y viceversa, todas las letras o sólo un rango de ellas, (%i1) /* todo a mayúsculas */ supcase(zz); ( %o1) LOS 3 CERDITOS (%i13) /* a minúsculas entre la 9a y 1a */ sdowncase(%,9,1); ( %o13) LOS 3 CErdiTOS (%i14) /* a minúsculas de la 13a hasta el final */ sdowncase(%,13); ( %o14) LOS 3 CErdiTos (%i15) /* todo a minúsculas */ sdowncase(%); ( %o7) los 3 cerditos (%i16) /* a mayúsculas entre la 1a y a */ supcase(%,1,); ( %o74) Los 3 cerditos Para comparar si dos cadenas son iguales disponemos de la función sequal, (%i17) tt: zz; ( %o17) Los 3 cerditos (%i18) sequal(tt,zz);
45 4 CAPÍTULO 4. ESTRUCTURAS DE DATOS ( %o18) true (%i19) tt: supcase(zz); ( %o19) LOS 3 CERDITOS (%i0) sequal(tt,zz); ( %o0) false Podemos buscar subcadenas dentro de otras cadenas mediante la función ssearch; los siguientes ejemplos pretenden aclarar su uso, (%i1) zz; ( %o1) Los 3 cerditos (%i) /* busca la subcadena "cer" */ ssearch("cer",zz); ( %o) 7 (%i3) /* busca la subcadena "Cer" */ ssearch("cer",zz); ( %o3) false (%i4) /* busca "Cer" ignorando tama~no */ ssearch("cer",zz, sequalignore); ( %o4) 7 (%i5) /* busca sólo entre caracteres 1 y 5 */ ssearch("cer",zz, sequalignore,1,5); ( %o5) false Igual que buscamos subcadenas dentro de otra cadena, también podemos hacer sustituciones; en este caso la función a utilizar es ssubst, (%i6) ssubst("mosqueteros","cerditos",zz);
46 4.3. CADENAS 43 ( %o6) Los 3 mosqueteros Esta función tiene como mínimo tres argumentos, el primero indica la nueva cadena a introducir, el segundo es la subcadena a ser sustituída y el tercero la cadena en la que hay que realizar la sustitución. Es posible utilizar esta función con más argumentos; como siempre, para más información,? ssubst. Con la función substring podemos obtener una subcadena a partir de otra, para lo cual debemos indicar los extremos de la subcadena, (%i7) substring(zz,7); ( %o7) cerditos (%i8) substring(zz,3,8); ( %o8) s 3 c Como se ve, el segundo argumento indica la posición a partir de la cual queremos la subcadena y el tercero la última posición; si no se incluye el tercer argumento, se entiende que la subcadena se extiende hasta el final de la de referencia. En algunos contextos el usuario tiene almacenada una expresión sintácticamente correcta de Maxima en formato de cadena y pretende evaluarla, en cuyo caso la función eval_string la analizará y evaluará, en tanto que parse_string la analiza pero no la evalúa, limitándose a transformarla a una expresión nominal de Maxima, (%i9) /* una cadena con una expresión de Maxima */ expr: "integrate(x^5,x)"; ( %o9) integrate(xˆ5,x) (%i30) /* se evalúa la cadena */ eval_string(expr); ( %o30) (%i31) /* se transforma en una expresión nominal */ parse_string(expr); ( %o31) integrate ( x 5, x ) (%i3) %; x 6 6
47 44 CAPÍTULO 4. ESTRUCTURAS DE DATOS ( %o3) x 6 6 Véase que el resultado de parse_string ya no es una cadena, sino una expresión correcta de Maxima sin evaluar, cuya evaluación se solicita a continuación con la doble comilla simple. Para más información sobre el procesamiento de cadenas de texto, se sugiere la lectura del capítulo stringproc del manual de referencia Conjuntos Se define a continuación un conjunto mediante la función set, (%i1) c1:set(a,[,k],b,sqrt(),a,set(a,b), 3,"Sancho",set(),b,sqrt(),a); ( %o1) { 3, }, {}, [, k], Sancho, a, {a, b}, b Como se ve, se admiten objetos de muy diversa naturaleza como elementos de un conjunto: números, expresiones, el conjunto vacío ({}), listas, otros conjuntos o cadenas de caracteres. Se habrá observado que al construir un conjunto se eliminan las repeticiones de elementos, cual es el caso del símbolo a en la lista anterior. Cuando se trabaja con listas, puede ser de utilidad considerar sus componentes como elementos de un conjunto, luego se necesita una función que nos transforme una lista en conjunto, (%i) [[,k],sqrt(),set(b,a),[k,],"panza"]; ( %o) [ [, k], ], {a, b}, [k, ], Panza (%i3) c:setify(%); ( %o3) {, [, k], Panza, {a, b}, [k, ] } el cambio en la naturaleza de estas dos colecciones de objetos se aprecia en la presencia de llaves en lugar de corchetes. De igual manera, podemos transformar un conjunto en lista, (%i4) listify(%o1); ( %o4) [ 3, ], {}, [, k], Sancho, a, {a, b}, b Comprobemos de paso que {} representa al conjunto vacío,
48 4.4. CONJUNTOS 45 (%i5) emptyp(%[3]); ( %o5) true Recuérdese que % sustituye a la última respuesta dada por Maxima, que en este caso había sido una lista, por lo que %[3] hace referencia a su tercera componente. Para comprobar si un cierto objeto forma parte de un conjunto hacemos uso de la instrucción elementp, (%i6) elementp(sqrt(),c1); ( %o6) true Es posible extraer un elemento de un conjunto y luego añadirle otro distinto (%i7) c1: disjoin(sqrt(),c1); /* sqrt() fuera */ ( %o7) {3, {}, [, k], Sancho, a, {a, b}, b} (%i8) c1: adjoin(sqrt(3),c1); /* sqrt(3) dentro */ ( %o8) { 3, } 3, {}, [, k], Sancho, a, {a, b}, b La sustitución que se acaba de realizar se pudo haber hecho con la función subst, (%i9) /* nuevamente a poner sqrt() */ subst(sqrt(),sqrt(3),c1); ( %o9) { 3, }, {}, [, k], Sancho, a, {a, b}, b La comprobación de si un conjunto es subconjunto de otro se hace con la fución subsetp, (%i10) subsetp(set([k,],"panza"),c); ( %o10) true A continuación algunos ejemplos de operaciones con conjuntos, (%i11) union(c1,c);
49 46 ( %o11) CAPÍTULO 4. ESTRUCTURAS DE DATOS { 3,, } 3, {}, [, k], Panza, Sancho, a, {a, b}, b, [k, ] (%i1) intersection(c1,c); ( %o1) {[, k], {a, b}} (%i13) setdifference(c1,c); ( %o13) { 3, } 3, {}, Sancho, a, b (%i14) cardinality(%); ( %o14) 6 Vemos aquí también cómo pedir el cardinal de un conjunto. Igual que se ha visto cómo aplicar una función a todos los elementos de una lista, podemos hacer lo mismo con los de un conjunto, (%i15) map(sin,set(1,,3,4,5)); ( %o15) {sin 1, sin, sin 3, sin 4, sin 5} Por último ya, el producto cartesiano de tres conjuntos, (%i16) cartesian_product(set(1,),set(a,b,c),set(x,y)); ( %o16) {[1, a, x], [1, a, y], [1, b, x], [1, b, y], [1, c, x], [1, c, y], [, a, x], [, a, y], [, b, x], [, b, y], [, c, x], [, c, y]} Ya por último, el conjunto de las partes de un conjunto, o conjunto potencia, lo podemos obtener con la función powerset (%i17) powerset({1,a,sqrt()}); ( %o17) { { {}, {1}, 1, } {, 1, } { } { } }, a, {1, a},,, a, {a}
50 4.5. GRAFOS Grafos El soporte para el desarrollo de algoritmos en el ámbito de la teoría de grafos se encuentra en el módulo graphs. Se trata de un paquete de funciones que aporta la infraestructura básica para el trabajo con grafos, orientados o no. No se admiten bucles ni aristas múltiples. (%i1) load(graphs)$ (%i) g: create_graph( [1,,3,4,5,6,7], [ [1,],[,3],[,4],[3,4], [4,5],[5,6],[4,6],[6,7] ]); ( %o) GRAPH Como vemos, la respuesta de la función create_graph no es muy informativa, pero nos podemos hacer una idea de la estructura del grafo haciendo (%i3) print_graph(g)$ Graph on 7 vertices with 8 edges. Adjacencies: 7 : 6 6 : : : : 4 : : Vemos que, por ejemplo, el vértice número 3 comparte aristas (no orientadas) con los vértices 4 y. En caso de que el grafo anterior fuese orientado (digrafo) su definición sería: (%i4) dg: create_graph( [1,,3,4,5,6,7], [ [1,],[,3],[,4],[3,4], [4,5],[5,6],[4,6],[6,7] ], directed = true)$ (%i5) print_graph(dg)$ Digraph on 7 vertices with 8 arcs. Adjacencies:
51 48 CAPÍTULO 4. ESTRUCTURAS DE DATOS a) b) c) d) 17 Figura 4.1: Grafos: a) no orientado; b) orientado; c) complemetario del no orientado; d) camino más corto. 7 : 6 : 7 5 : 6 4 : : 4 : : En este caso, print_graph nos muestra los nodos hijos a la derecha de los padres. En particular, el nodo 7 no es origen de ningún arco, pero es hijo del nodo 6. Otra manera alternativa de mostrar la estructura de un grafo es mediante la función draw_graph (%i6) draw_graph(g, show_id=true, terminal=eps)$ (%i7) draw_graph(dg, show_id=true, head_length=0.05, terminal=eps)$ Cuyos aspectos se pueden observar en los apartados a) y b) de la Figura 4.1. Continuemos con el grafo no orientado g definido más arriba y veamos cómo extraer de él cierta información: el número cromático, matriz de adyacencia, grafo complementario, si es conexo, o planar, su diámetro, etc.
52 4.5. GRAFOS 49 (%i8) chromatic_number(g); ( %o8) 3 (%i9) adjacency_matrix(g); ( %o9) (%i10) gc:complement_graph(g); ( %o10) GRAPH (%i11) print_graph(%)$ Graph on 7 vertices with 13 edges. Adjacencies: 1 : : : : : : : Pedimos ahora a Maxima que nos represente gráficamente el complementario de g, cuyo aspecto es el del apartado c) de la Figura 4.1. (%i1) draw_graph(gc, show_id=true, terminal=eps)$ Continuamos con algunas otras propiedades: (%i13) is_connected(g); ( %o13) true
53 50 CAPÍTULO 4. ESTRUCTURAS DE DATOS (%i14) is_planar(g); ( %o14) true (%i15) diameter(g); ( %o15) 4 Junto con la definición manual, también es posible generar determinados grafos mediante funciones, como el cíclico, el completo, el dodecaedro, un grafo aleatorio, etc. (%i16) print_graph(g1: cycle_graph(5))$ Graph on 5 vertices with 5 edges. Adjacencies: 4 : : 4 : : 0 0 : 4 1 (%i17) print_graph(g: complete_graph(4))$ Graph on 4 vertices with 6 edges. Adjacencies: 3 : 1 0 : : : 3 1 (%i18) print_graph(g3: dodecahedron_graph())$ Graph on 0 vertices with 30 edges. Adjacencies: 19 : : : : : :
54 4.5. GRAFOS : : : : : : : : : : : 17 4 : : : (%i19) /* 0.4 de probabilidad para cada arco */ print_graph(g4: random_graph(10,0.4))$ Graph on 10 vertices with 15 edges. Adjacencies: 9 : : : : : : 3 : : : : Por último ya, un ejemplo de cálculo de la ruta más corta. Se genera un grafo orientado aleatorio, a continuación se calcula la ruta más corta entre dos vértices dados y finalmente se hace una representación gráfica del cálculo realizado, la cual se puede ver en el apartado d) de la Figura 4.1. (%i0) r : random_digraph(0, 0.5)$ (%i1) /* Una lista de los vértices del digrafo */ vertices(r); ( %o1) [19, 18, 17, 16, 15, 14, 13, 1, 11, 10, 9, 8, 7, 6, 5, 4, 3,, 1, 0]
55 5 CAPÍTULO 4. ESTRUCTURAS DE DATOS (%i) sh: shortest_path(0,14,r); ( %o) [0, 4, 1, 17, 14] (%i3) draw_graph( r, show_edges show_id head_length = 0.05, terminal = eps)$ = vertices_to_path(sh), = true, El paquete graphs también permite el uso de etiquetas en los nodos y ponderaciones en las aristas y arcos. El lector interesado podrá consultar la documentación correspondiente en el manual de referencia.
56 Capítulo 5 Álgebra 5.1. Transformaciones simbólicas Expresiones algebraicas Sin duda, una de las capacidades más destacables de Maxima es su habilidad para manipular expresiones algebraicas. Desarrollemos un ejemplo que empieza por asignar a la variable q una expresión literal: (%i1) q: (x+3)^5-(x-a)^3+(x+b)^(-1)+(x-1/4)^(-5); ( %o1) 1 x + b (x a)3 + (x + 3) ( x 1 4 ) 5 Se observa que en principio Maxima no realiza ningún cálculo. La función expand se encarga de desarrollar las potencias y los productos algebraicos, (%i) expand(q); ( %o) 1 x 5 5 x x3 8 5 x x x+b + x x x a x + 70 x 3 a x x + a No obstante es posible que no nos interese desplegar toda la expresión, entre otras cosas para evitar una respuesta farragosa y difícil de interpretar; en tal caso podemos utilizar expand añadiéndole dos argumentos y operar de la manera siguiente (%i3) expand(q,3,); ( %o3) 1 x + b + (x + 3)5 x a x 3 a x + 1 ( x 1 4 ) 5 + a 3 53
57 54 CAPÍTULO 5. ÁLGEBRA Con los argumentos adicionales indicamos que queremos la expansión de todas aquellas potencias con exponente positivo menor o igual a 3 y de las que teniendo el exponente negativo no excedan en valor absoluto de. Dada una expresión con valores literales, podemos desear sustituir alguna letra por otra expresión; por ejemplo, si queremos hacer los cambios a =, b = c en el último resultado obtenido, (%i4) %,a=,b=*c; ( %o4) 1 x + c + (x + 3)5 x x 1 x + 1 ( x 1 4 ) En el ejemplo %i4 se presenta una sentencia en la que hay elementos separados por comas (,). Esta es una forma simplificada de utilizar la función ev, que evalúa la primera expresión asignando los valores que se le van indicando a continuación; por ejemplo, (%i4) se podría haber escrito como ev(%,a=,b=*c). El uso de la variante con ev está más indicado para ser utilizado dentro de expresiones más amplias. Obsérvese el resultado siguiente (%i5) 3*x^ + ev(x^4,x=5); ( %o5) 3 x + 65 donde la sustitución x = 5 se ha realizado exclusivamente dentro del contexto delimitado por la función ev. De forma más general, la función subst sustituye subexpresiones enteras. En el siguiente ejemplo, introducimos una expresión algebraica y a continuación sustituimos todos los binomios x+y por la letra k, (%i6) 1/(x+y)-(y+x)/z+(x+y)^; ( %o6) y + x z (%i7) subst(k,x+y,%); + (y + x) + 1 y + x ( %o7) k z + k + 1 k No obstante, el siguiente resultado nos sugiere que debemos ser precavidos con el uso de esta función, ya que Maxima no siempre interpretará como subexpresión aquella que para nosotros sí lo es: (%i8) subst(sqrt(k),x+y,(x+y)^+(x+y));
58 5.1. TRANSFORMACIONES SIMBÓLICAS 55 ( %o8) y + x + k La función subst realiza una simple sustitución sintáctica que, si bien es suficiente en la mayor parte de los casos, a veces se queda corta como en la situación anterior; la función ratsubst es capaz de utilizar información semántica y obtene mejores resultados, (%i9) ratsubst(sqrt(k),x+y,(x+y)^+(x+y)); ( %o9) k + k Como una aplicación práctica de subst, veamos cómo podemos utilizarla para obtener el conjugado de un número complejo, (%i10) subst(-%i,%i,a+b*%i); ( %o10) a i b La operación inversa de la expansión es la factorización. Expandamos y factoricemos sucesivamente un polinomio para comprobar los resultados, (%i11) expand((a-)*(b+1)^*(a+b)^5); ( %o11) a b 7 b a b 6 8 a b 6 4 b a 3 b 5 10 a b 5 19 a b 5 b a 4 b 4 35 a b 4 10 a b a 5 b a 4 b 3 30 a 3 b 3 0 a b 3 + a 6 b + 8 a 5 b 10 a 4 b 0 a 3 b + a 6 b + a 5 b 10 a 4 b + a 6 a 5 (%i1) factor(%); ( %o1) (a ) (b + 1) (b + a) 5 El máximo común divisor de un conjunto de polinomios se calcula con la función gcd y el mínimo común múltiplo con lcm (%i13) p1: x^7-4*x^6-7*x^5+8*x^4+11*x^3-4*x^-5*x; ( %o13) x 7 4 x 6 7 x x x 3 4 x 5 x (%i14) p: x^4-*x^3-4*x^+*x+3; ( %o14) x 4 x 3 4 x + x + 3 (%i15) gcd(p1,p);
59 56 CAPÍTULO 5. ÁLGEBRA ( %o15) x 3 + x x 1 (%i16) load(functs)$ (%i17) lcm(p1,p); ( %o17) (x 5) (x 3) (x 1) x (x + 1) 3 En (%i13) y (%i14) definimos los polinomios p1 y p, a continuación calculamos su máximo común divisor (mcd) en (%i15) y antes de pedir el mínimo común múltiplo en (%i17) cargamos el paquete functs en el que se encuentra definida la función lcm. Es posible que deseemos disponer del mcd factorizado, por lo que hacemos (%i18) factor(%o15); ( %o18) (x 1) (x + 1) Hay varias funciones en Maxima que se pueden utilizar para extraer coeficientes de expresiones algebraicas; una de ellas es ratcoef. Dado el polinomio (x y) 4, podemos querer extraer el coeficiente de x y, (%i19) ratcoef((x-y)^4, x^*y^); ( %o19) 6 Cuando tenemos expandido un polinomio de varias variables, podemos querer transformarlo en otro univariante con coeficientes literales; aquí puede ser de ayuda la función rat, (%i0) q: m^*x^+x^+*b*m*x-6*m*x-4*x+b^-6*b+13; ( %o0) m x + x + b m x 6 m x 4 x + b 6 b + 13 (%i1) rat(q, x); ( %o1) ( m + 1 ) x + (( b 6) m 4) x + b 6 b + 13 Llamando nuevamente a ratcoef podemos extraer el coeficiente del término de primer grado respecto de x, (%i) ratcoef(%,x); ( %o) ( b 6) m 4 Otro ejemplo ilstrativo de rat en acción,
60 5.1. TRANSFORMACIONES SIMBÓLICAS 57 (%i3) rat((x^b+a*x^b+x)^, x^b); ( %o3) ( a + a + 1 ) ( x b) + ( a + ) x x b + x Si en algún momento queremos inhibir las simplificaciones, podemos hacer uso de la variable global simp, cuyo valor por defecto es true, (%i4) simp: false$ (%i5) + + a + a; ( %o5) + + a + a (%i6) simp: true$ (%i7) %o5; ( %o7) a + 4 Veamos ahora algunas transformaciones que se pueden hacer con fracciones algebraicas. (%i8) expr:(x^-x)/(x^+x-6)-5/(x^-4); ( %o8) x x x + x 6 5 x 4 Podemos factorizar esta expresión, (%i9) factor(expr); ( %o9) (x 3) ( x + 4 x + 5 ) (x ) (x + ) (x + 3) La expansión devuelve el siguiente resultado, (%i30) expand(expr); ( %o30) x x + x 6 x x + x 6 5 x 4 Ahora descomponemos la fracción algebraica expr en fracciones simples, (%i31) partfrac(expr,x);
61 58 CAPÍTULO 5. ÁLGEBRA 1 ( %o31) 5 (x + 3) (x + ) 17 0 (x ) + 1 Por último, transformamos la misma expresión a su forma canónica CRE (Cannonical Rational Expression), que es un formato que Maxima utiliza para reducir expresiones algebraicas equivalentes a una forma única, (%i3) radcan(expr); ( %o3) x 3 + x 7 x 15 x x 4 x 1 La función ratsimp también simplifica cualquier expresión racional, así como las subexpresiones racionales que son argumentos de funciones cualesquiera. El resultado se devuelve como el cociente de dos polinomios. En ocasiones no es suficiente con una sola ejecución de ratsimp, por lo que será necesario aplicarla más veces, esto es lo que hace precisamente la función fullratsimp; concretemos esto con un ejemplo: (%i33) (x^(a/)-1)^*(x^(a/)+1)^ / (x^a-1); ( %o33) ( ) x a ( ) 1 x a + 1 x a 1 (%i34) ratsimp(%); /* simplificamos una vez */ ( %o34) x a x a + 1 x a 1 (%i35) ratsimp(%); /* simplificamos otra vez */ ( %o35) x a 1 (%i36) fullratsimp(%o3); /* simplificamos todo de una vez! */ ( %o36) x a 1 Dada una fracción algebraica, podemos obtener separadamente el numerador y el denominador, (%i37) fr: (x^3-4*x^+4*x-)/(x^+x+1); ( %o37) x 3 4 x + 4 x x + x + 1
62 5.1. TRANSFORMACIONES SIMBÓLICAS 59 (%i38) num(fr); ( %o38) x 3 4 x + 4 x (%i39) denom(fr); ( %o39) x + x + 1 Maxima nunca puede adivinar a priori las intenciones del usuario para con las expresiones con las que está trabajando, por ello sigue la política de no tomar ciertas decisiones, por lo que a veces puede parecer que no hace lo suficiente o que no sabe hacerlo. Esta forma vaga de proceder la suple con una serie de funciones y variables globales que permitirán al usuario darle al motor simbólico ciertas directrices sobre qué debe hacer con las expresiones, cómo reducirlas o transformarlas. Sin ánimo de ser exhaustivos, siguen algunos ejemplos en los que se controlan transformaciones de tipo logarítmico y trigonométrico. Expresiones logarítmicas Las expresiones que contienen logaritmos suelen ser ambiguas a la hora de reducirlas. Transformamos a continuación la misma expresión según diferentes criterios, (%i1) log(x^r)-log(x*y) + a*log(x/y); ( %o1) log (x y) + a log ( ) x + r log x y En principio, ha hecho la transformación log x r r log x. Esto se controla con la variable global logexpand, cuyo valor por defecto es true, lo cual sólo permite la reducción recién indicada. Si también queremos que nos expanda los logaritmos de productos y divisiones, le cambiaremos el valor a esta variable global, tal como se indica a continuación, (%i) logexpand: super$ (%i3) log(x^r)-log(x*y) + a*log(x/y); ( %o3) log y + a (log x log y) + r log x log x (%i4) /* pedimos que nos simplifique la expresión anterior */ ratsimp(%);
63 60 CAPÍTULO 5. ÁLGEBRA ( %o4) ( a 1) log y + (r + a 1) log x En caso de no querer que nos haga transformación alguna sobre la expresión de entrada, simplemente haremos (%i5) logexpand: false$ (%i6) log(x^r)-log(x*y) + a*log(x/y); ( %o6) log (x y) + a log ( ) x + log x r y Devolvemos ahora la variable logexpand a su estado por defecto, (%i7) logexpand: true$ La función logcontract es la que nos va a permitir compactar expresiones logarítmicas, (%i8) logcontract(*(a*log(x) + *a*log(y))); ( %o8) a log ( x y 4) Por defecto, Maxima contrae los coeficientes enteros; para incluir en este caso la variable a dentro del logaritmo le asignaremos la propiedad de ser un número entero, (%i9) declare(a, integer)$ (%i10) %i8; /* fuerza reevaluación de expresión nominal */ ( %o10) log ( x a y 4 a) Esta es una forma de hacerlo, pero Maxima permite un mayor refinamiento sobre este particular haciendo uso de la variable global logconcoeffp. Expresiones trigonométricas Toda esta casuística de diferentes formas de representar una misma expresión aumenta considerablemente con aquéllas en las que intervienen las funciones trigonométricas e hiperbólicas. Empecemos con la función trigexpand, cuyo comportamiento se controla con las variables globales trigexpand (que le es homónima), trigexpandplus y trigexpandtimes, cuyos valores por defectos son, respectivamente, false, true y true. Puestos a experimentar, definamos primero una expresión trigonométrica e interpretamos algunos resultados,
64 5.1. TRANSFORMACIONES SIMBÓLICAS 61 (%i1) expr: x+sin(3*x+y)/sin(x); ( %o1) sin (y + 3 x) sin x + x (%i) trigexpand(expr); /* aquí trigexpand vale false */ ( %o) cos (3 x) sin y + sin (3 x) cos y sin x + x Vemos que sólo se desarrolló la suma del numerador; ahora cambiamos el valor lógico de trigexpand, (%i3) trigexpand(expr), trigexpand=true; ( cos 3 x 3 cos x sin x ) sin y + ( 3 cos x sin x sin 3 x ) cos y ( %o3) sin x + x Cuando la asignación de la variable global se hace como se acaba de indicar, sólo tiene efecto temporal durante esa ejecución, sin haberse alterado a nivel global. Podemos ordenar a Maxima que simplifique la expresión anterior, (%i4) ratsimp(%); ( %o4) ( 3 cos x sin x cos 3 x ) sin y + ( sin 3 x 3 cos x sin x ) cos y x sin x sin x Podemos inhibir el desarrollo de la suma, (%i5) trigexpand(expr), trigexpandplus=false; ( %o5) sin (y + 3 x) sin x + x Otra variable global de interés, que no tiene nada que ver con la función trigexpand, es halfangles, con valor false por defecto, que controla si se reducen los argumentos trigonométricos con denominador dos, (%i6) sin(x/); ( x ) ( %o6) sin (%i7) sin(x/), halfangles=true;
65 6 CAPÍTULO 5. ÁLGEBRA ( %o7) 1 cos x La función trigsimp fuerza el uso de las identidades fundamentales sin x + cos x = 1 y cosh x sinh x = 1 para simplificar expresiones, (%i8) 5*sin(x)^ + 4*cos(x)^; ( %o8) 5 sin x + 4 cos x (%i9) trigsimp(%); ( %o9) sin x + 4 Como un último ejemplo, veamos cómo reducir productos y potencias de funciones trigonométricas a combinaciones lineales, (%i10) -sin(x)^+3*cos(x)^*tan(x); ( %o10) 3 cos x tan x sin x (%i11) trigreduce(%); ( %o11) 3 sin ( x) + cos ( x) 1 Y si queremos transformar expresiones trigonométricas en complejas, (%i1) exponentialize(%); ( %o1) e i x + e i x 4 3 i ( e i x e i x) 1 4 Otras funciones útiles en este contexto y que el usuario podrá consultar en la documentación son: triginverses, trigsign y trigrat.
66 5.. ECUACIONES Ecuaciones Almacenar una ecuación en una variable es tan simple como hacer (%i1) ec: 3 * x = 1 + x; ( %o1) 3 x = x + 1 A partir de ahí, aquellas operaciones en las que intervenga la variable serán realizadas a ambos miembros de la igualdad; restamos x en los dos lados y a continuación dividimos lo que queda entre, (%i) %-x; ( %o) x = 1 (%i3) %/; ( %o3) x = 1 obteniendo de esta manera la solución de la ecuación como si hubiésemos operado manualmente. Ni qué decir tiene que la ecuación anterior se pudo haber resuelto de un modo más inmediato, (%i4) solve(ec); ( %o4) [ x = 1 ] La instrucción solve puede admitir como segundo argumento la incógnita que se pretende calcular, lo que resultará de utilidad cuando en la ecuación aparezcan constantes literales, (%i5) solve((-a)/x-3=b*x+1/x,x); ( %o5) [ (4 4 a) b x =, x = b ] (4 4 a) b b Las soluciones de las ecuaciones serán probablemente utilizadas en cálculos posteriores, por lo que nos interesará poder extraerlas de la lista anterior; a continuación tomamos el primer resultado calculado por Maxima mediante la función part y después asignamos a la variable sol el resultado numérico, (%i6) part(%,1);
67 64 CAPÍTULO 5. ÁLGEBRA ( %o6) x = (%i7) sol: rhs(%); ( %o7) (4 4 a) b b (4 4 a) b b La función rhs devuelve el miembro derecho de la igualdad, mientras que lhs haría lo propio con el miembro izquierdo. Es posible resolver ecuaciones polinómicas de grado 4, pero desgraciadamente, como es de esperar, Maxima no dispone de un método general algebraico que permita resolver ecuaciones polinómicas de grado mayor que cuatro, (%i8) solve(x^5-3*x^4 + *x^3 -*x^ - x + 4 = 0); ( %o8) [ 0 = x 5 3 x 4 + x 3 x x + 4 ] por lo que solve devolverá la misma ecuación sin resolver. También solve se puede utilizar para la resolución de sistemas, en cuyo caso las ecuaciones deben ir agrupadas en una lista, así como las incógnitas; nos planteamos la resolución del siguiente sistema no lineal de incógnitas x e y, { 3 x y = 6 x = y + a (%i9) solve([3*x^-y^=6,x=y+a],[x,y]); ( %o9) [[ x = [ x = 3 a +4+a, y = 3 a +4 a, y = ], ]] 3 a +4+3 a 3 a +4 3 a Una alternativa al uso de solve es la función algsys. Veamos cómo algsys trata la resolución de la ecuación polinómica anterior %o8, (%i10) algsys([x^5-3*x^4 + *x^3 -*x^ - x + 4 = 0],[x]); ( %o10) [[x = ], [x = i], [x = i ], [x = ], [x = ]] Como se ve, al no ser capaz de resolverla algebraicamente, nos brinda la oportunidad de conocer una aproximación numérica de la solución. La función algsys reconoce la variable global realonly, que cuando toma el valor true, hará que se ignoren las soluciones complejas,
68 5.. ECUACIONES 65 (%i11) realonly:true$ (%i1) %i10; /* recalcula la entrada %i10 */ ( %o1) [[x = ], [x = ], [x = ]] (%i13) realonly:false$ /* le devolvemos el valor por defecto */ Un sistema no lineal con coeficientes paramétricos y complejos 3u av = t +i u+t = 3v + u t u = 1 (%i14) algsys([3*u-a*v=t,(+%i)/(u+t)=3*v+u,t/u=1],[u,v,t]); ( %o14) [[ [ u = u = i a a+6 + a a+6, v = i+ a a+6, t = i a ] i+ a a+6, ]] a+6 + a a+6, v = i+ a a+6, t = i+ a a+6 Veamos cómo trata algsys las ecuaciones indeterminadas, devolviendo la solución en términos paramétricos, (%i15) algsys([3*x^-5*y=x],[x,y]); ( %o15) [[ x = %r 1, y = 3 %r 1 %r 1 5 ]] (%i16) %, %r1:1; ( %o16) [[ x = 1, y = ]] 5 Maxima nombra los parámetros siguiendo el esquema %rn, siendo n un número entero positivo. En la entrada %i16 pedimos que en %o15 se sustituya el parámetro por la unidad. En ocasiones conviene eliminar las variables de los sistemas de forma controlada; en tales casos la función a utilizar será eliminate, (%i17) eliminate([3*u-a*v=t,(+%i)/(u+t)=3*v+u,t/u=1],[u]);
69 66 ( %o17) CAPÍTULO 5. ÁLGEBRA [ t a v, ( a + 9 a ) v (5 a + 36) t v 4 t + 9 i + 18 ] Cuando de lo que se trata es de resolver un sistema de ecuaciones lineales, la mejor opción es linsolve, cuya sintaxis es similar a la de las funciones anteriores. En el siguiente ejemplo, el objetivo es x 4y + z = 1 3x + y + 9z = x + y 4x + y + z = 3y (%i18) linsolve( [ * x - 4 * y + * z = -, 1/3* x + * y + 9 * z = x + y, -4 * x + sqrt() * y + z = 3 * y], [x,y,z]); ( %o18) [ x = , y = , z = 53 ] Cuando la matriz de coeficientes del sistema es dispersa, la función fast_linsolve será preferible, ya que aprovechará tal circunstancia para encontrar las soluciones de forma más rápida. Si los coeficientes del sistema se tienen en formato matricial, quizás sea más apropiado el uso de la función linsolve_by_lu, tal como se indica en la sección dedicada a las matrices. No todo es resoluble simbólicamente. Existen en Maxima varios procedimientos cuya naturaleza es estrictamente numérica. Uno de ellos es realroots, especializado en el cálculo de aproximaciones racionales de las raíces reales de ecuaciones polinómicas; el segundo parámetro, que es opcional, indica la cota de error. (%i19) realroots(x^8+x^3+x+1, 5e-6); ( %o19) [ x = 1, x = ] 5488 En cambio, allroots obtiene aproximaciones en formato decimal de coma flotante de todas las raíces de las ecuaciones polinómicas, tanto reales como complejas, (%i0) allroots(x^8+x^3+x+1);
70 5.. ECUACIONES 67 ( %o0) [x = i ,, x = i, x = , x = i , x = i ,, x = i , x = i, x = ] Más allá de las ecuaciones algebraicas, find_root utiliza el método de bipartición para resolver ecuaciones en su más amplio sentido, (%i1) f(x):=144 * sin(x) + 1*sqrt(3)*%pi - 36*x^ - 1*%pi*x$ (%i) find_root(f(z),z,1,); ( %o) El uso del método de Newton requiere cargar en memoria el módulo correspondiente. Veamos como ejemplo la ecuación u u 5 = u (%i3) load(mnewton)$ /* carga el paquete */ (%i4) mnewton([*u^u-5],[u],[1]); 0 errors, 0 warnings ( %o14) [[u = ]] y el sistema { x + 3 log(x) y = 0 x xy 5x + 1 = 0 (%i5) mnewton([x+3*log(x)-y^, *x^-x*y-5*x+1],[x, y], [5, 5]); ( %o5) [[x = , y = ]] (%i6) mnewton([x+3*log(x)-y^, *x^-x*y-5*x+1],[x, y], [1, -]); ( %o6) [[x = , y = ]] En los anteriores ejemplos, el primer argumento es una lista con la ecuación o ecuaciones a resolver, las cuales se suponen igualadas a cero; el segundo argumento es la lista de variables y el último el valor inicial a partir del cual se generará el algoritmo y cuya elección determinará las diferentes soluciones del problema.
71 68 CAPÍTULO 5. ÁLGEBRA 5.3. Inecuaciones Las rutinas para resolver inecuaciones las guarda Maxima en el paquete de funciones fourier_elim, por lo que primero las debemos cargar en memoria, (%i1) load("fourier_elim")$ Ciertos mensajes de avisos que nos devuelve Maxima los podemos olvidar, ya que solo nos indican que algunas macros del sistema han tenido que ser redefinidas. Vamos a resolver una inecuación de primer grado, (%i) ine : (5*x-16)/6 + (x+8)/1 < (x+1)/3; ( %o) 5 x x < x La función clave a llamar se llama precisamente fourier_elim y la sintaxis es similar a la de solve; el primer argumento a pasarle es la lista de inecuaciones, de momento una sola, y el segundo la lista de incógnitas, también una sola, (%i3) fourier_elim([ine], [x]); ( %o3) [x < 4] Vamos a ver ahora en qué intervalos de la recta el polinomio x 4 + 5x 3 + 5x 5x 6 toma valores positivos, (%i4) p : x^4+5*x^3+5*x^-5*x-6 $ (%i5) fourier_elim([p > 0],[x]); ( %o5) [1 < x] [ < x, x < 1] [x < 3] Tal como se ha visto hasta ahora, la función fourier_elim nos devuelve los resultados como disyunciones de listas, interpretándose cada una de ellas como una conjunción lógica de proposiciones; así, el último resultado habría que interpretarlo como (1 < x) ( < x x < 1) (x < 3). El resultado anterior cambia sensiblemente cuando admitimos la igualdad en la inecuación,
72 5.3. INECUACIONES 69 (%i6) fourier_elim([p >= 0],[x]); ( %o6) [x = 3] [x = ] [x = 1] [x = 1] [1 < x] [ < x, x < 1] [x < 3] Abordamos ahora la resolución de un sistema de inecuaciones lineales con una incógnita, (%i7) fourier_elim( [(x+)/4 <= x/-3, (8-x)/3 < (1+x)/-1], [x]); ( %o7) [x = 14] [14 < x] Ahora con tres, (%i8) fourier_elim( [0 < x, x < 1, 0 < y+x*y, y-*z+x < 1, 0 < z, x+y+z < 4], [z,x,y] ); [( %o8)( max 0, y + x 1 ) ] < z, z < y x + 4, 0 < x, x < min (1, 3 y), 0 < y, y < 3 Otros ejemplos de inecuaciones no lineales, (%i9) fourier_elim([max(x,y) < 1, min(x,y) > -1],[x,y]); ( %o19) [ 1 < x, x < 1, 1 < y, y < 1] (%i10) fourier_elim([abs(x) + abs(x/) + abs(x/3) # 1],[x]); ( %o10) [x = 0] [ 0 < x, x < 6 ] [ 611 ] [ 11 < x, x < 0 x < 6 ] [ ] < x (%i11) fourier_elim([log(x) < log(a)],[x,a]); ( %o11) [0 < x, x < a, 0 < a]
73 70 CAPÍTULO 5. ÁLGEBRA 5.4. Matrices La definición de una matriz es extremadamente simple en Maxima, (%i1) m1: matrix([3,4,0],[6,0,-],[0,6,a]); ( %o1) a También es posible definir una matriz de forma interactiva tal y como muestra el siguiente ejemplo, (%i) entermatrix(,3); Row 1 Column 1: 4/7; Row 1 Column : 0; Row 1 Column 3: %pi; Row Column 1: sqrt(); Row Column : log(3); Row Column 3: -9; Matrix entered. ( %o) ( 4 ) 7 0 π log 3 9 Existe un tercer método para construir matrices que es útil cuando el elemento (i, j)-ésimo de la misma es función de su posición dentro de la matriz. A continuación, se fija en primer lugar la regla que permite definir un elemento cualquiera y luego en base a ella se construye una matriz de dimensiones 5 (%i3) a[i,j]:=i+j$ (%i4) genmatrix(a,,5); ( %o4) ( ) Obsérvese que el símbolo de asignación para el elemento genérico es :=. Podemos acceder a los diferentes elementos de la matriz haciendo referencia a sus subíndices, indicando primero la fila y después la columna:
74 5.4. MATRICES 71 (%i5) m1[3,1]; ( %o5) 0 Se puede extraer una submatriz con la función submatrix, teniendo en cuenta que los enteros que preceden al nombre de la matriz original son las filas a eliminar y los que se colocan detrás indican las columnas que no interesan; en el siguiente ejemplo, queremos la submatriz que nos queda de m1 después de extraer la primera fila y la segunda columna, (%i6) submatrix(1,m1,); ( %o6) ( ) 6 0 a Otro ejemplo es el siguiente, (%i7) submatrix(1,,m1,3); ( %o7) ( 0 6 ) en el que eliminamos las dos primeras filas y la última columna, se pilla el truco? Al igual que se extraen submatrices, es posible añadir filas y columnas a una matriz dada; por ejemplo, (%i8) addrow(m1,[1,1,1],[,,]); ( %o8) a (%i9) addcol(%,[7,7,7,7,7]); ( %o9) a La matriz identidad es más fácil construirla mediante la función ident, (%i10) ident(3);
75 7 CAPÍTULO 5. ÁLGEBRA ( %o10) y la matriz con todos sus elementos iguales a cero, (%i11) zeromatrix(,4); ( %o11) ( ) También, una matriz diagonal con todos los elementos de la diagonal principal iguales puede construirse con una llamada a la función diagmatrix, (%i1) diagmatrix(4,%e); ( %o1) e e e e En todo caso, debe tenerse cuidado en que si la matriz no se construye de forma apropiada, Maxima no la reconoce como tal. Para saber si una expresión es reconocida como una matriz se utiliza la función matrixp; la siguiente secuencia permite aclarar lo que se pretende decir, (%i13) matrix([[1,,3],[4,5,6]]); /* construcción correcta */ ( %o13) ( [1,, 3] [4, 5, 6] ) (%i14) matrixp(%); /* es la anterior realmente una matriz? */ ( %o14) true (%i15) [[7,8,9],[0,1,]]; /* otra matriz */ ( %o15) [[7, 8, 9], [0, 1, ]] (%i16) matrixp(%); /* será una matriz? */ ( %o16) false Casos particulares de submatrices son las filas y las columnas; los ejemplos se explican por sí solos:
76 5.4. MATRICES 73 (%i17) col(m1,3); ( %o17) 0 a (%i18) row(m1,); ( %o18) ( 6 0 ) Con las matrices se pueden realizar múltiples operaciones. Empezamos por el cálculo de la potencia de una matriz: (%i19) m1^^; ( %o19) a 36 6 a a 1 Nótese que se utiliza dos veces el símbolo ^ antes del exponente; en caso de escribirlo una sola vez se calcularían las potencias de cada uno de los elementos de la matriz independientemente, como se indica en el siguiente ejemplo, (%i0) m:m1^; ( %o38) a Para la suma, resta y producto matriciales se utilizan los operadores +, - y., respectivamente, (%i1) m1+m; ( %o1) a + a (%i) m1-m; ( %o) a a
77 74 CAPÍTULO 5. ÁLGEBRA (%i3) m1.m; ( %o3) a a a Sin embargo, tanto el producto elemento a elemento de dos matrices, como la multiplicación por un escalar se realizan mediante el operador *, como indican los siguientes dos ejemplos, (%i4) m1*m; ( %o4) a 3 (%i5) 4*m1; ( %o5) a Otros cálculos frecuentes con matrices son la transposición, el determinante, la inversión, el polinomio característico, así como los valores y vectores propios; para todos ellos hay funciones en Maxima: (%i6) transpose(m1); ( %o6) (%i7) determinant(m1); /*la transpuesta*/ a /*el determinante*/ ( %o7) 36 4 a (%i8) invert(m1); /*la inversa*/ ( %o8) a 4 a 36 4 a 8 6 a 36 4 a 3 a a 36 4 a 36 4 a a a a (%i9) invert(m1),detout; /*la inversa, con el determinante fuera*/
78 5.4. MATRICES 75 ( %o9) (%i30) charpoly(m1,x); /*pol. 1 4 a 8 6 a 3 a a caract. con variable x*/ ( %o30) (3 x) (1 (a x) x) 4 (a x) (%i31) expand(%); /*pol. caract. expandido*/ ( %o31) x 3 + a x + 3 x 3 a x + 1 x 4 a + 36 Vamos a suponer ahora que a vale cero y calculemos los valores propios de la matriz, (%i3) m1,a=0; ( %o3) (%i33) eigenvalues(%); ( %o33) [[ ] ] 15 i i 3,, 6, [1, 1, 1] El resultado que se obtiene es una lista formada por dos sublistas, en la primera se encuentran los valores propios, que en este caso son λ 1 = 3 15 i, λ = i y λ 3 = 6, mientras que en la segunda sublista se nos indican las multiplicidades de cada una de las λ i. Para el cálculo de los vectores propios, (%i34) eigenvectors(%o3); ( %o34) [[[ [ 1, ] 15 i 3, 6 15 i+9 8, i i+3, ], [1, 1, 1] ], [ 1,, 15 i 9 8, i+1 8 ], [ 1, 3 4, 3 ] ] 4 Lo que se obtiene es, en primer lugar, los valores propios junto con sus multiplicidades, el mismo resultado que se obtuvo con la función eigenvalues, y a continuación los vectores propios de la matriz asociados a cada uno de los valores propios. A veces interesa que los vectores sean unitarios, de norma 1, para lo que será de utilidad la función uniteigenvectors, que se encuentra definida en el paquete eigen.lisp, lo que significa que antes de hacer uso de ella habrá que ejecutar la orden load(eigen). También podemos solicitar los vectores propios unitarios por medio de la función uniteigenvectors,
79 76 CAPÍTULO 5. ÁLGEBRA (%i35) uniteigenvectors(%o3); ( %o35) [[[ ] ] 15 i+3, 15 i 3, 6, [1, 1, 1], 3, 15 i+9, i 1 ], ], [ [ 8 3 3, 15 i , i [ 17, 3, ]] El rango, el menor de un elemento y una base del espacio nulo de cierta matriz pueden calcularse, respectivamente, con las funciones rank, minor y nullspace, (%i36) rank(%o3); /* el rango de la matriz*/ ( %o36) 3 (%i37) minor(%o3,,1); /* el menor de un elemento */ ( %o37) ( ) (%i38) nullspace(%o3); /* base del núcleo */ 0 errors, 0 warnings ( %o38) span () que es la respuesta que se obtiene cuando el espacio nulo está formado por un único elemento. De forma similar a como la instrucción map aplica una función a todos los elementos de una lista, matrixmap hace lo propio con los elementos de una matriz, (%i39) matrixmap(sin,%o3); ( %o39) sin 3 sin 4 0 sin 6 0 sin 0 sin 6 0 Avancemos un poco más en otros aspectos del álgebra lineal. Empezamos con el cálculo de la matriz de Jordan J de cierta matriz dada A, esto es, la matriz que verifica A = SJS 1, para cierta S. Para ello es necesario cargar previamente el paquete diag, se introduce la matriz A y la función jordan se encargará de obtener lo que buscamos en un formato que luego utilizamos para obtener las matrices J y S. Finalmente, comprobamos los resultados:
80 5.4. MATRICES 77 (%i40) A: matrix([,4,-6,0],[4,6,-3,-4],[0,0,4,0],[0,4,-6,]); ( %o40) (%i41) load(diag)$ jl: jordan(a); ( %o4) [[6, 1], [, ], [4, 1]] (%i43) J: dispjordan(jl); ( %o43) (%i44) S: ModeMatrix (A,jl); ( %o44) (%i45) is(a = S.J.(S^^-1)); ( %o45) true Cuando los coeficientes de un sistema lineal se tengan en forma matricial y los términos independientes en un vector, el sistema se podrá resolver con la función linsolve_by_lu, (%i46) A : matrix([sqrt(),4],[a,4/5]); ( %o46) ( ) 4 4 a 5 (%i47) B : [1/3,a];
81 78 ( %o47) [ ] 1 3, a CAPÍTULO 5. ÁLGEBRA Resolvemos ahora AX = B y luego simplificamos algo el resultado (%i48) linsolve_by_lu(a,b); ( %o48) (%i49) ratsimp(%); ( ) 4 a a a a a a, false ( %o49) 15 a 15 a 6 (15 5) a 60 a 1, false El resultado que se obtiene es una lista que contiene la solución y el símbolo false. Cuando la resolución del sistema sea de índole numérica, esta segunda componente será sustituida por el número de condición, una medida de la bondad del resultado numérico. Veamos de aplicar esto a cierta ecuación matricial del estilo MX = MM (%i50) M : matrix([sqrt(3),/5],[3,sin()]); ( %o50) ( ) sin (%i51) linsolve_by_lu(m,m.m, floatfield); ( %o51) [( ) ] , Veamos ahora cómo obtener algunas matrices especiales: la hessiana de una función, la de Hilbert, la de Toeplitz, la de Vandermone y la de Hankel, (%i5) hessian(exp(x^*y+a*z),[x,y,z]); ( %o5) 4 x y ea z+x y + y ea z+x y x3 y ea z+x y + x ea z+x y a x y ea z+x y x 3 y e a z+x y + x e a z+x y x 4 e a z+x y a x e a z+x y a x y e a z+x y a x e a z+x y a e a z+x y
82 5.4. MATRICES 79 (%i53) hilbert_matrix(5); ( %o53) (%i54) toeplitz([1,,3,4],[t,x,y,z]); ( %o54) 1 x y z 1 x y 3 1 x (%i55) vandermonde_matrix([u,v,x,y,z]); 0 errors, 0 warnings ( %o55) 1 u u u 3 u 4 1 v v v 3 v 4 1 x x x 3 x 4 1 y y y 3 y 4 1 z z z 3 z 4 (%i56) hankel ([v,x,y,z],[p,q,r,s]); ( %o56) v x y z x y z q y z q r z q r s El lector interesado puede obtener información sobre otras funciones matriciales accediendo a la documentación sobre el paquete linearalgebra. Veamos ahora cómo trabajar con matrices por bloques. Puesto que los elementos de una matriz pueden ser objetos de cualquier tipo, definamos los elementos como matrices, (%i57) m: matrix([matrix([1,0],[0,-1]),matrix([0,1],[1,0])], [matrix([0,1],[1,0]),matrix([1,0],[0,-1])]) ;
83 80 ( %o57) ( ) ( ) ( ) ( ) CAPÍTULO 5. ÁLGEBRA Podemos ahora hacer algunas operaciones básicas con esta matriz, (%i58) m + m ; ( %o58) ( ) ( ) 0 0 ( 0 0 ) ( ) (%i59) m - m ; ( %o59) ( ) ( ) ( ) ( ) (%i60) 3 * m ; ( %o60) ( ) ( ) ( ) ( ) Veamos lo que pasa con la multiplicación matricial, (%i61) m. m ; ( %o61) ( ) ( ) ( ) ( ) Este resultado es incorrecto porque en la multiplicación matricial, los productos entre elementos se hacen con el operador *, que es el correcto cuando los elementos son números; pero cuando los elementos son también matrices, la multiplicación de elementos debe realizarse con el operador del producto matricial. En el ejemplo anterior, los bloques se multiplicaron elemento a elemento. Para indicarle a Maxima la operación correcta, debemos hacer uso de la variable global matrix_element_mult,
84 5.5. PATRONES Y REGLAS 81 (%i6) matrix_element_mult : "." $ (%i63) m. m ; ( %o63) ( ) ( ) ( ) ( ) Patrones y reglas En un programa de cálculo simbólico, éste no sólo debe tener información sobre una función u operación a partir de su definición, sino que también habrá propiedades y reglas de transformación de expresiones de las que un programa como Maxima debe tener noticia. Nos planteamos la siguiente situación: necesitamos trabajar con una función G(x, y) que, independientemente de que esté definida o no, sabemos que es igual a la expresión H(x,y) x en todo su dominio, siendo H(x, y) otra función; queremos que la primera expresión sea sustituída por la segunda y además queremos tener bajo control estas sustituciones. Todo ello se consigue trabajando con patrones y reglas. Antes de definir la regla de sustitución es necesario saber a qué patrones será aplicable, para lo cual admitiremos que los argumentos de G(x, y) pueden tener cualquier forma: (%i1) matchdeclare ([x,y], true); ( %o1) done En este caso, las variables patrón serán x e y, siendo el segundo argumento una función de predicado 1 que devolverá true si el patrón se cumple; en el caso presente, el patrón es universal y admite cualquier formato para estas variables. Si quisiésemos que la regla se aplicase sólo a números enteros, se debería escribir matchdeclare ([x,y], integerp); si quisiésemos que la regla se aplicase siempre que x sea un número, sin imponer restricciones a y, escribiríamos matchdeclare (x, numberp, y, true). Como se ve, los argumentos impares son variables o listas de variables y los pares las condiciones de los patrones. Se define ahora la regla de sustitución indicada más arriba, a la que llamaremos regla1, 1 Se denominan así todas aquellas funciones que devuelven como resultado de su evaluación true o false.
85 8 CAPÍTULO 5. ÁLGEBRA (%i) defrule (regla1, G(x,y), H(x,y)/x); ( %o) regla 1 : G (x, y) H (x, y) x Aplicamos ahora la regla a las expresiones G(f, +k) y G(G(4, 6), +k), (%i3) apply1(g(f,+k),regla1); ( %o3) H (f, k + ) f (%i4) apply1(g(g(4,6),+k),regla1); ( %o4) 4 H ( ) H(4,6) 4, k + H (4, 6) Como se ve en este último ejemplo, la regla se aplica a todas las subexpresiones que contengan a la función G. Además, apply1 aplica la regla desde fuera hacia adentro. La variable global maxapplydepth indica a Maxima hasta qué nivel puede bajar en la expresión para aplicar la regla; su valor por defecto es 10000, pero se lo podemos cambiar, (%i5) maxapplydepth; ( %o5) (%i6) maxapplydepth:1; ( %o6) 1 (%i7) apply1(g(g(4,6),+k),regla1); ( %o7) H (G (4, 6), k + ) G (4, 6) Quizás sea nuestro deseo realizar la sustitución desde dentro hacia fuera, controlando también a qué niveles se aplica; applyb1 y maxapplyheight son las claves ahora. (%i8) maxapplyheight:1; ( %o8) 1
86 5.5. PATRONES Y REGLAS 83 (%i9) applyb1(g(g(4,6),+k),regla1); ( ) H (4, 6) ( %o9) G, k + 4 Obsérvese que hemos estado controlando el comportamiento de las funciones G y H sin haberlas definido explícitamente, pero nada nos impide hacerlo, (%i10) H(u,v):= u^v+1; v (%o10) H(u, v) := u + 1 (%i11) applyb1(g(g(4,6),+k),regla1); ( ) 4097 ( %o11) G 4, k + Continuemos esta exposición con un ejemplo algo más sofisticado. Supongamos que cierta función F verifica la igualdad F (x 1 + x x n ) = F (x 1 ) + F (x ) F (x n ) + x 1 x... x n y que queremos definir una regla que realice esta transformación. Puesto que la regla se aplicará sólo cuando el argumento de F sea un sumando, necesitamos una función de predicado que defina este patrón, (%i1) esunasuma(expr):= not atom(expr) and op(expr)="+" $ (%i13) matchdeclare(z,esunasuma)$ En la definición de la nueva regla, le indicamos a Maxima que cada vez que se encuentre con la función F y que ésta tenga como argumento una suma, la transforme según se indica: map(f, args(z)) aplica la función a cada uno de los sumandos de z, sumando después estos resultados con la función apply, finalmente a este resultado se le suma el producto de los sumandos de z, (%i14) defrule(regla, F(z), apply("+",map(f, args(z))) + apply("*", args(z))); ( %o14) regla : F (z) apply (+, map (F, args (z))) + apply (*, args (z)) Veamos unos cuantos resultados de la aplicación de esta regla, (%i15) apply1(f(a+b), regla);
87 84 CAPÍTULO 5. ÁLGEBRA ( %o15) F (b) + a b + F (a) (%i16) apply1(f(a+b), regla) - apply1(f(a+c), regla) ; ( %o16) F (c) a c + F (b) + a b (%i17) apply1(f(a+b+c), regla); ( %o17) F (c) + a b c + F (b) + F (a) (%i18) apply1(f(4), regla); ( %o18) F (4) (%i19) apply1(f(4+5), regla); ( %o19) F (9) (%i0) simp:false$ /* inhibe la simplificación de Maxima */ (%i1) apply1(f(4+5), regla); ( %o1) (F (4) + F (5)) (%i) simp:true$ /* restaura la simplificación de Maxima */ (%i3) %o1; ( %o3) F (5) + F (4) + 0 En los ejemplos recién vistos, hemos tenido que indicar expresamente a Maxima en qué momento debe aplicar una regla. Otro mecanismo de definición de reglas es el aportado por tellsimp y tellsimpafter; el primero define reglas a aplicar antes de que Maxima aplique las suyas propias, y el segundo para reglas que se aplican después de las habituales de Maxima. Otra función útil en estos contextos es declare, con la que se pueden declarar propiedades algebraicas a nuevos operadores. A continuación se declara una operación de nombre o como infija (esto es, que sus operandos se colocan a ambos lados del operador); obsérvese que sólo cuando la operación se declara como conmutativa Maxima considera que a o b es lo mismo que b o a, además, no ha sido necesario definir qué es lo que hace la nueva operación con sus argumentos, (%i4) infix("o");
88 5.5. PATRONES Y REGLAS 85 ( %o4) o (%i5) is(a o b = b o a); ( %o5) false (%i6) declare("o", commutative)$ (%i7) is(a o b = b o a); ( %o7) true
89 86 CAPÍTULO 5. ÁLGEBRA
90 Capítulo 6 Cálculo 6.1. Funciones matemáticas En Maxima están definidas un gran número de funciones, algunas de las cuales se presentan en la Figura 6.1. A continuación se desarrollan algunos ejemplos sobre su uso. Empecemos por algunos ejemplos, (%i1) log(%e); ( %o1) 1 (%i) abs(-3^-x); ( %o) (%i3) signum(-3^-x); 1 3 x ( %o3) 1 La función genfact(m,n,p) es el factorial generalizado, de forma que genfact(m,m,1) coincide con m!, (%i4) genfact(5,5,1)-5!; ( %o4) 0 y genfact(m,m/,) es igual a m!!, (%i5) genfact(5,5/,)-5!!; 87
91 88 CAPÍTULO 6. CÁLCULO abs(x) abs(x) min(x1,x,...) mín(x 1, x,...) max(x1,x,...) máx(x 1, x,...) 1 si x < 0 signum(x) signo(x) = 0 si x = 0 1 si x > 0 x! x! x!! x!! ( binomial(m,n) m n) = m(m 1)...[m (n 1)] n! genfact(m,n,p) m(m p)(m p)... [m (n 1)p] sqrt(x) x exp(x) e x log(x) ln(x) sin(x) sin(x) cos(x) cos(x) tan(x) tan(x) csc(x) csc(x) sec(x) sec(x) cot(x) cot(x) asin(x) arcsin(x) acos(x) arc cos(x) atan(x) arctan(x) atan(x,y) arctan( x y ) ( π, π) sinh(x) cosh(x) sinh(x) = 1 (ex e x ) cosh(x) = 1 (ex + e x ) tanh(x) tanh(x) = sinh(x) cosh(x) asinh(x) arcsinh(x) acosh(x) arccosh(x) atanh(x) arctanh(x) gamma(x) Γ(x) = 0 e u u x 1 du, x > 0 gamma_incomplete(a,x) Γ(a, x) = x e t t a 1 dt beta(a,b) B(a, b) = Γ(a)Γ(b) Γ(a+b) beta_incomplete(a,b,x) B(a, b, x) = x 0 (1 t)b 1 t a 1 dt erf(x) erf(x) = x 0 π e u du Figura 6.1: Algunas funciones de Maxima.
92 6.1. FUNCIONES MATEMÁTICAS 89 ( %o5) 0 Maxima siempre devuelve resultados exactos, que nosotros podemos solicitar en formato decimal, (%i6) asin(1); ( %o6) (%i7) float(%); π ( %o7) Pero si damos el argumento en formato decimal, Maxima también devuelve el resultado en este mismo formato, sin necesidad de hacer uso de float, (%i8) asin(1.0); ( %o8) Recordemos que el formato decimal lo podemos pedir con precisión arbitraria, (%i9) fpprec:50$ bfloat(%o8); ( %o9) B 10 0 La función de error está relacionada con la función de distribución de la variable aleatoria normal X N (0, 1) de la forma Φ(x) = Pr(X x) = erf ( x ), por lo que la probabilidad de que esta variable tome un valor menor que 1.5 es (%i10) *erf(1.5/sqrt()),numer; ( %o10) Una forma más elegante de hacer lo anterior es definir nuestra propia función de distribución a partir de la de error, para lo que se hace uso del símbolo :=,
93 90 CAPÍTULO 6. CÁLCULO (%i11) F(x):=1/+erf(x/sqrt())/ $ (%i1) F(1.5),numer; ( %o1) No terminan aquí las funciones de Maxima; junto a las ya expuestas habría que incluir las funciones de Airy, elípticas y de Bessel, sobre las que se podrá obtener más información ejecutando la instrucción?? y utilizando como argumento airy, elliptic o bessel, según el caso. Por ejemplo, (%i13)?? airy 0: airy_ai (Funciones y variables para las funciones especiales) 1: airy_bi (Funciones y variables para las funciones especiales) : airy_dai (Funciones y variables para las funciones especiales) 3: airy_dbi (Funciones y variables para las funciones especiales) Enter space-separated numbers, all or none : 0 -- Función: airy_ai (<x>) Función Ai de Airy, tal como la definen Abramowitz y Stegun, Handbook of Mathematical Functions, Sección La ecuación de Airy diff (y(x), x, ) - x y(x) = 0 tiene dos soluciones linealmente independientes, y = Ai(x) y y = Bi(x). La derivada diff (airy_ai(x), x) es airy_dai(x). Si el argumento x es un número decimal real o complejo, se devolverá el valor numérico de airy_ai siempre que sea posible. Véanse airy_bi, airy_dai y airy_dbi. Maxima reconoce los dominios en el plano complejo de las funciones; los siguientes ejemplos lo demuestran: (%i14) sin(%i); ( %o14) i sinh 1 (%i15) log(-3.0); ( %o15) i (%i16) asin(.0); ( %o16) i
94 6.. LÍMITES Límites Sin más preámbulos, veamos algunos ejemplos de cómo calcular límites con la asistencia de Maxima. En primer lugar vemos que es posible hacer que la variable se aproxime al infinito (x ) haciendo uso del símbolo inf, o que se aproxime al menos infinito (x ) haciendo uso de minf, (%i1) limit(1/sqrt(x),x,inf); ( %o1) 0 (%i) limit((exp(x)-exp(-x))/(exp(x)+exp(-x)),x,minf); ( %o) 1 que nos permite calcular lím x 1 x y lím x e x e x e x +e x, respectivamente. Los siguientes ejemplos muestran límites en los que la variable x se aproxima a puntos de discontinuidad, (%i3) limit((x^-x/3-/3)/(5*x^-11*x+6),x,1); ( %o3) 5 3 (%i4) limit(1/(x-1)^,x,1); ( %o4) donde hemos obtenido los resultados y lím x 1 x x 3 3 5x 11x + 6 = lím x 1 (x 1) =. Sin embargo, ciertos límites no se pueden resolver sin aportar información 1 adicional, tal es el caso de lím x 1 x 1, para el que podemos hacer (%i5) limit(1/(x-1),x,1); ( %o5) infinity donde Maxima nos responde con el símbolo und de undefined o indefinido. En tales situaciones podemos indicarle al asistente que la variable x se aproxima a 1 por la derecha (x 1 + ) o por la izquierda (x 1 ), (%i6) limit(1/(x-1),x,1,plus); ( %o6) (%i7) limit(1/(x-1),x,1,minus); ( %o7)
95 9 CAPÍTULO 6. CÁLCULO 6.3. Derivadas Maxima controla el cálculo de derivadas mediante la instrucción diff. A continuación se presentan algunos ejemplos sobre su uso, (%i1) diff(x^log(a*x),x); /* primera derivada */ ( %o1) x log(a x) ( log (a x) x + log x ) x (%i) diff(x^log(a*x),x,); /*derivada segunda*/ ( %o) x log(a x) ( log (a x) x + log x ) ( +x log(a x) log (a x) x x log x x + ) x (%i3) factor(%); /* ayudamos a Maxima a mejorar el resultado */ ( %o3) x log(a x) ( log (a x) + log x log (a x) log (a x) + log x log x + ) donde se han calculado la primera y segunda derivadas de la función y = x ln(ax). Nótese que en la entrada (%i3) se le ha pedido al programa que factorizase la salida (%o). Pedimos ahora a Maxima que nos calcule el siguiente resultado que implica derivadas parciales, ϑ 10 ϑx 3 ϑy 5 ϑz (ex sin(y) tan(z)) = e x cos(y) sec (z) tan(z). (%i4) diff(exp(x)*sin(y)*tan(z),x,3,y,5,z,); ( %o4) e x cos y sec z tan z Maxima también nos puede ayudar a la hora de aplicar la regla de la cadena en el cálculo de derivadas de funciones vectoriales con variable también vectorial. Supónganse que cierta variable z depende de otras dos x y y, las cuales a su vez dependen de u y v. Veamos cómo se aplica la regla de la cadena para obtener ϑz ϑv, ϑz ϑyϑv o ϑz ϑuϑv. (%i5) depends(z,[x,y],[x,y],[u,v]); ( %o5) [z (x, y), x (u, v), y (u, v)] (%i6) diff(z,v,1);
96 6.3. DERIVADAS 93 ( %o6) ( ) d d d v y d y z + d ( ) d d v x d x z (%i7) diff(z,y,1,v,1); ( %o7) ( ) d d d v y d y z + d ( ) d d v x d x d y z (%i8) diff(z,u,1,v,1); ( %o8) ( ( d d u y d d v y d ) z d y ) ( ( d d u x d d v x d z d x + d d v x ( d d x d y z )) + d d u d v y ( d d y z ) + + d d v y ( d d x d y z )) + d d u d v x ( d d x z) La función depends admite un número par de argumentos, los cuales pueden ser variables o listas de variables. Las que ocupan los lugares impares son las variables dependientes y las que están en las posiciones pares las independientes; denominación ésta que tiene aquí un sentido relativo, ya que una misma variable, como la x de la entrada %i5, puede actuar como independiente respecto de z y como dependiente respecto de u y v. En cualquier momento podemos solicitarle a Maxima que nos recuerde el cuadro de dependencias, (%i9) dependencies; ( %o9) [z (x, y), x (u, v), y (u, v)] También podemos eliminar algunas de la dependencias previamente especificadas, (%i10) remove(x,dependency); ( %o10) done (%i11) dependencies; ( %o11) [z (x, y), y (u, v)] (%i1) diff(z,y,1,v,1); ( %o1) ( ) d d d v y d y z Veamos cómo deriva Maxima funciones definidas implícitamente. En el siguiente ejemplo, para evitar que y sea considerada una constante, le declararemos una dependencia respecto de x,
97 94 CAPÍTULO 6. CÁLCULO (%i13) depends(y,x)$ (%i14) diff(x^+y^3=*x*y,x); ( ) ( ) d d ( %o14) 3 y d x y + x = x d x y + y Cuando se solicita el cálculo de una derivada sin especificar la variable respecto de la cual se deriva, Maxima utilizará el símbolo del para representar las diferenciales, (%i15) diff(x^); ( %o15) x del (x) lo que se interpretará como xdx. Si en la expresión a derivar hay más de una variable, habrá diferenciales para todas, (%i16) diff(x^+y^3=*x*y); ( %o16) ( ( ) ) ( ( ) ) d d 3 y d x y + x del (x) = x d x y + y del (x) Recuérdese que durante este cálculo estaba todavía activa la dependencia declarada en la entrada (%i13). Finalmente, para acabar esta sección, hagamos referencia al desarrollo de Taylor de tercer grado de la función en el entorno de x = 1, y = x ln x x 1 (%i17) taylor((x*log(x))/(x^-1),x,1,3); ( %o17) (%i18) expand(%); 1 (x 1) (x 1) ( %o18) x 3 1 x x A continuación un ejemplo de desarrollo multivariante de la función y = exp(x sin(xy)) alrededor del punto (, 0) hasta grado respecto de cada variable,
98 6.4. INTEGRALES 95 (%i19) taylor(exp(x^*sin(x*y)),[x,,],[y,0,]); ( %o19) 1+8 y+3 y + + ( 1 y + 96 y + ) (x )+ ( 6 y + 10 y + ) (x ) + (%i0) expand(%); ( %o0) 10 x y 384 x y + 30 y + 6 x y 1 x y + 8 y + 1 En ocasiones resulta necesario operar con funciones cuyas derivadas son desconocidas, quizás porque la propia función lo es. Esta situación puede llevar a que Maxima devuelva expresiones realmente complicadas de manipular. Un ejemplo es el siguiente, en el que trabajamos con una función f arbitraria (%i1) taylor(f(x + x^),x,1,1); ( d ( %o1) f () + d x f ( x + x ) ) (x 1) + x=1 Podemos facilitar las cosas si le definimos una función derivada a f, a la que llamaremos df, (%i) gradef(f(x),df(x))$ (%i3) taylor(f(x+x^),x,1,1); ( %o3) f () + 3 df () (x 1) + El paquete pdiff, que se encuentra en la carpeta share/contrib/, aporta una solución alternativa a este problema Integrales La función de Maxima que controla el cálculo de integrales es integrate, tanto para las definidas como indefinidas; empecemos por estas últimas, (%i1) integrate(cos(x)^3/sin(x)^4,x); ( %o1) 3 sin x 1 3 sin 3 x (%i) integrate(a[3]*x^3+a[]*x^+a[1]*x+a[0],x);
99 96 CAPÍTULO 6. CÁLCULO ( %o) a 3 x a x a 1 x que nos devuelve, respectivamente, las integrales + a 0 x y cos 3 x sin 4 x dx (a 3 x 3 + a x + a 1 x + a 0 )dx. Además, este último ejemplo nos ofrece la oportunidad de ver cómo escribir coeficientes con subíndices. Ahora un par de ejemplos sobre la integral definida, (%i3) integrate(*x/((x-1)*(x+)),x,3,5); ( log 7 ( %o3) 3 log log 4 3 log ) 3 (%i4) float(%); /*aproximación decimal*/ ( %o4) (%i5) integrate(asin(x),x,0,u); Is u positive, negative, or zero? positive; ( %o5) u arcsin u + 1 u 1 esto es, y u x dx (x 1)(x + ) arcsin(x)dx = u arcsin(u) + 1 u 1, u > 0. Nótese en este último ejemplo cómo antes de dar el resultado Maxima pregunta si u es positivo, negativo o nulo; tras contestarle escribiendo positive; (punto y coma incluido) obtenemos finalmente el resultado. En previsión de situaciones como esta, podemos darle al sistema toda la información relevante sobre los parámetros utilizados antes de pedirle el cálculo de la integral, (%i6) assume(u>0);
100 6.4. INTEGRALES 97 ( %o6) [u > 0] (%i7) integrate(asin(x),x,0,u); ( %o7) u arcsin u + 1 u 1 Cuando Maxima no puede resolver la integral, siempre queda el recurso de los métodos numéricos. El paquete quadpack, escrito inicialmente en Fortran y portado a Lisp para Maxima, es el encargado de estos menesteres; dispone de varias funciones, pero nos detendremos tan sólo en dos de ellas, siendo la primera la utilizada para integrales definidas en intervalos acotados, (%i8) /* El integrador simbólico no puede con esta integral */ integrate(exp(sin(x)),x,,7); ( %o8) 7 e sin x dx (%i9) /* Resolvemos numéricamente */ quad_qag(exp(sin(x)),x,,7,3); ( %o9) [ , , 31, 0 ] La función quad_qag tiene un argumento extra, que debe ser un número entero entre 1 y 6, el cual hace referencia al algoritmo que la función debe utilizar para la cuadratura; la regla heurística a seguir por el usuario es dar un número tanto más alto cuanto más oscile la función en el intervalo de integración. El resultado que obtenemos es una lista con cuatro elementos: el valor aproximado de la integral, la estimación del error, el número de veces que se tuvo que evaluar el integrando y, finalmente, un código de error que será cero si no surgieron problemas. La otra función de integración numérica a la que hacemos referencia es quad_qagi, a utilizar en intervalos no acotados. En el siguiente ejemplo se pretende calcular la probabilidad de que una variable aleatoria χ de n = 4 grados de libertad, sea mayor que la unidad (Pr(χ 4 > 1)), (%i10) n:4$ (%i11) integrate(x^(n/-1)*exp(-y/)/^(n/)*gamma(n/),x,1,inf); Integral is divergent -- an error. To debug this try debugmode(true); (%i1) quad_qagi(x^(n/-1)*exp(-x/)/^(n/)*gamma(n/),x,1,inf);
101 98 ( %o1) CAPÍTULO 6. CÁLCULO [ , , 165, 0 ] El integrador simbólico falla emitiendo un mensaje sobre la divergencia de la integral. La función quad_qagi ha necesitado 165 evaluaciones del integrando para alcanzar una estimación numérica de la integral, la cual se corresponde aceptablemente con la estimación que hacen los algoritmos del paquete de distribuciones de probabilidad (ver Sección 8.1). Otras funciones de cuadratura numérica son quad_qags, quad_qawc, quad_qawf, quad_qawo y quad_qaws, cuyas peculiaridades podrá consultar el lector interesado en el manual de referencia. La transformada de Laplace de una función f(x) se define como la integral L(p) = 0 f(x)e px dx, siendo p un número complejo. Así, la transformada de Laplace de f(x) = ke kx es (%i13) laplace(k*exp(-k*x),x,p); ( %o13) k p + k y calculando la transformada inversa volvemos al punto de partida, (%i14) ilt(%,p,x); ( %o14) k e k x La transformada de Fourier de una función se reduce a la de Laplace cuando el argumento p toma el valor it, siendo i la unidad imaginaria y t R, F (t) = 0 f(x)e itx dx. De esta manera, la transformada de Fourier de f(x) = ke kx es (%i15) laplace(k*exp(-k*x),x,-%i*t); ( %o15) k k i t Nótese que si x > 0, la f(x) anterior es precisamente la función de densidad de una variable aleatoria exponencial de parámetro k, por lo que este último resultado coincide precisamente con la función característica de esta misma distribución.
102 6.5. ECUACIONES DIFERENCIALES 99 Puesto que hablamos de la transformada de Fourier, quizás sea este un buen lugar para mencionar la transformada discreta de Fourier. Su cálculo se hace con el algoritmo de la transformada rápida, implementado en la función fft, que admite como argumento una lista o array de datos cuya longitud debe ser una potencia de dos. La función inversa está programada en la función inverse_fft. Para más información, consúltese la documentación correspondiente a estas funciones. Como comentario final, queda hacer notar que la función integrate admite integrandos matriciales, tal como indica el siguiente ejemplo. (%i16) M: matrix([sin(x), x^+1],[cos(x), exp(x)]); ( %o16) ( sin x x ) + 1 cos x (%i17) integrate(m, x) = integrate(m, x); ( %o17) ( sin x x ) ( + 1 cos x x 3 cos x e x dx = 3 + x ) sin x e x 6.5. Ecuaciones diferenciales Con Maxima se pueden resolver analíticamente algunas ecuaciones diferenciales ordinarias de primer y segundo orden mediante la instrucción ode. Una ecuación diferencial de primer orden tiene la forma general F (x, y, y ) = 0, donde y = dy dx. Para expresar una de estas ecuaciones se hace uso de diff, (%i1) /* ecuación de variables separadas */ ec:(x-1)*y^3+(y-1)*x^3* diff(y,x)=0; ( ) d ( %o1) x 3 (y 1) d x y + (x 1) y 3 = 0 siendo obligatorio el uso de la comilla simple ( ) antes de diff al objeto de evitar el cálculo de la derivada, que por otro lado daría cero al no haberse declarado la variable y como dependiente de x. Para la resolución de esta ecuación tan solo habrá que hacer (%i) ode(ec,y,x); ( %o) e x y 1 y = %c x 1 x donde %c representa una constante, que se ajustará de acuerdo a la condición inicial que se le imponga a la ecuación. Supóngase que se sabe que cuando x =, debe verificarse que y = 3, lo cual haremos saber a Maxima a través de la función ic1,
103 100 CAPÍTULO 6. CÁLCULO (%i3) ic1(%o,x=,y=-3); ( %o3) y 1 y = x + 7 x 36 7 x Veamos ejemplos de otros tipos de ecuaciones diferenciales que puede resolver Maxima, (%i4) /* ecuacion homogénea */ ode(x^3+y^3+3*x*y^* diff(y,x),y,x); ( %o4) 4 x y 3 + x 4 4 = %c En este caso, cuando no se incluye el símbolo de igualdad, se da por hecho que la expresión es igual a cero. (%i5) /* reducible a homogénea */ ode( diff(y,x)=(x+y-1)/(x-y-1),y,x); ( %o5) log ( y + x x + 1 ) ( ) + arctan x 1 y = %c 4 (%i6) /* ecuación exacta */ ode((4*x^3+8*y)+(8*x-4*y^3)* diff(y,x),y,x); ( %o6) y x y + x 4 = %c (%i7) /* Bernoulli */ ode( diff(y,x)-y+sqrt(y),y,x); ( %o7) log ( y 1) = x + %c (%i8) solve(%,y); ( %o8) ] [y = e x+ %c + e x + %c + 1 En este último caso, optamos por obtener la solución en su forma explícita. Una ecuación diferencial ordinaria de segundo orden tiene la forma general F (x, y, y, y ) = 0, siendo y la segunda derivada de y respecto de x. Como ejemplo, (%i9) diff(y,x)=x+ diff(y,x,);
104 6.5. ECUACIONES DIFERENCIALES 101 ( %o9) d d x y = d d x y + x (%i10) ode(%,y,x); ( %o10) y = %k 1 e x + x + x + + %k Maxima nos devuelve un resultado que depende de dos parámetros, %k1 y %k, que para ajustarlos necesitaremos proporcionar ciertas condiciones iniciales; si sabemos que cuando x = 1 entonces y = 1 y y = dy dx =, x=1 haremos uso de la instrucción ic, (%i11) ic(%,x=1,y=-1,diff(y,x)=); ( %o11) y = x + x + En el caso de las ecuaciones de segundo orden, también es posible ajustar los parámetros de la solución especificando condiciones de contorno, esto es, fijando dos puntos del plano por los que pase la solución; así, si la solución obtenida en (%o10) debe pasar por los puntos ( 1, 3) y (, 5 3 ), hacemos (%i1) bc(%o10,x=-1,y=3,x=,y=5/3); 7 35 ex+1 ( %o1) y = 6 e x + x e e 3 6 Nótese que este cálculo se le solicita a Maxima con bc. La resolución de sistemas de ecuaciones diferenciales se hace con llamadas a la función desolve. En este contexto es preciso tener en cuenta que se debe utilizar notación funcional dentro de la expresión diff; un ejemplo aclarará este punto, resolviendo el sistema { df(x) dx = 3f(x) g(x) dg(x) dx = f(x) g(x) (%i13) desolve([ diff(f(x),x)=3*f(x)-*g(x), diff(g(x),x)=*f(x)-*g(x)], [f(x),g(x)]); ( %o13) [f (x) = g (x) = ( g(0) f(0)) e x ( g(0) 4 f(0)) e x 3 3, ] (4 g(0) f(0)) e x (g(0) f(0)) e x 3 3
105 10 CAPÍTULO 6. CÁLCULO Como se ve, las referecias a las funciones deben incluir la variable independiente y las ecuaciones estarán acotadas entre corchetes, así como los nombres de las funciones. Observamos en la respuesta que nos da Maxima la presencia de f(0) y g(0), lo cual es debido a que se desconocen las condiciones de contorno del sistema. En este último ejemplo, supongamos que queremos resolver el sistema de ecuaciones diferenciales df(x) dx = f(x) + g(x) + 3h(x) dg(x) dx = g(x) h(x) dh(x) dx = f(x) + h(x) bajo las condiciones f(0) = 1, g(0) = 3 y f(0) = 1. En primer lugar introduciremos estas condiciones con la función atvalue, para posteriormente solicitar la resolución del sistema, (%i14) atvalue(f(x),x=0,-1)$ (%i15) atvalue(g(x),x=0,3)$ (%i16) atvalue(h(x),x=0,1)$ (%i17) desolve([ diff(f(x),x)=f(x)+g(x)+3*h(x), diff(g(x),x)=g(x)-*h(x), diff(h(x),x)=f(x)+h(x)], [f(x),g(x),h(x)]); ( %o17) [ f (x) = x e x + e x e x, g (x) = x e x + e x + e x, h (x) = x e x + e x] La función desolve también nos va a permitir resolver ecuaciones de orden mayor que dos. En el siguiente ejemplo, abordamos la resolución de la ecuación d 3 y(x) dx 3 + d y(x) dx + dy(x) dx + y(x) = 0 bajo las condiciones que se describen a continuación, (%i18) atvalue( diff(y(x),x,),x=0,v)$ (%i19) atvalue( diff(y(x),x),x=0,u)$ (%i0) atvalue(y(x),x=0,w)$ Ahora resolvemos, (%i1) desolve( diff(y(x),x,3)+ diff(y(x),x,)+ diff(y(x),x)+y(x)=0, y(x)); ( %o1) y (x) = (w + v + u) sin x + (w v) cos x + (w + v) e x El paquete plotdf permite generar campos de direcciones, bien de ecuaciones diferenciales de primer orden dy = F (x, y), dx
106 6.5. ECUACIONES DIFERENCIALES 103 bien de sistemas { dx dt dy dt = G(x, y) = F (x, y) Los argumentos a pasar a la función plotdf son la función F, en el primer caso, y una lista con las funciones F y G en el segundo. Las variables serán siempre x e y. Como ejemplo, pidamos a Maxima que genere el campo de direcciones de la ecuación diferencial dy dx = 1 + y + y (%i) load(plotdf)$ (%i3) plotdf(1 + y + y^); El gráfico que resulta es el de la Figura 6. a), en el que además se observan dos trayectorias que se dibujaron de forma interactiva al hacer clic sobre dos puntos del plano. La función plotdf admite varias opciones, algunas de las cuales aprovechamos en el siguiente ejemplo. Supongamos el modelo predador-presa de Lotka-Volterra, dependiente de dos parámetros h y k, { dx dt dy dt = x + hxy = x + kxy El siguiente código permite generar el campo de direcciones correspondiente, dándoles inicialmente a h y k los valores -1. y 0.9, respectivamente; además, aparecerán sobre la ventana gráfica dos barras de deslizamiento para alterar de forma interactiva estos dos parámetros y ver los cambios que se producen en el campo. (%i4) plotdf([*x+k*x*y, -y+h*x*y], [parameters,"k=-1.,h=0.9"], [sliders,"k=-:,h=-:"]); En la Figura 6. b) se observa el gráfico obtenido despúes de pedir trayectorias concretas que pasan por varios puntos del plano (en el archivo gráfico no aparecen las barras de deslizamiento). Cuando Maxima no es capaz de resolver la ecuación propuesta, se podrá recurrir al método numérico de Runge-Kutta, el cual se encuentra programado en el paquete diffeq. Como primer ejemplo, nos planteamos la resolución de la ecuación dy dt = y + exp( 3t), con la condición y(0) = 1. La función ode es incapaz de resolverla: (%i5) ec: diff(y,t)+*y^-exp(-3*t)=0; ( %o5) d d t y + y e 3 t = 0
107 104 CAPÍTULO 6. CÁLCULO a) b) Figura 6.: Campos de direcciones creados con la función plotdf : a) campo de la ecuación dy dx = 1+y+y ; b) campo correspondiente al modelo predadorpresa.
108 6.5. ECUACIONES DIFERENCIALES 105 (%i6) ode(ec,y,t); ( %o6) false Abordamos ahora el problema con un enfoque numérico, para lo cual definimos la expresión en Maxima, f(t, y) = dy dt = y + exp( 3t) (%i7) load(diffeq)$ (%i8) f(t,y):= -*y^+exp(-3*t) $ (%i9) res: runge1(f,0,5,0.5,1); ( %o9) [[0.0, 0.5, 1.0, 1.5,.0,.5, 3.0, 3.5, 4.0, 4.5, 5.0], [1.0, , , , , , , , , , ], [ 1.0, , , , , , , , , , ]] (%i30) drawd( points_joined = true, point_type = dot, points(res[1],res[]), terminal = eps)$ La función runge1 necesita cinco argumentos: la función derivada dy dt, los valores inicial, t 0, y final, t 1, de la variable independiente, la amplitud de los subintervalos y el valor que toma y en t 0. El resultado es una lista que a su vez contiene tres listas: las abscisas t, las ordenadas y y las correspondientes derivadas. Al final del ejemplo anterior, se solicita la representación gráfica de la solución, cuyo aspecto es el mostrado por la Figura 6.3 a). (Consúltese la sección 7 para una descripción de la función drawd.) Nos planteamos ahora la resolución de la ecuación diferencial de segundo orden d y dt = 0.(1 y ) dy dt y, con y(0) = 0.1 y y (0) = 1 para lo cual definimos en Maxima la función g(t, y, y ) = d y dt = 0.(1 y )y y,
109 106 CAPÍTULO 6. CÁLCULO a) b) c) d) Figura 6.3: Resolución numérica de ecuaciones diferenciales con diffeq : a), solución de la ecuación de primer orden dy dt = y + exp( 3t), y(0) = 1; b), solución de la ecuación de segundo orden d y = 0.(1 y ) dy dt dt y, con las condiciones y(0) = 0.1 y y (0) = 1; c), diagrama (y, y ); d), diagrama (y, y ). (%i31) g(t,y,yp) := 0.*(1-y^)*yp - y $ (%i3) res: runge(g,0,100,0.1,0.1,1)$ (%i33) midibujo(i,j):= drawd(points_joined = true, (%i34) midibujo(1,)$ (%i35) midibujo(,3)$ (%i36) midibujo(,4)$ point_type = dot, points(res[i],res[j]), terminal = eps)$ La función runge necesita seis argumentos: la función derivada d y dt, los valores inicial, t 0, y final, t 1, de la variable independiente, la amplitud de los subintervalos y los valores que toman y y su primera derivada en t 0. El resultado es una lista que a su vez contiene cuatro listas: las abscisas t, las ordenadas y, las correspondientes primeras derivadas y, por último, las segundas derivadas. Al final de este último ejemplo se solicitan algunos gráficos asociados a la solución, los formados con los pares (t, y), (y, y )
110 6.6. VECTORES Y CAMPOS VECTORIALES 107 y (y, y ), que son los correspondientes a los apartados b), c) y d) de la Figura 6.3. El paquete dynamics dispone de otra rutina para el método de Runge- Kutta, rk, que permite la resolución de sistemas de ecuaciones difereciales. Para resolver el sistema { dx dt = 4 x 4y dy dt = y x + 1 cargamos el paquete y ejecutamos la sentencia correspondiente, junto con el gráfico asociado que no mostramos. (%i37) load(dynamics)$ (%i38) rk([4-x^-4*y^,y^-x^+1],[x,y],[-1.5,0.75],[t,0,4,0.0])$ (%i39) draw3d( points_joined = true, point_type = dot, points(%), terminal = eps)$ Para más información sobre la sintaxis de esta función, ejecútese? rk. El paquete dynamics contiene otras funciones relacionadas con los sistemas dinámicos y los fractales Vectores y campos vectoriales En Maxima, los vectores se introducen como simples listas (Sección 4.1), siendo el caso que con ellas se pueden realizar las operaciones de adición, producto por un número y producto escalar de vectores, (%i1) [1,,3]+[a,b,c]; ( %o1) [a + 1, b +, c + 3] (%i) s*[a,b,c]; ( %o) [a s, b s, c s] (%i3) [1,,3].[a,b,c]; /* producto escalar */ ( %o3) 3 c + b + a El cálculo del módulo de un vector se puede hacer mediante la definición previa de una función al efecto:
111 108 CAPÍTULO 6. CÁLCULO (%i4) modulo(v):= if listp(v) then sqrt(apply("+",v^)) else error("mucho ojito: ", v, " no es un vector!!!!")$ (%i5) xx:[a,b,c,d,e]$ (%i6) yy:[3,4,-6,0,4/5]$ (%i7) modulo(xx-yy); ( %o7) ( e 4 ) + d 5 + (c + 6) + (b 4) + (a 3) Los operadores diferenciales que son de uso común en el ámbito de los campos vectoriales están programados en el paquete vect, lo que implica que debe ser cargado en memoria antes de ser utilizado. Sigue a continuación una sesión de ejemplo sobre cómo usarlo. Partamos de los campos escalares φ(x, y, z) = x+y +z y ψ(x, y, z) = 4x + log(y + z ) y demostremos que sus superficies de nivel son ortogonales probando que φ ψ = 0, siendo el operador gradiente, (%i8) /* Se carga el paquete */ load(vect)$ (%i9) /* Se definen los campos escalares */ phi: y^+z^-x$ psi:log(y^+z^)+4*x$ (%i11) grad(phi). grad(psi); ( %o11) grad ( z + y x ) grad ( log ( z + y ) + 4 x ) Como se ve, Maxima se limita a devolvernos la misma expresión que le introducimos; el estilo de trabajo del paquete vect requiere el uso de dos funciones: express y ev, la primera para obtener la expresión anterior en términos de derivadas y la segunda para forzar el uso de éstas. (%i1) express(%); ( %o1) d d z d d y d d x (%i13) ev(%,diff); ( z + y x ) ( d d z ( z + y x ) ( d d y ( z + y x ) ( d d x ( log ( z + y ) + 4 x )) + ( ( log z + y ) + 4 x )) + ( ( log z + y ) + 4 x )) ( %o13) 4 z z + y + 4 y z + y 4
112 6.6. VECTORES Y CAMPOS VECTORIALES 109 (%i14) ratsimp(%); ( %o14) 0 Al final, hemos tenido que ayudar un poco a Maxima para que terminase de reducir la última expresión. Sea ahora el campo vectorial definido por F = xyi + x zj e x+y k y pidámosle a Maxima que calcule su divergencia, ( F), rotacional ( F)y laplaciano ( F) (%i15) F: [x*y,x^*z,exp(x+y)]$ (%i16) div (F); /* divergencia */ ( %o16) div ([ x y, x z, e y+x]) (%i17) express (%); ( %o17) d d y ( x z ) + d d z ey+x + d d x (x y) (%i18) ev (%, diff); ( %o18) y (%i19) curl (F); /* rotacional */ ( %o19) curl ([ x y, x z, e y+x]) (%i0) express (%); ( %o0) [ d d y ey+x d d z ( x z ), (%i1) ev (%, diff); d d z d d (x y) d x ey+x, d x ( x z ) d d y ] (x y) ( %o1) [ e y+x x, e y+x, x z x ] (%i) laplacian (F); /* laplaciano */ ( %o) laplacian ([ x y, x z, e y+x]) (%i3) express (%);
113 110 ( %o3) CAPÍTULO 6. CÁLCULO d d z [ x y, x z, e y+x] + d d y [ x y, x z, e y+x] + d d x [ x y, x z, e y+x] (%i4) ev (%, diff); ( %o4) [ 0, z, e y+x ] Nótese en todos los casos el usos de la secuencia express - ev. Si el usuario encuentra incómoda esta forma de operar, siempre podrá definir una función que automatice el proceso; por ejemplo, (%i5) milaplaciano(v):= ev(express(laplacian(v)), diff) $ (%i6) milaplaciano(f); ( %o6) [ 0, z, e y+x ] Por último, el paquete vect incluye también la definición del producto vectorial, al cual le asigna el operador ~, (%i7) [a, b, c] ~ [x, y, z]; ( %o7) [a, b, c] [x, y, z] (%i8) express(%); ( %o8) [b z c y, c x a z, a y b x]
114 Capítulo 7 Gráficos Maxima no está habilitado para realizar él mismo gráficos, por lo que necesitará de un programa externo que realice esta tarea. Nosotros nos limitaremos a ordenar qué tipo de gráfico queremos y Maxima se encargará de comunicárselo a la aplicación gráfica que esté activa en ese momento, que por defecto será Gnuplot. La otra herramienta gráfica es Openmath, un programa Tcl-tk que se distribuye conjuntamente con Maxima El módulo plot Las funciones plotd y plot3d son las que se utilizan por defecto. Veamos un ejemplo de cada una de ellas. (%i1) xy:[[10,.6], [0,.9], [30,1.1], [40,1.3], [50,1.4]]$ (%i) plotd([[discrete,xy], *%pi*sqrt(u/980)], [u,0,50], [style, [points,5,,6], [lines,1,1]], [legend,"datos experimentais","predicion da teoria"], [xlabel,"lonxitude do pendulo (cm)"], [ylabel,"periodo (s)"], [gnuplot_preamble, "set terminal postscript eps;set out plotd.eps "])$ Se ha definido en el ejemplo una lista de puntos empíricos a representar, junto con su función teórica, pidiendo su representación gráfica en el dominio [0, 50]. El resto de código hace referencia a diversas opciones: style, legend, xlabel, ylabel y gnuplot_preamble, esta última permite escribir código de Gnuplot para que éste lo ejecute; aquí se le pide que devuelva el gráfico en formato Postscript. El ejemplo siguiente genera una superficie explícita. Las dos salidas gráficas se representan en la Figura 7.1 (%i3) plot3d (^(-u^ + v^), [u, -3, 3], [v, -, ], [gnuplot_preamble, "set terminal postscript eps;set out plot3d.eps "])$ 111
115 11 CAPÍTULO 7. GRÁFICOS 1.6 Datos experimentais Predicion da teoria periodo (s) Lonxitude do pendulo (cm) ^(v^-u^) Figura 7.1: Gráficos generados por las funciones plotd y plot3d.
116 7.. EL MÓDULO DRAW 113 El control de las opciones gráficas se consigue manipulando la variable global plot_options, cuyo estado por defecto es (%i4) plot_options; (%o4) [[x, e+305, e+305], [y, e+305, e+305], [t, - 3, 3], [GRID, 30, 30], [VIEW_DIRECTION, 1, 1, 1], [COLOUR_Z, FALSE], [TRANSFORM_XY, FALSE], [RUN_VIEWER, TRUE], [PLOT_FORMAT, GNUPLOT], [GNUPLOT_TERM, DEFAULT], [GNUPLOT_OUT_FILE, FALSE], [NTICKS, 10], [ADAPT_DEPTH, 10], [GNUPLOT_PM3D, FALSE], [GNUPLOT_PREAMBLE, ], [GNUPLOT_CURVE_TITLES, [DEFAULT]], [GNUPLOT_CURVE_STYLES, [with lines 3, with lines 1, with lines, with lines 5, with lines 4, with lines 6, with lines 7]], [GNUPLOT_DEFAULT_TERM_COMMAND, ], [GNUPLOT_DUMB_TERM_COMMAND, set term dumb 79 ], [GNUPLOT_PS_TERM_COMMAND, set size 1.5, 1.5;set term postsc# ript eps enhanced color solid 4]] Para mayor información sobre el significado de cada uno de los elementos de esta lista, así como de las funciones plotd y plot3d, consúltese la documentación. Por defecto, Maxima invocará al programa Gnuplot para la generación de gráficos, pero quizás prefiramos el programa Openmath, que forma parte de la distribución de Maxima; en tal caso tendríamos que modificar previamente las opciones guardadas en plot_options y a continuación solicitar el gráfico deseado, como en este caso en el que se representa la función gamma y su inversa. (%i5) set_plot_option([plot_format, openmath])$ (%i6) plotd([gamma(x),1/gamma(x)],[x,-4.5,5],[y,-10,10])$ El resto de esta sección lo dedicaremos a hacer una somera descripción del paquete draw, un proyecto consistente en el desarrollo de un interfaz que permita aprovechar al máximo las habilidades gráficas de Gnuplot. 7.. El módulo draw Aquí, las funciones a utilizar son drawd, draw3d y draw, para escenas en d, 3d y para gráficos múltiples y animaciones, respectivamente. Puesto que se trata de un módulo adicional, será necesario cargarlo en memoria antes de hacer uso de él. Empecemos por un ejemplo comentado. (%i1) load(draw)$ (%i) drawd(
117 114 CAPÍTULO 7. GRÁFICOS key = "Cubic poly", explicit(%pi*x^3+sqrt()*x^+10,x,0,), color = blue, key = "Parametric curve", line_width = 3, nticks = 50, parametric(*cos(rrr)+3, rrr, rrr, 0, 6*%pi), line_type = dots, points_joined = true, point_type = diamant, point_size = 3, color = red, line_width = 1, key = "Empiric data", points(makelist(random(40.0),k,1,5)), title = "DRAWING CURVES", terminal = eps )$ Los argumentos de drawd se dividen en tres grupos: objetos gráficos (en el ejemplo, explicit, parametric y points), opciones locales (que afectan directamente a la representación de los objetos gráficos, como key, color, line_width, nticks, line_type, points_joined y line_width) y opciones globales (que hacen referencia a aspectos generales del gráfico, como title y terminal). Las opciones se indican como igualdades, escribiendo a la izquierda el nombre de la opción y a la derecha el valor que se le quiera asignar, a la vez que los objetos gráficos tienen el formato igual que las llamadas a funciones, esto es, el nombre del objeto a dibujar seguido, entre paréntesis, de los parámetros que definen al objeto. Todos estos argumentos se interpretan secuencialmente, de forma que al asignar un cierto valor a una opción local, ésta afectará a todos los objetos gráficos que le sigan. En este ejemplo, line_width comienza teniendo valor 1, que es el asignado por defecto, luego se le da el valor 3 para la representación de la curva paramétrica y finalmente se le devuelve el valor original antes de representar los segmentos que unen los puntos aleatoriamente generados. En cambio, las dos opciones globales, title y terminal, aunque se colocaron al final, podrían haberse ubicado en cualquier otro lugar. Siguiendo con los gráficos en dos dimensiones, el siguiente ejemplo muestra una escena en la que intervienen los objetos ellipse, image, label, vector y, ya lo conocemos, explicit. Antes de ejecutar esta instrucción es necesario leer el fichero gráfico gatos.xpm 1 (%i3) cats: read_xpm("gatos.xpm")$ (%i4) drawd( 1 El formato gráfico XPM es el único que puede leer Maxima.
118 7.. EL MÓDULO DRAW 115 terminal = eps, yrange = [-4,10], ellipse(5,3,8,6,0,360), image(cats,0,0,10,7), line_width =, head_length = 0.3, color = blue, label(["this is Francisco",-1,-0.5]), vector([-1,0],[,4]), color = green, vector([11,7],[-,-1]), label(["this is Manolita",11,8]), explicit(sin(x)-,x,-4,15) )$ Junto con los objetos gráficos introducidos en los ejemplos, cuyos resultados se pueden ver en los apartados a) y b) de la Figura 7., también existen polygon, rectangle, polar, implicit y geomap, este último para mapas cartográficos. Mostramos a continuación algunas escenas tridimensionales. En primer lugar, el valor absoluto de la función Γ de Euler, junto con sus líneas de contorno. (%i5) gamma(x,y):= block([re,im,g:gamma(x+%i*y)], re:realpart(g), im:imagpart(g), sqrt(re^+im^))$ (%i6) draw3d( zrange = [0,6], xu_grid = 50, yv_grid = 50, surface_hide = true, contour = surface, contour_levels = [0,0.5,6], /* de 0 a 6 en saltos de 0.5 */ color = cyan, terminal = eps, explicit(gamma(x,y),x,-4,4,y,-,))$ Una superficie paramétrica, que junto con la anterior escena, se pueden ver ambas en la Figura 7., apartados c) y d). (%i7) draw3d( terminal = eps, title = "Figure 8 - Klein bottle", rot_horizontal = 360,
119 116 CAPÍTULO 7. GRÁFICOS xrange = [-3.4,3.4], yrange = [-3.4,3.4], zrange = [-1.4,1.4], xtics = none, ytics = none, ztics = none, axis_3d = false, surface_hide = true, parametric_surface((+cos(u/)*sin(v)-sin(u/)*sin(*v))*cos(u), (+cos(u/)*sin(v)-sin(u/)*sin(*v))*sin(u), sin(u/)*sin(v) + cos(u/)*sin(*v), u, -%pi, 360*%pi/180-%pi, v, 0, *%pi) )$ En el gráfico siguiente se combinan varios objetos: una superficie explícita, dos curvas paramétricas y dos etiquetas, cuyo aspecto es el mostrado en el apartado e) de la Figura 7.. (%i8) draw3d( color = green, explicit(exp(sin(x)+cos(x^)),x,-3,3,y,-3,3), color = blue, parametric(cos(5*u)^,sin(7*u),u-,u,0,), color = brown, line_width =, parametric(t^,sin(t),+t,t,0,), surface_hide = true, title = "Surface & curves", color = red, label(["up",-,0,3]), label(["down",,0,-3]), rot_horizontal = 10, rot_vertical = 84, terminal = eps )$ En el siguiente ejemplo hacemos una proyección esférica del hemisferio sur, que se ve en el apartado f) de la Figura 7.. El paquete worldmap carga en memoria las coordenadas latitud longitud de las líneas fronterizas y costeras de todos los países del mundo, así como algunas funciones necesarias para su manipulación y procesamiento, (%i9) load(worldmap)$ (%i10) draw3d( surface_hide = true, rot_horizontal = 60, rot_vertical = 131,
120 7.. EL MÓDULO DRAW 117 color = cyan, parametric_surface( cos(phi)*cos(theta), cos(phi)*sin(theta), sin(phi), theta,-%pi,%pi, phi,-%pi/,%pi/), color = red, geomap([south_america,africa,australia], [spherical_projection,0,0,0,1]), color = blue, geomap([south_america,africa,australia], [cylindrical_projection,0,0,0,1,]), terminal = eps)$ Además de los objetos gráficos tridimensionales ya vistos, también se hayan definidos points, vector e implicit. También es posible generar múltiples gráficos en un mismo fichero o hacer animaciones en formato GIF o en la ventana gráfica. Para ver más ejemplos de gráficos generados con el paquete draw, se recomienda acceder a la dirección o consultar el sistema de ayuda de Maxima. Ya como colofón, un ejemplo de gráfico múltiple en el que se muestra también cómo dar sombreado a las superficies tridimensionales. El resultado en el apartado g) de la Figura 7.. (%i11) escena1: gr3d(surface_hide = true, enhanced3d = true, palette = gray, explicit(sin(sqrt(x^+y^)),x,-5,5,y,-5,5))$ (%i1) escena: gr3d(surface_hide = true, enhanced3d = true, palette = gray, user_preamble = "set pm3d map", explicit(sin(sqrt(x^+y^)),x,-5,5,y,-5,5))$ (%i13) draw( columns =, terminal = eps_color,
121 CAPÍTULO 7. GRÁFICOS DRAWING CURVES Cubic poly Parametric curve Empiric data 8 This is Manolita This is Francisco a) b) Figure 8 - Klein bottle c) d) Surface & curves 4 UP DOWN e) f) g) Figura 7.: Gráficos generados con el paquete draw: a) y b) con drawd; c), d), e) y f) con draw3d; g), un gráfico múltiple.
UNIVERSIDADE DA CORUÑA Departamento de Computación FACULTADE DE INFORMÁTICA PROYECTO FIN DE CARRERA INGENIERÍA TÉCNICA EN INFORMÁTICA DE GESTIÓN Herramienta de análisis de redes en un sistema de información
Eugenia Bahit Curso: Python para Principiantes www.eugeniabahit.com Eugenia Bahit 2012. Creative Commons Atribución-NoComercial 3.0 Comparte el conocimiento en libertad! Python para Principiantes de Eugenia

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 resolución