intro.html

============================================

Sistemas de Computación Algebraica

como Ambientes de Cálculo Científico

============================================

Por Héctor Hernández 29/09/04

>

Algebra Computacional ?

En el campo de la computación científica los métodos y herramientas de análisis numérico son tradicionalmente más comunes que los Sistemas de Computación Algebraica ( SCA ).

La expresión: " cálculo por computadora " generalmente se entiende como Computación Numérica

- FORTRAN, Basic, C, C++, Java => Precisión fija (Punto Flotante)

>

Los SCA son programas que pueden operar con símbolos que representan Objetos Matemáticos :

- Números (Enteros, racionales, reales, complejos...)

- Polinómios, Funciones Racionales, Sistemas de Ecuaciones,

- Grupos, Anillos, Algebras . . .

>

Estas herramientas trabajan de la forma:

Input : solve(problema);

Output : respuesta

>

Sistemas de Propósito Especial: (Creados para hacer calculos en un area específica)

FORM, GAP, CAMAL, SHEEP, STENSOR, LiE, KANT

Sistemas de Propósito General: (Especies de navajas suizas!)

Axiom, Derive, Reduce, Maple, MatLab, Mathematica, Macsyma, MuPAD

Los Sistemas de Computación Algebraica (SCA).

Propiedades de los SCA

Son programas Interactivos
Precisión arimética arbitraria
Cálculos exactos con expresiones matemáticas simbólicas
Manipulación de expresiones y subexpresiones
Aproximación analítica y numérica
Extensión a la programación

>

Sistemas de Computación Algebraica Modernos de Propósito General

Un SCA de propósito general es un ambiente completamente integrado de computación para la investigación y la educación:

Interfaz gráfica (worksheet) o ambiente interactivo:

* Procesador de texto, de fórmulas y de gráficas.

     * Con salidas en Latex, RTF, HTML, FORTRAN y C; o hyperlinks a otros documentos.
     * Manuales en línea.
     * Enlaces a otros programas y bibliotecas

Capacidades para cálculo numérico
Capacidades para visualización, con salidas gráficas en diferentes formatos: PostScript, GIF,JPG, . . .
Pensado para usuarios no especializados en computación.

>

Los más populares SCA :

Derive, Reduce , Mathematica , Axiom , Maple, MuPAD, Matlab

>

Principal Ventaja:

Enorme capacidad para realizar cálculos algebraicos largos y tediosos.

>

Por ejemplo, demostrar que la función:

es solución de la Ecuación Diferencial:

Puede tomarle a un PC un tiempo de CPU relativamente corto:

tiempo de cpu = 1,365 seg

>

Un Ejemplo: Maple

MapleV es una herramienta de computación científica con las siguientes caracteristicas generales:

° Manipulador Simbólico

° Gran colección de Funciones Numéricas

° Capacidad gráfica en 2D y 3D

° Lenguaje de programación Avanzada

° Sintaxis similar al FORTRAN, PASCAL o C

° Hoja de Cálculo y Editor de Texto con salidas en Latex y HTML

° Maple está dispononible para DOS, Windows, MacOS, UNIX, Linux...
° La gran ventaja : una hoja de cálculo puede ser utilizada sobre cualquier plataforma sin necesidad de ser alterada.

>

1 Sintaxis Básica:

1.- Zona de Texto (En color Negro)

2.- Zona de comandos (En Color Rojo) (Cada comando debe finalizar con ; o con : )

3.- Zona de respuestas (En color Azul)

> 25! + 13^23;

Las líneas con comentarios se colocan a continuación del símbolo #

> c^2 = a^2 + b^2 ; # Teorema de Pitágoras

El siguiente comando reinicia la hoja de cálculo

> restart;

>

2 Primeros Pasos

Maple hace cálculos tanto con numeros enteros como en Punto Flotante

> 15 + 5^20;

> 15. + 5^20;

Pero el énfasis radica en los cálculos exáctos

> cos(Pi/12)^2 + ln(2/3+5)/7;

> evalf(%); # EVALuate using Floating-point

Con % se puede utilizar la última salida de Maple como entrada en el comando siguiente.

Existe un orden de precedencia para los operadores: + - * / ^

> 1+2*3^2;

> (1+2)*3^2;

Por lo tanto, su uso descuidado produce errores

> 2^3^5;

Error, `^` unexpected

> (2^3)^5;

> 2^-5;

Error, `-` unexpected

> 2^(-5);

Constantes numéricas y letras griegas. Maple hace distinción entre mayúsculas y minúsculas

> evalf(Pi);

