>

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

Maple V : un Ambiente de Cálculo Científico

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

28/07/99

>

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

Plataformas:

° Maple V está dispononible para DOS, Windows, MacOS, UNIX
° 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

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

3.- Zona de respuestas

> 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 V 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 V 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;

`^` unexpected

> 2^-5;

`-` unexpected

Constantes numéricas y letras griegas. MapleV 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);

> 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 simbolo = , se utliza para mostrar una relación entre variables

> h = x + y;

> h;

Por lo tanto, existen nombres que ya estan definidos

> 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(%);

> ?gamma

> ?Psi

Integración:

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

> Int( g(x) , x ): % = value(%);

Evalua integrales definidas tanto analítica como numéricamente

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

> Int( f(x) , x=1..3): %= value(%);

> Int( f(x) , x=1..3 ): %= evalf(%);

Evalua integrales impropias y las funciones que de ella emergen

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

> Int(h(t),t=0..infinity): %= value(%);

Definite integration: Can't determine if the integral is convergent.

Need to know the sign of --> z

Will now try indefinite integration and then take limits.

> assume(z>0):

> Int(h(t),t=0..infinity): %= value(%);

> GAMMA(0.5); # Funcion Gamma

> infolevel[int]:=5:

> Int(x/(x^2+2*x+2),x): %= value(%);

int/indef1:   first-stage indefinite integration

int/ratpoly:   rational function integration

> diff(rhs(%),x);# Chequeamos...

> 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);

> allvalues(%);# evalua expresiones con 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 únicamnte 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, method=mgear,
value=array([0,0.6,1.1,1.5,2.3,2.5]));

> with(plots):

> 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 Formatos de Salidas

MapleV permite pasar las expresiones a LATEX

> alias(f=f(x,y,z)):# Otra manera de definir funciones

> f:=sin(sqrt(x^2+y^2+z^2)/(y^2+z^2));

> Diff(f,x):%= diff(f,x);

> latex(%);

{\frac {\partial }{\partial x}}\sin({\frac {\sqrt {{x}^{2}+{y}^{2}+{z}

^{2}}}{{y}^{2}+{z}^{2}}})=\cos({\frac {\sqrt {{x}^{2}+{y}^{2}+{z}^{2}}

}{{y}^{2}+{z}^{2}}})x{\frac {1}{\sqrt {{x}^{2}+{y}^{2}+{z}^{2}}}}

\left ({y}^{2}+{z}^{2}\right )^{-1}

> latex(%,`c:/usuarios/hector/archi1.txt`);

...o en código FORTRAN (en precisión doble y optimizado si el caso)

> h:= Pi*f/2 + 3*diff(f,x)+ 5*diff(f,y) + 4*diff(f,z);

> precision := double:

> fortran([fun=h],optimized);

      t1 = x**2

      t2 = y**2

      t3 = z**2

      t5 = dsqrt(t1+t2+t3)

      t6 = t2+t3

      t7 = 1.D0/t6

      t8 = t5*t7

      t9 = dsin(t8)

      t12 = dcos(t8)

      t13 = 1.D0/t5

      t18 = t13*t7

      t20 = t6**2

      t22 = t5/t20

      fun = 0.3141592653589793D1*t9/2.D0+3.D0*t12*t13*t7*x+5.D0*t12*(t18

     #*y-2.D0*t22*y)+4.D0*t12*(t18*z-2.D0*t22*z)

> fortran([fun=h],optimized,filename=`c:/usuarios/hector/archi2.txt`);

> restart:

7 Manipulación algebraica

Maple V , 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 V

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

> value(%);

> factor(%);

> factor(%+1/25);

> restart;

8 Maple V como lenguaje de programación

> isprime(37);

> isprime(3778);

> for i to 8 do

> n := 2^(2*i+1)-1;

> if isprime(n)

> then print( n, `Es primo !`);

> else print( n, `No es primo`)

> fi;

> od;

>

> f := proc(x,y);

> 2*x-2*y;

> end;

> f(2,5); f(w,z+1);

>

> h := proc(x);

> if x > 5 then x^2 else 5-2*x fi;

> end;

> h(-2); h(20);

> restart:

9 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

- La biblioteca de programas desarrollados por usuarios (Share Library)

> ?index,library

> ?index,packages

> ?index,libmisc

> restart:

10 Graficos y Visualización

Una mayor capacidad grafica es obtenida con el paquete plots

> with(plots);

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]});

>

11 Procesos en UNIX

Es posible utilizar Maple V en UNIX por medio del comado: maple

nodo5%> maple

| \ ^ / | Maple V Release 3 (IBM De Venezuela)
._| \ | | / |_. Copyright (c) 1981-1994 by Waterloo Maple Software
\ MAPLE / University of Waterloo. All rights reserved.
<____ ____ > Waterloo Maple Software.
| Type ? for help.

>

> evalf(sqrt(Pi),100);

1.772453850905516027298167483341145182797549456122387\
12821380778985291128459132181374950656738544665

>

> Int( tan(x), x ) : "= value(");

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