> evalf(pi);

La precisión de su evaluación puede ser controlada:

> sqrt( 68 ); sqrt( 68. ); evalf(sqrt(68),50 );

Se puede redondear a una fracción

> convert(%,fraction);

El simbolo % permite que la ultima salida se asigne a la nueva instrucción!

> restart;

Para asignar un objeto a una variable se hace con el operador :=

> f := arctan( (2*x^2-1)/(2*x^2+1) );

> f;

Mientras que el símbolo = , se utliza para mostrar una relación entre variables

> h = x + y;

> h;

Por lo tanto, existen nombres que ya estan definidos y por lo tanto no puedes ser utilizados

> abs:= f*sin(x);

Error, attempting to assign to `abs` which is protected

abs es la función valor absoluto

> abs(x-a);

> restart;

3 Cálculo Elemental

Se definen, se evalúan y se derivan funciones abstractas

> f:= (x,y) -> exp(x^2+y^2)/(x-y);

> f(a^(1/2),b^(1/2));

> f(2,3);

> evalf(%);

> diff( f(x,y) , x) + diff( f(x,y) , y) ;

Conserva las expresiones para la derivada interna

> f:= (t,r) -> g(c*t+r) + h(c*t-r);

> diff(f(t,r),t$3);

Para funciones compuestas se utiliza el operador @ y para derivar el operador D

> (cos@ln@sqrt)(x);

> D(cos@ln@sqrt)(x);

> diff(cos(ln(sqrt(x))),x);

Límites:

> f:= (x) -> 1/(1+r/x)^x;

> Limit(f(x),x=infinity) = limit(f(x), x=infinity);

Las mayúsculas a la izquierda ( Limit ) proveen una representación del Límite. Las minúsculas ( limit ) su ejecución.

También podemos calcular límites a la izquierda y a la derecha.

> g:= (x)-> tan(x + Pi) ;

> Limit(g(x), x=Pi/2,right)=limit(g(x),x=Pi/2,right);

Sumatorias:

> Sum((1+i)/i^2, i = 1..n) ;

> value(%);

Si se tienen dudas se invoca la ayuda

> ?gamma

> ?Psi

Integración:

> g:=(x)-> 1/( x*( b*x+c*x^2)^(1/2) );

> Int( g(x) , x ) = int( g(x) , x );

Evalua integrales definidas tanto analítica como numéricamente

> f:=(x)-> exp(Pi*x);

> Int( f(x) , x=1..3) = int( f(x) , x=1..3);

Evalua integrales impropias y las funciones que de ella emergen

> h:= (t)-> exp(-t)*t^(z-1);

> Int(h(t),t=0..infinity) = int(h(t),t=0..infinity);

> ?GAMMA

> GAMMA(0.5); # Funcion Gamma

> Int(x/(x^2+2*x+2),x) = int(x/(x^2+2*x+2),x);

Comprobemos que la solucion es la correcta

> diff(rhs(%),x);

> simplify(%);

> restart;

4 Resolución de sistemas de ecuaciones algebraicas

Resuelve sistemas de ecuaciones algebráicas. El resultado queda almacenado en una lista, y luego hay que asignar esos valores al "nombre" de las variables correspondientes para seguir operando.

> ecu1:= R[1]*i[1]+R[3]*i[3]+R[4]*i[4]-V[1]=0;

> ecu2:= R[2]*i[2]-V[2]-R[3]*i[3]=0;

> ecu3:= i[1]-i[2]-i[3]= 0;

> solve({ecu1,ecu2,ecu3},{i[1],i[2],i[3]} ); assign(%);

> i[1],i[2],i[3];

Ahora al evaluar x, y, z en las ecuaciones ecu1 ... ecu3

> ecu1; ecu2; ecu3;

> simplify(ecu1); simplify(ecu2); simplify(ecu3);

Seguidamente procedemos a "limpiar" el contenido de las variables para poder utilizarlas en otras expresiones.

> i[1]:='i[1]'; i[2]:='i[2]'; i[3]:='i[3]';

Soluciones aproximadas para encontrar las raíces de un polinomio

> p:= x^3 - 3*x^2 = 17*x - 51;

> sol:=solve(p,x);

> evalf(%);

> q:= x^5 +3*x^4 - 6*x^3 + x^2 + 12*x - 14=0;

> solve(q , x);

> ?RootOf

> fsolve(q,x); # fsolve utiliza punto-flotante

> restart:

5 Ecuaciones Diferenciales

> ED1:=diff(y(x),x)+y(x)/x=alpha/(x*y(x)^2);

> dsolve(ED1,y(x));

> sols:=dsolve({ED1,y(1)=5},y(x));

Resuelve ecuaciones diferenciales, con y sin valores iniciales por varios métodos. Series . . .

> dsolve({ED1,y(1)=5},y(x),series);

Sistemas de Ecuaciones Difereciales. Mediante transformaciones de Laplace. . .

> sys:= diff(y(x),x)=z(x) , diff(z(x),x)=y(x);

> fcns:={y(x),z(x)};

> s:=dsolve({sys,y(0)=0,z(0)=1},fcns,method=laplace);

> s[1]; # Selecciono una de las soluciones

Puedo tomar únicamente el lado derecho de la solución para una posterior manipulación

> 2*( rhs(s[1] )); # right hand sides

Ecuaciones diferenciales famosas: Bessel

> ED2:=x^2*diff(y(x),x$2)+x*diff(y(x),x)+(x^2-mu)*y(x)=0;

> dsolve(ED2,y(x)); assign(%);

> simplify(ED2);

Se pueden hacer graficas para un conjunto de condiciones iniciales

> y:=subs(mu=10, _C1=1/2, _C2=1/3 ,y(x));

> plot(y, x=0..20,-2..1 ,title=`y(x)`);

Soluciones Numéricas:

> de1 :={diff(g(t),t$2)=-f(t)-g(t),diff(f(t),t$2)= diff(g(t),t)+f(t)};

> init1 := {g(0)=1, D(g)(0)=0, f(0)=0, D(f)(0)=1}:

> F := dsolve(de1 union init1, {g(t),f(t)},type=numeric);

> F(0.0);

> F:=dsolve(de1 union init1,{g(t),f(t)},type=numeric,output=array([0,1,2,3,4]));

> with(plots):

Warning, the name changecoords has been redefined

> odeplot(F,[t,g(t)],1..10,labels=[t,g]);

> odeplot(F,[t,f(t)],1..10,labels=[t,f]);

> with(DEtools):

> S:=cos(xi)*diff(h(xi),xi$3)-diff(h(xi),xi$2)+Pi*diff(h(xi),xi)=h(xi)-xi;

> DEplot(S,h(xi),xi=-2.5..1.4,[[h(0)=1,D(h)(0)=2,(D@@2)(h)(0)=1]],h=-4..5,stepsize=.05);

> restart;

6 Manipulación algebraica

Maple , automáticamente evalua y simplifica expresiones

> p:=(a+b)^3*(a+b)^2;

> expand(p);

> factor(p);

> t:=cos(x)^5+sin(x)^4+2*cos(x)^2-2*sin(x)^2-cos(2*x);

> simplify(t);

> r:=alpha*(x^3-y^3)/(beta*(x^2+x-y-y^2));

> normal(r);

> n:=numer(r); d:=denom(r);

> n*d*delta;

> expand(%);

> collect(%,x);

> coeff(%,x^3);

> factor(%);

Puede convertir muchos tipos de expresiones en otras más especificas

> expr1:=(a*x^2+b)/(x*(-3*x^2-x+4));

> expr2:=convert(expr1,parfrac,x);

> expr3:=sin(x)*cos(x);

> expr4:=convert(%,exp);

Se pueden aislar los operandos de las expresiones

> op(3,expr2); op(3,expr4);

> op(3,expr2)*op(3,expr4);

> simplify(%);

No siempre todo funciona de la manera correcta. La siguiente intregal es trivial, pero veamos lo que sucede con Maple

> Int(2*x*(x^2+1)^24,x);

> value(%);

> factor(%);

> factor(%+1/25);

> restart;

7 Bibliotecas

No todos los comandos son cargados en la memoria cuando Maple V es iniciado. Sólo los comandos estandard son cargados automáticamente

> P:= x^5-24*x^4+85*x^3-108*x^2+166*x-120 ;

> factor(P);

> realroot(P);

> readlib(realroot):

> realroot(P,1);

Las Bibliotecas de Maple V :

- La biblioteca estandard

- La biblioteca de misceláneos

- Paquetes

> ?index

> ?index,packages

> ?index,misc

> restart:

8 Graficos y Visualización

Una mayor capacidad grafica es obtenida con el paquete plots

> with(plots);

Warning, the name changecoords has been redefined

Graficos básicos en 2 y 3 Dimensiones:

> f:=x-> exp(-x^2)*sin(Pi*x^3);

> plot( f(x),x=-2..2,title=`f(x)`);

>

> g:=(x,y)-> f(x)*f(y);