>

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

nodo5%> maple < archimap

Maple V ejecutará todos los comandos que se encuentran en el archivo archimap.

>

nodo5%> maple < archimap > archimap.out

Maple V ejecutará todos los comandos que se encuentran en el archivo archimap y colocará los resultados en el archivo llamado archimap.out

>

nodo5%> maple < archimap | more

Aqui todos los comandos del archivo serán ejecutados para luego ser enviados al terminal

>

Maple V 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 V tiene errores de sintáxis

nodo5%> mint < archimap

>

12 Aplicaciones

Algebra Lineal

> restart:

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

Warning, new definition for norm

Warning, new definition for trace

>

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] );

Se puede leer datos de un archivo e introducirlos en una matriz

> readlib(readdata): # Libreria para leer datos

> mat:=matrix(readdata(`c:/usuarios/hector/datos`,3)):

> mat[1,3]; mat[3,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;

enter element 1,3 > 2*beta;

enter element 2,3 > 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});

>

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);

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);

> K1:=map(evalf,K);

> eigenvects(K1);

Error, (in linalg/evalf) matrix entries must all evaluate to float

> K2:=map(convert, K1 ,rational);

> eigenvects(K2);

> map(evalf,[%]);

>

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

Lectura de datos de un archivo para una posterior visualizacion

> Mdat:=matrix(readdata(`c:/usuarios/hector/setdatos`,5));

> x:=col(Mdat,3);

> y1:=col(Mdat,4); y2:=col(Mdat,5);

> lny1:=map(ln,y1); lny2:=map(ln,y2);

> numrows:=rowdim(Mdat);

Podemos hacer simples graficos de x vs y

> set1:= [seq( [x[i] , y1[i]] , i=1..numrows)]:

> set2:= [seq( [x[i] , y2[i]] , i=1..numrows)]:

> p1:=plot(set1,0..1,style=point,labels=[x,y],color=red):

> p2:=plot(set2,0..1,style=point,labels=[x,y],color=green):

> datos:=plots[display]({p1,p2}):

> datos;

Se quiere graficar de Q = x*ln(y1) + (1-x)*ln(y2) , lo siguiente: Q vs x

> Q:=[seq(x[i]*lny1[i]+(1-x[i])*lny2[i], i=1..numrows)]:

> Qsx:=[seq([ x[i], Q[i] ] , i=1..numrows)]:

> Qx:= plot({Qsx}, 0..1,0..0.1, style=point):

> Qfig:=plots[display]( {Qx } ):

> Qfig;

>

>

Astronomía

Datos Meteorológicos

Display of Viking Lander 2 meteorological data

Sols 51 to 53 (Summer 1976)

for Release 5

Maple procedures written by
Alexander F. Walz

alexander.f.walz@t-online.de
July 10, 1999

> restart:

> mean := proc(l::listlist)
local vals;
vals := op(map(x -> op(2, x), l));
(min(vals)+max(vals))/2
end:

> adapt := proc(l::listlist, scale::numeric, f::numeric)
op(map(x -> [op(1, x), op(2, x) * scale - f], l))
end:

> filein := `c:/usuarios/hector/segment1.htm`:

> data := NULL:

> alt := readline(filein):

> tbl := table([]):

> i := 0:

> while alt <> 0 do
if not(member(substring(convert(alt, name), 1 .. 1), {`C`, `<`})) then
i := i + 1:
tbl[i] := sscanf(alt, `%*f %*f %f %f %d %f %*f %f`)
fi:
alt := readline(filein):
od:

> data := [seq(tbl[j], j=1 .. i)]: i := 'i': j := 'j':

> temp := map(x -> [op(1, x), op(5, x)], data):

> p1 := plot(temp, sol=51 .. 53, `deg C`=-90 .. -10, axes=frame):

> wind := map(x -> [op(1, x), op(2, x)], data):

> p2 := plot([adapt(wind, -mean(temp)/mean(wind), -1.9*mean(temp) )],
sol=51 .. 53, axes=frame, color=navy, linestyle=3):

>

> plots[display](p1, p2, labels=[`sol`, `deg C`],
title=`Viking Lander 2 meteorological data Sols 51 to 53`);

>

Fractales

Nombre Fórmula de iteración Inicialización de z Inicialización de c
Julia zn+1 := zn + c z0 := x + I*y c = a + I*b

http://SunSite.informatik.rwth-aachen.de/maple/frame03.htm

Julia sets, written by John Oprea. oprea@math.csuohio.edu

> restart; with(plots):

> julia := proc(c,x, y)
local z, m;
z := evalf(x+y*I);
for m from 0 to 30 while abs(z) < 3 do
z := z^2 + c
od;
m
end:

> J := proc(d)
global phonyvar;
phonyvar := d;
(x, y) -> julia(phonyvar, x, y)
end:

> plot3d(0, -2 .. 2, -1.3 ..1.3, style=patchnogrid,
orientation=[-90,0], grid=[250, 250],
scaling=constrained, color=J(-1.25));

>

>

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