> plot3d(g(x,y),x=-1..1,y=-1..1);

>

> h:=(x) -> cos(sqrt(x^2+3*y^2))/(1+x^2/8);

> j:= (x) -> 1/3 -(2*x^2+y^2)/19;

> plot3d({h(x),j(x)},x=-3..3,y=-3..3,title=`Interseccion de h(x) y j(x)`);

Mapas Topol ó gicos

> with(plots):

> contourplot(sin(x*y),x=-10..10,y=-10..10);

Graficando Campos Vectoriales :

> phi[1]:=(x,y)->ln((x^2+2*x+1+y^2)^(1/2));

> phi[2]:=(x,y)->-ln((x^2-2*x+1+y^2)^(1/2));

> fieldplot([phi[1](x,y),phi[2](x,y)],x=-2..2,y=-2..2,arrows=SLIM,grid=[8,8]);

> gradplot((-phi[1](x,y)-phi[2](x,y)),x=-2..2,y=-1..1,arrows=SLIM,grid=[10,10]);

Animaciones:

> animate3d(cos(sqrt(t*x^2+t*y^2)),x=-5..5,y=-5..5,t=1..8);

Con un poco más de complejidad

> tubeplot({[10*cos(t),10*sin(t),0,t=0..2*Pi,radius=2+cos(7*t),numpoints=120,tubepoint=24],[0,10+5*cos(t),5*sin(t),t=0..2*Pi,radius=1.5,numpoint=50,tubepoint=18]});

>

9 Aplicaciones

Algebra Lineal

> restart:

> with(linalg):# Leemos el paquete de A.L.

Warning, the protected names norm and trace have been redefined and unprotected

>

Creando Matrices y Vectores

> A:=array( [ [a, b, c], [d, e, f ], [g, h, i ] ] );

> B:=array(antisymmetric,1..3,1..3):

> B[1,2]:=x+a; B[1,3]:=y+b; B[2,3]:=z+c;

> print(B);

> v:= vector([1,2,3] );

> C:=matrix(4,4, (i,j) -> i^(j-1));

> M:=array(antisymmetric,1..3,1..3):

> entermatrix(M);

enter element 1,2 > 2*alpha;

> 2*beta;

> 2*gamma:

> N:=diag(mu,nu,lambda);

> G:=diag(M,N);

>

La función "solve" para resolver un sistema de ecuaciones:

> ec1:= (c^2-1)*x+(c-1)*y = (1-c)^2 ;

> ec2:=(c-1)*x+(c^2-1)*y=c-1;

> sis:=({ec1,ec2}):

> solve(sis,{x,y});

${y = 2/c/(c+2), x = (c^2-2)/c/(c+2)}$

>

Algebra de matrices para resolver el sistema de ecuaciones

> H:=genmatrix(sis,[x,y],'flag');

> A:=genmatrix(sis,[x,y]);

> B:=delcols(H,1..2);

La ecuacion A X = B se puede resolver con:

> sol:= linsolve(A,B);

> x:=sol[1,1]; y:=sol[2,1];

> simplify(ec1); simplify(ec2);

> x:='x': y:='y': # Limpio las variables x y

>

Operaciones con matrices

> S:=evalm(M+N);

> evalm(3*M + 5/3*N);

> evalm(M*N);

Error, (in evalm/evaluate) use the &* operator for matrix/vector multiplication

> evalm(M &* N);

> evalm(N &* M);

> V:= vector([x,y,z]);

> evalm(M &* V);

> S1:= evalm( (M)^(3) );

> S1:= map(factor,S1);

> print(S);

> zeta:=inverse(S);

> map(simplify,evalm(S &* zeta));

> trace(zeta);

> det(zeta);

> normal( trace(zeta) / det(zeta) );

>

Cálculo de Autovectores y Autovalores

> K:=matrix([[sqrt(2),alpha,beta],[alpha,sqrt(2),alpha],[beta,alpha,sqrt(2)]]);

> pc:=charpoly(K, lambda); # Polinomio Caracteristico

> factor(%);

> solve(pc=0,lambda);

> eigenvals(K);

> eigenvects(K);

$[1/2*beta+2^(1/2)+1/2*(beta^2+8*alpha^2)^(1/2), 1, {vector([1/2*(1/2*beta+1/2*(beta^2+8*alpha^2)^(1/2))/alpha, 1, 1/2*(1/2*beta+1/2*(beta^2+8*alpha^2)^(1/2))/alpha])}], [1/2*beta+2^(1/2)-1/2*(beta^2+8*...$
$[1/2*beta+2^(1/2)+1/2*(beta^2+8*alpha^2)^(1/2), 1, {vector([1/2*(1/2*beta+1/2*(beta^2+8*alpha^2)^(1/2))/alpha, 1, 1/2*(1/2*beta+1/2*(beta^2+8*alpha^2)^(1/2))/alpha])}], [1/2*beta+2^(1/2)-1/2*(beta^2+8*...$

El resultado es una lista de la forma:

[ e_i, m_i, { v[1,i],... v[n_i,i] } ]

donde:

                e_i = son los autovalores
    m_i = multiplicidad
            {v[1,i], ..., v[n_i,i]} = conjunto de vectores bases

> print(K);

>

Solución de B x = v mediante Gauss-Jordan

> B:=matrix(3,3,[19,-50,88, 53,85,-49, 78, 17, 72] );

> v1:= vector( [3 ,5,-2] ):

> v2:= vector( [4 ,-5,9] ):

> v3:= vector( [21 ,4,7] ):

> augment(B,v1,v2,v3 );

Construimos la matriz aumentada, e invocamos la solución

> gaussjord(%);

> leastsqrs(B, v1); leastsqrs(B, v2);leastsqrs(B, v3);

Funciones test

> beta := sqrt(2)*sqrt(5-sqrt(5));

> A:= matrix([[(sqrt(5)+beta +1)/4, -beta/2 ],[ beta/4 ,( sqrt(5) - beta+1)/4]]);

> B:=matrix( [[(sqrt(5)+1)/4,-beta/4 ],[ beta/4,(sqrt(5)+1)/4 ] ]);

Queremos verificar si A y B son similares, es decir, si existe una matriz T, no singular, de manera que

B = T A T^(-1).

> issimilar(A,B,'T');

> map(simplify,T);

>

Cálculo Vectorial.

Gradiente

> f:= 4*x*y*z - 5*y*x^3;

> gradf:= grad(f,[x,y,z]);

Vectores

> v:=vector([4*x-3*x^3*y,7*x*y*z^2+5*y^3,4*x^2*y^2+2*x]);

Rotores

> curlv:= curl(v,[x,y,z] );

Jacobianos

> jacobian(%,[x,y,z]);

Laplacianos

> laplacf:= laplacian(f,[x,y,z] );

Evaluación boleana de expresiones

> evalb(laplacf = diverge(gradf,[x,y,z]));

Diferentes tipos de coordenadas

> f := r*sin(theta)*z^2: v := [r, theta, z]:

> gr:=grad(f, v, coords=cylindrical);

> diverge(gr,v, coords=spherical);

> laph= laplacian(f,v);

> laph= laplacian(f,v,coords=bipolarcylindrical);

> ?coords

>

10 Procesos en UNIX

Es posible utilizar Maple en un computador que funcione bajo alguno de los diferentes versiones de sistemas operativos tipo UNIX, como por ejemplo el Linux. Cuando no queremos utilizar el ambiente grafico, podemos recurrir al ambiente de texto escribiendo el comado maple

nodo5%> maple

| \ ^ / | Maple 9 (IBM INTEL LINUX)

<____ ____> Waterloo Maple Inc.
| Type ? for help.

> evalf(sqrt(Pi),100);

1.7724538509055160272981674833411451827975494561223871282138077898529112845\
91032181374950656738544665

> Int( tan(x), x ) = int( tan(x), x );

                             /
                           |
                           | tan(x) dx = -ln(cos(x))
                           |
                          /

>

Maple sobre UNIX utiliza los archivos de entradas y salidas estandard para leer e imprimir información: < , > , |

nodo5%> maple < archimap.txt

Con esta instrucción Maple ejecutará todos los comandos que se encuentran en el archivo archimap.txt

>

En la siguiente instrucción Maple ejecutará todos los comandos que se encuentran en el archivo archimap.txt y colocará los resultados en el archivo llamado archimap.out

nodo5%> maple < archimap.txt > archimap.out

>

También se puede hacer que todos los comandos del archivo sean ejecutados para luego ser enviados al terminal

nodo5%> maple < archimap.txt | more

>

Maple puede ser detenido temporalmente con el comando Control_Z (^Z), de manera que para colocar procesos en "background" se procede de la manera usual:

nodo5%> maple -q < archimap > archimap.out

^Z
Suspended
nodo5%> bg
[2] maple -q < archimap > archimap.out &

nodo5%>

>

Programa Mint:
Es un programa que permite verificar si un archivo fuente de Maple tiene errores de sintáxis

nodo5%> mint < archimap.txt

>

============================================