# Cálculos por temas

## Cálculos básicos

(2+%i)^5$

expand(%);

2/5+3/4;

%pi;

float(%pi);

bfloat(%pi);

fpprec:500;

bfloat(%pi);

2^2000;

set_display('ascii)$

bfloat(%pi);

2^2000;

antonio:binomial(49,7);

antonio;

b:3^15;

a/b;

float(%);

expand((1+%i)^10);

set_display('xml)$

integrate(sin(x),x);

f(x):=sin(x)+exp(x);

integrate(f(x),x);

diff(f(x),x);

diff(sin(x),x,2);

suma:0$
producto:1$
for j:1 thru 50 do (suma:suma+j,producto:producto*j)$
suma;
producto;

## Cálculos con matrices

Ejercicio 2.7 (encontramos algunas matrices que verifican la propiedad del ejercicio)

A:matrix([0,b],[1/b,1]);
A^^2-A-ident(2);

Ejercicio 2.8

A:matrix([a^2,a*b,a*c],[a*b,b^2,b*c],[a*c,c*b,c^2]);

A.A;

Ejercicio 2.9

A:matrix([1,2,3,4],[2,0,3,4],[2,3,0,4],[2,3,4,0]);

determinant(%);

invert(A);

echelon(A);

Ejercicio 2.13

A:matrix([1,2,3],[0,1,2],[0,0,1])$
PotAPorBinomio(n):=matrix([1,2*n,2*n^2+n],[0,1,2*n],[0,0,2*n+1])$
PotA(n):=block(
    [aux,j],
    aux:ident(3),
    for j:1 thru n do aux:aux.A,
    aux)$

PotA(7);
PotAPorBinomio(7);
A^^7;

## Sistemas de ecuaciones

## Espacios vectoriales

Un ejercicio sobre matrices de cambio de base

/*Ejercicio sobre cambios de base*/
/*Introducimos dos bases, beta y betap y calculamos las matrices de cambio de base entre ellas*/
basebeta:matrix([1,0,3,4],[1,3,4,2],[2,5,8,11],[4,6,8,10]);
rank(basebeta);
basebetap:matrix([2,1,3,1],[1,4,5,1],[7,6,5,4],[1,2,1,1]);
rank(basebetap);
mbcb:transpose(basebeta);
mbcbp:transpose(basebetap);
mbbp:invert(mbcb).mbcbp;
mbpb:invert(mbbp);
print("La matriz de cambio de base de beta a betaprima es",mbbp);
print("La matriz de cambio de base de betap a beta es",mbpb);

/*Determinamos si el conjunto de vectores S={(1,1,1,1)_b,(1,1,1,1)_bp,(2,1,1,2)_b,(4,3,2,1)_bp} es o no
una base. Para ello los expresamos todos en la base beta*/

v:transpose(mbbp.transpose([1,1,1,1]));

w:transpose(mbbp.[4,3,2,1]);

Senmatrizporfilas:matrix([1,1,1,1],v[1],[2,1,1,2],w[1]);

rank(%);

Ejercicio de cambio de base 4.26

mbcb1:transpose(matrix([0,1,0,0,0],[1,0,1,0,0],[0,0,0,1,2],[0,1,-1,1,0],[0,1,1,0,0]));

mbcb2:transpose(matrix([1,0,0,0,2],[-1,0,0,1,0],[0,2,0,1,0],[0,0,1,0,0],[0,0,-1,1,0]));

/*solución del apartado a*/
mb1b2:invert(mbcb1).mbcb2;

/*solucion del apartado b. p es el polinomio con coordenadas en la base canónica*/
p:[1,1,1,1,1];
penb1:invert(mbcb1).p;
penb2:invert(mbcb2).p;

Un cambio de base

mb2b3:matrix([1,0,0,0],[2,1,0,0],[-8,2,1,0],[20,-4,2,1]);
invert(mb2b3).[1,1,1,1];
invert(mb2b3);

## Aplicaciones lineales

### Ejercicio 5.07

diff((1+x)^3*x,x);

### Ejercicio 5.10

Vamos a comprobar que los vectores de los que me dan la imagen son una base

m:matrix([-1,1,3],[3,2,-1],[-2,1,1]);
/*Esta es la transpuesta de la matriz de cambio de base de la canónica a la base beta, formada por los vectores de los que me dan la imagen*/;

determinant(m);

Como el determinante de m es diferente de 0  tenemos que los vectores fila de beta son una base. Así que la aplicación g está bien definida

mbcb:transpose(m);

mbbcg:matrix([6,1,-2],[-4,14,-5],[16,-12,1]);

mbcbcg:mbbcg.invert(mbcb);

Comprobamos que el resultado es consistente calculando la imagen de los vectores de beta

mbcbcg.m[1]; /*m[1] saca la primera fila de la matriz m*/
mbcbcg.m[2];
mbcbcg.m[3];

Calculamos la dimensión de la imagen

dimIm:rank(mbcbcg);

Con lo cual sabemos que la dimensión del kernel es 1 porque ambas deben
sumar 3. Para calcular una base del kernel usamos nullspace

nullspace(mbcbcg);

### Ejercicio 5.11

Tomamos como base los vectores de los que le conocemos sus imágenes

b:matrix([1,1,1,1],[1,1,1,0],[1,0,1,0],[-1,-2,0,0]);

/*Efectivamente es una base*/
determinant(b);

mbbc3f:transpose(matrix([0,0,1],[0,0,-1],[1,1,-1],[1,1,1]));

c3

mbcb1:transpose(matrix([1,1,1,0],[1,1,0,0],[1,0,0,0],[1,1,1,1]));

mbcb2:transpose(matrix([1,1,1],[0,1,1],[0,0,1]));

mb2bc:invert(mbcb2);

mbbc4:invert(transpose(b));

mbc4bc3f: mbbc3f.mbbc4;

mb1b2:mb2bc.mbc4bc3f.mbcb1;

### Ejercicio 5.12

  Mf:matrix([1,0,0,0],[1,-1,1,2])$
  f(x,y,z,t):=Mf.[x,y,z,t];
  baseker:nullspace(Mf);
  baseim:columnspace(Mf);
  dimker:4-rank(Mf);
  dimim:4-dimker;
f(x,y,z,t)[1][1];
mbcbp:transpose(matrix([1,1],[0,2]));
mbcb:transpose(matrix([1,1,1,0],[1,1,0,0],[1,0,0,0],[1,1,1,1]));
mfbbp:invert(mbcbp).Mf.mbcb;

## Diagonalización

Ejercicio de diagonalización de los entregables

A:matrix([5,-3,-9,15],[3,-1,-9,15],[3,-3,-7,15],[2,0,-12,18]);
polcaracteristico:factor(determinant(A-x*ident(4)));
p:factor(charpoly(A,x));
solve(p=0,x);
multiplicities;
B2:A-2*ident(4);
B3:A-3*ident(4);
B8:A-8*ident(4);
dimV2:4-rank(B2);
nullspace(B2);
nullspace(B3);
nullspace(B8);

Queremos encontrar una matriz que sea diagonalizable y que tenga por
matriz diagonal a D

a:1$
b:3$
D:matrix([a,0,0,0],[0,a,0,0],[0,0,b,0],[0,0,0,b])$
P:matrix([1,0,0,1],[1,1,0,1],[1,1,1,0],[1,1,1,1])$
P:P.P$
A:P.D.invert(P)$
display(D,P,A);

## Polinomio de Taylor y funciones de una variable

### Polinomio de Taylor

T:taylor(g(x),x,0,10);

f:ratdisrep(T);

g(x):=cos(x);
a:0;
orden:10;
T:taylor(g(x),x,a,orden);
Ttrun:ratdisrep(T);
TaylorDef(xx):=ev(Ttrun,x=xx);
display(TaylorDef(x),g(x));

FuncionVsTaylor(f,a,b,orden):=block(
g(x):=f,
T:taylor(g(x),x,a,orden),
Ttrun:ratdisrep(T),
TaylorDef(xx):=ev(Ttrun,x=xx),
display(TaylorDef(x),g(x)),
plot2d([TaylorDef(x),g(x)],[x,a,b])
)$

for orden:1 thru 8 do FuncionVsTaylor(cos(x),0,%pi,orden);

plot2d([cos(x),sin(x)],[x,0,4]);

powerdisp;
powerdisp:true$
ratdisrep(taylor(sin(x)/x,x,0,10));
powerdisp:false$
ratdisrep(taylor(sin(x)/x,x,0,10));

### Polinomio de Taylor manual

kill(c)$
FormulaTaylor(f,x0,n):=block([w,PoliTay,R],
            PoliTay:ev(f,x=x0)+sum(diff(f,x,i)/i!*(z-x)^i,i,1,n),
            R:ev( diff(f,x,n+1),x=c ),
            ev(PoliTay,x=x0,z=x)+R*x^(n+1) /(n+1)!)$
PolinomioTaylor(f,x0,n):=block([w,PoliTay,R],
            PoliTay:ev(f,x=x0)+sum(diff(f,x,i)/i!*(z-x)^i,i,1,n),
            ev(PoliTay,x=x0,z=x))$
RestoDeLagrange(f,x0,n):=block([w,R],
            R:subst(c,x,diff(f,x,n+1)),R*x^(n+1)/(n+1)!)$

h1(x):=sin(3*x)$
h2(x):=cos(2*x)$
h3(x):=log(1+x)$
h(x):=h1(x)+h2(x)+h3(x)$
orden:2$
PolinomioTaylor(h(x),0,orden);
g(u):=subst(u,x,PolinomioTaylor(h(x),0,orden));

resto(u):=subst(u,x,RestoDeLagrange(h(x),0,orden));
resto1(u):=subst(u,x,RestoDeLagrange(h1(x),0,orden));
resto2(u):=subst(u,x,RestoDeLagrange(h2(x),0,orden));
resto3(u):=subst(u,x,RestoDeLagrange(h3(x),0,orden));

g(u);
resto(u);

### Discontinuidades

Discontinuidad de segunda especie

a:.5$
wxplot2d(sin(1/x),[x,-a,a]);

## Límites de funciones de varias variables

plot3d(2^(-x^2+y^2),[x,-3,3],[y,-2,2])$

### Ejercicio 11

plot3d(x*y*(x^2-y^2)/(x^2+y^2),[x,-2,2],[y,-2,2])$

plot3d((x^6+y^6)/(x^3-y^3),[x,-2,2],[y,-2,2])$

plot3d((x^6+y^6)/(x^3-y^3),[x,-10^-2,10^2],[y,10^-2,10^2])$

plot2d([x,x+x^5],[x,-2,2],[y,-2,2]);

g(x):=(x^6+(x+x^5)^6)/(x^3-(x+x^5)^3);

plot2d(g(x),[x,0,0.01]);

## Diferenciabilidad de funciones de varias variables

### Ejemplo

plot2d([1+2*(x-1),x^2],[x,-2,2]);

plot2d([[sin(x),[x,0,2*%pi]],[cos(x),[x,0,%pi]]]);

plot3d(y^2/x,[x,-10^-2,10^-2],[y,-10^-2,10^-2])$

a:1$plot3d(y^2/x,[x,-a,a],[y,-a,a])$

a:3$plot3d(sin(x+y),[x,-a,a],[y,-a,a])$

### Plano tangente y superficie

define(parcialx(x,y),diff(superficie(x,y),x))

Si queremos exportar la gráfica a un fichero (en este caso la he exportado a eps). El ordenador guardará la gráfica en el directorio que contenga la variable maxima_tempdir

maxima_tempdir:"/home/gabi/Sotad/CopyGSL/Docencia/Dibujos"$

plot3d([[superficie(x,y),[x,x0-a1,x0+a1],[y,y0-a1,y0+a1]],
    [ptangente(x,y),[x,x0-a2,x0+a2],[y,y0-a2,y0+a2]]],[gnuplot_term,ps],
  [gnuplot_out_file,"PlanoTangenteEsferaRadio2.eps"],[legend,"Plano tangente a una superficie"])$

define(parcialx(x,y),diff(superficie(x,y),x));

parcialx(0,0);

### Extremos relativos

plot3d(exp(x*y+50/x+20/y),[x,4.8,5.2],[y,1.8,2.2],[plot_format,gnuplot_pipes])$

f(x,y):=x*y*log(x^2+y^2);
a:1/sqrt(2*exp(1));
float(a);
amplitud:.5;

plot3d(f(x,y),[x,-1,1],[y,-1,1],[plot_format,gnuplot_pipes])$

#### Ejercicio 9.23.a y observaciones de composición con funciones estrictamente decrecientes

f(x,y):=x^2+(y-1)^2;
plot3d(f(x,y),[x,-1,1],[y,-0,2],[plot_format,gnuplot_pipes])$
g(x,y):=1/(f(x,y)+3)$
plot3d(g(x,y),[x,-1,1],[y,-0,2],[plot_format,gnuplot_pipes])$

s(x,y):=-f(x,y)$
plot3d(s(x,y),[x,-1,1],[y,-0,2],[plot_format,gnuplot_pipes])$

hessian(g(x,y),[x,y]);

### Extremos condicionados

plot3d (2^(-u^2 + v^2), [u, -3, 3], [v, -2, 2])$

#### Sobre el problema de números positivos que sumen una cantidad fijada

x:0.1$
y:26$
z:27-x-y$
V:x*y*z$
display(V);

V(x):=float(x*26*(1-x))$

V(0.1);
V(10^-6);
V(10^-12);

float(2+sqrt(8));

## Integración

### Sumas de Riemann

/* Two Maxima functions for introductory integral calculus
(c) 2016, themaximalist.org  */
/* RiemannSum
fn is the function to be integrated on the interval [a,b]
n rectangles
opt specifies:
0: left endpoints
1: midpoints
2: right endpoints
*/
RiemannSum(fn,a,b,n,opt):=
block([xx,s],
xx(i,n):= a+(i+opt/2)*(b-a)/n,
s: sum(ev(fn,x=xx(i,n)),i,0,n-1)*(b-a)/n,
float(s)
)$
/* RiemannRectangles to Draw rectangles
fn is the function to be integrated:
on the interval [a,b]
n rectangles
opt specifies:
0: left endpoints
1: midpoints
2: right endpoints
*/
load(draw)$
RiemannRectangles(fn,a,b,n,opt):=
block([xi,xx,rects,wxd],
xi(i,n):= a+i*(b-a)/n,
xx(i,n):= a+(i+opt/2)*(b-a)/n,
rects(n):=makelist(rectangle([xi(k,n),0], [xi(k+1,n),ev(fn,x=xx(k,n))]),k,0,n-1),
wxd(n):=apply(wxdraw2d,
append([xrange=[a-(b-a)/4, b+(b-a)/4],color=blue],
rects(n),
[transparent=true, explicit(fn,x,a, b)]
)
),
wxd(n)
)$

RiemannRectangles(sin(x),0,%pi,10,0);

### Ejercicios de los apuntes

## Bucles for y condicionales, ejemplos

/*NÚMEROS PRIMOS ENTRE M Y N*/
M:2000;
N:2050;
lista:[];
listanoprimos:[];
suma:0;
for j:M thru N do ( 
    if primep(j) then (lista:append(lista,[j]),suma:suma+j) 
                    else (listanoprimos:append(listanoprimos,[j]))
);
print("La lista de primos entre",M," y ",N," es: ",lista);
print("La lista de no primos entre",M," y ",N," es: ",listanoprimos);
print("La suma de la lista de primos es: ",suma);

/*SUMA DE NÚMEROS ENTRE LOS NÚMEROS NInicial Y NFinal*/
SumasEntreNumeros(NInicial,NFinal):=
block(
sumaimpares:0,
numeroimpares:0,
sumapares:0,
numeropares:0,
sumatotal:0,
NumeroFinal:801,
NumeroInicial:1, /*Esta variable blablabla*/
numerostotales:0,
for j:NumeroInicial thru NumeroFinal do(
        if remainder(j,2)=1 then (sumaimpares:sumaimpares+j,numeroimpares:numeroimpares+1)
                            else (sumapares:sumapares+j,numeropares:numeropares+1),
        sumatotal:sumatotal+j,
        numerostotales:numerostotales+1
),
print("La suma de los impares entre ",NumeroInicial," y ",NumeroFinal," es: ",sumaimpares),
print("La suma de los pares entre ",NumeroInicial," y ",NumeroFinal," es: ",sumapares),
print("La suma total de los números entre ",NumeroInicial," y ",NumeroFinal," es: ",sumatotal),
print("Entre ",NumeroInicial," y ",NumeroFinal," hay: ",numeropares, "números pares"),
print("Entre ",NumeroInicial," y ",NumeroFinal," hay: ",numeroimpares, "números impares"),
print("Entre ",NumeroInicial," y ",NumeroFinal," hay: ",numerostotales, "números")
);

SumasEntreNumeros(1,801);

/*LISTADO DE MÚLTIPLOS DE SEIS*/
j:6;
suma:0;¦
listamultiplosseis:[];
while (j<=43) do (  suma:suma+j,
                   listamultiplosseis:append(listamultiplosseis,[j]) ,
        j:j+6);
suma;
listamultiplosseis;
listamultiplosseis[3];
cardinality(setify(listamultiplosseis));
length(listamultiplosseis);

## Terraplanismo

### Calculando el hundimiento aparente

La función hundimiento da lo que se hunde un terreno al recorrer una distancia geodésica en la superficie terrestre. Ambas cantidades se miden en kilómetros

R:6371$
hundimiento(distancia):=float(R*(1-cos(distancia/R)))$
distancia(hundimiento):=float(R*acos(1-hundimiento/R))$

hundimiento(68);

hundimiento(25);

hundimiento(50);

hundimiento(100);

hundimiento(200);

distancia(0.1);

distancia(.5);

distancia(1);

distancia(2);

### Semejanza de triángulos

altura:float(68*410/31);

altura+hundimiento(68)*1000;

# Ejercicios entregables

## Ejercicios entregables del primer cuatrimestre

Ejercicio 16

f(x):=(x-1)*(x-3)*(x-12)$
g(x):=(x-1)*(x-3)$

A:abs(integrate(f(x)-g(x),x,1,3))+
abs(integrate(g(x)-f(x),x,3,13));

float(A);

h(x):=x^4/4-17/3*x^3+55/2*x^2-39*x;

h(3)-h(1)+(-h(13)+h(3));

Ejercicio 15

M:matrix([7*%i+2,3*%i+3],[7*%i,3*%i+7]);
N:matrix([7*%i+1,3*%i+1],[7*%i,7*%i+7]);
P:matrix([1,%i+2],[7,3]);

rectform(%i*M.N);

Ejercicio 16

expand(P.P);

i:%i;

Ejercicio 17

i*(M+N).P*i;

expand(%);

Ejercicio 18

expand(invert(P));

Ejercicio 19

invert(P).M;

rectform(%);

Ejercicio 20

A:matrix([2,7,3],[3,14,1],[7,4,3]);
B:matrix([7,3],[7,7],[3,1]);
C:matrix([3,8,4],[4,15,2],[8,5,4]);

A.B;

Ejercicio 21 

transpose(B).transpose(A);

Ejercicio 22

(A+ident(3))^^2;

Ejercicio 23

A^^2+2*A+ident(3);

Ejercicio 24

A.C;

Ejercicio 25

C.A;

Ejercicio 26

(A+C)^^2;

Ejercicio 27

A^^2+2*A.C+C^^2;

Ejercicio 28

A^^2+A.C+C.A+C^^2;

Ejercicio 29

H:matrix([2,7],[0,%i/2]);
invert(H);

Ejercicio 50

/*Datos del problema*/
beta4:matrix([1,1,1,1],[1,1,1,0],[2,1,0,0],[1,0,0,0])$
mbcb4:transpose(beta4)$
beta5:matrix([1,1,1,1,1],[1,1,2,1,0],[7,1,1,0,0],[0,1,0,0,0],[1,0,0,0,0])$
mbcb5:transpose(beta5)$
mbcbcf:matrix([1,1,2,1],[0,1,7,2],[0,0,1,3],[0,0,0,0],[0,0,0,0])$

mb4b5f:invert(mbcb5).mbcbcf.mbcb4;

Ejercicio 62

A:matrix([5,-1,-15,21],[3,1,-15,21],[3,-1,-13,21],[2,2,-18,24]);

p(x):=determinant(A-x*ident(4));

p(x);

expand(p(x));

factor(p(x));

 solve(p(x)=0,x);

rank(A-2*ident(4));

echelon(A-2*ident(4));

echelon(A-3*ident(4));

nullspace(A-2*ident(4));

nullspace(A-10*ident(4));

nullspace(A-3*ident(4));

A.[10,9,7,4];

A.[6,6,5,3];

P:transpose(matrix([0,3,4,3],[6,3,1,0],[1,1,1,1],[3,3,3,2]));

invert(P);

invert(P).A.P;

## Ejercicios de prácticas

Ejercicio 13

A:matrix([172,28,-7,0],[-972,-158,42,0],[-24,-4,11,0],[0,0,0,4]);
charpoly(A,x);
expand(%);
solve(%=0,x);
multiplicities;

Ejercicio 14

nullspace(A-10*ident(4));

nullspace(A-11*ident(4));

nullspace(A-4*ident(4));

Ejercicio 17

B:matrix([1,2,3],[2,4,6],[3,4,5]);
poli(x):=charpoly(B,x);
poli(x);

B:matrix([1,2,3],[2,4,6],[3,4,5]);
poli(x):=charpoly(B,x);
expand(poli(x));

(B-4*ident(3))^^2;

Ejercicio 21

f(x):=x^2-12*x+27;
g(x):=x^3-22*x^2+147*x-270;

solve(f(x)=g(x),x);

wxplot2d([f(x),g(x)],[x,3,11]);

area:abs(integrate(g(x)-f(x),x,3,9))+abs(integrate(f(x)-g(x),x,9,11));

Ejercicio 22

f(x):=x^4-35*x^3+433*x^2-2181*x+3510;
g(x):=x^3-22*x^2+147*x-270;
wxplot2d([f(x),g(x)],[x,3,14]);
solve(f(x)=g(x),x);
integrate(g(x)-f(x),x,3,9)+integrate(f(x)-g(x),x,9,10)+integrate(g(x)-f(x),x,10,14);

Ejercicio 25

f(x):=sqrt(6^2*(1-x^2/3^2));

volumen:%pi*integrate(f(x)^2,x,-3,3);

Ejercicio 27

ymas(x):=20+sqrt(6^2-x^2);
ymenos(x):=20-sqrt(6^2-x^2);

volumen:quad_qags(%pi*(ymas(x)^2-ymenos(x)^2),x,-6,6);

Ejercicio 28

ymas(x):=21+sqrt(6^2-x^2);
ymenos(x):=21-sqrt(6^2-x^2);

area:quad_qags(2*%pi*ymas(x)*sqrt(1+diff(ymas(x),x)^2)+2*%pi*ymenos(x)*sqrt(1+diff(ymenos(x),x)^2),x,-6,6);

## Ejercicios entregables del segundo cuatrimeste

Ejercicios 1 a 4

f(x,y):=[3*y^12+3*x^4+3,x^3*y^11+2,2*x^2*y^2];
g(x,y,z):=[3*x*y*z+y+x];
jacobianog:jacobian(g(x,y,z),[x,y,z]);
jacobianof:jacobian(f(x,y),[x,y]);
f(1,1);

aux1:ev(jacobianof,x=1,y=1);
aux2:ev(jacobianog,x=9,y=3,z=2);

En la solución siguiente tenemos el jacobiano de g compuesto con f en el punto (1,1) con el que puedes 
sacar la ecuación del plano tangente.

aux2.aux1;

f(1,1);

g(9,3,2);

Ejercicio 18

vi:integrate(integrate(integrate(14*log(z)*r+sin(r^2)*r,r,4,7),t,0,2*%pi),z,1,12);

float(vi);

/*Mucha atención, los argumentos de los cosenos y los senos van en radianes*/

# Curso del Profesor David López Medina con el que puedes 
aprender los comandos básicos estructurados por temas. 

## Álgebra lineal

### Vectores

Definición: v=(2,-3), w=(1,2)

v:[2,-3];
w:[1,2];

Componentes: v_2

v[2];

Suma: v+w

v+w;

Producto por escalares: -3·v

-3*v;

Producto escalar: v·w

v.w;

Norma: ||v||

sqrt(v.v);

### Matrices

Definición: A=(1 2), B=(0  1), C=(1 2 3)
              (3 4)    (-2 1)    (4 5 6)
                                 (7 8 9)

A:matrix([1,2],[3,4]);
B:matrix([1,1],[-1,1]);
C:matrix([1,2,3],[4,5,6],[7,8,9]);

Coeficientes: c_{2,3}

C[2][3];

Suma: A+B

A+B;

Producto por escalares: 3·A

3*A;

Producto: A·B (OJO: no es A*B)

A.B;
A*B;

Matriz por vector: A·v

A.v;

Traspuesta: A^t

transpose(A);

Inversa: A^-1 (OJO: no es A^-1)

invert(A);
A^-1;

Determinante: |A|

determinant(A);

Rango: rg(C)

rank(C);

Núcleo: Ker(C)

nullspace(C);

args(%);

Imagen: Im(C)

columnspace(C);
args(%);

### Sistemas lineales

SCD: x+y=2,2x+3y=2

linsolve([x+y=2,2*x+3*y=2],[x,y]);

SI: x+y=2,2x+2y=2

linsolve([x+y=2,2*x+2*y=2],[x,y]);

SCI (solución paramétrica): x+y=2,2x+2y=4

linsolve([x+y=2,2*x+2*y=4],[x,y]);

SCI (solución implícita): x+y=2,2x+2y=4

linsolve_params:false$
linsolve([x+y=2,2*x+2*y=4],[x,y]);

### Diagonalización

Polinomio característico: p_A(x)

charpoly(A,x);

expand(%);

Autovalores: lambda_i

eigenvalues(A);

solve(charpoly(A,x),x);

Autovectores: v_i

eigenvectors(A);

## Cálculo

### Cálculo diferencial

Límites finitos: lim(x->0) sen x/x

limit(sin(x)/x,x,0);

Límites infinitos: lim(x->-inf) (2x-1)^3/(x^2+1)

limit((2*x-1)^3/(x^2+1),x,-inf);

Series numéricas: sum 1/n^2

simpsum;
sum(1/n^2, n, 1, inf);
simpsum:true;
sum(1/n^2, n, 1, inf);
reset(simpsum)$

Derivadas: d^3 sin 2xy/dx^3

diff(sin(2*x*y),x,3);

Jacobiana: J(x+y^2,sin(x·y),x^y)

jacobian([x+y^2,sin(x*y),x^y],[x,y]);

Hessiana: H(sin(x+2y))

hessian(sin(x+2*y),[x,y]);

Taylor: sin x/x, n=10, x0=0

taylor(sin(x)/x,x,0,10);

ratdisrep(%);

powerdisp;
powerdisp:true$
ratdisrep(taylor(sin(x)/x,x,0,10));
powerdisp:false$
ratdisrep(taylor(sin(x)/x,x,0,10));

### Cálculo integral

Primitiva: int x/(x^3-2x^2+1)

integrate(x/(x^3-2*x^2+1),x);

Definida: int[0,pi] cos^2 x

integrate(cos(x)^2,x,0,%pi);

Impropia: int[-inf,inf] e^(-x^2/2)

integrate(exp(-x^2/2),x,-inf,inf);

rootscontract(%);

### Cálculo vectorial

Gradiente: grad x·y^2·sin(z)

load(vect)$
express(grad(x*y^2*sin(z)));

ev(%,diff);

Divergencia: div (x,z·y^2,sin(x+z))

express(div([x,z*y^2,sin(x+z)]));
ev(%,diff);

Rotacional: rot (x^2,y^3,z)

express(curl([x^2,y^3,z]));
ev(%,diff);

Laplaciano: lap x·y^2

express(laplacian(x*y^2));
ev(%,diff);

## Variable compleja

### Operaciones

Definición: z=1+i, omega=2·e^(pi·i/6) (OJO: Im z antes de Re z)

z:1+%i;
%omega:2*exp(%pi*%i/8);

powerdisp:true$
z;

Operaciones: z+omega, z*omega

z+%omega;
z*%omega;

Parte real: Re(omega)

realpart(%omega);

Parte imaginaria: Im(omega)

imagpart(%omega);

Forma binómica: Re(omega)+i·Im(omega)

rectform(%omega);

Módulo: |z|

abs(z);

Argumento: arg(z)

carg(z);

Forma exponencial: |z|·e^(i·arg(z))

polarform(z);

Conjugado: z^*,omega^*

conjugate(z);
conjugate(%omega);

### Series

Definición: f(z)=(z+1)/(z^3+4z^2-3z-18) (OJO: z definido)

f(z):=(z+1)/(z^3+4*z^2-3*z-18)$
f(z);

kill(z)$
f(z):=(z+1)/(z^3+4*z^2-3*z-18)$
f(z);

Numerador: p(z)

num(f(z));

Denominador: q(z)

denom(f(z));

Singularidades: q(z)=0

solve(denom(f(z)),z);

factor(denom(f(z)));

Serie de potencias: f(z), z0=0

powerseries(f(z),z,0);

niceindices(%);

Serie de Laurent: f(z), z0=-3

niceindicespref;
niceindicespref:[n,m]$
niceindices(powerseries(f(z),z,-3));

Residuo: Res(f,2)

residue(f(z),z,2);

## Ecuaciones diferenciales

### EDO de primer orden

Separadas: y'=sen x/sqrt(y)

ode2('diff(y,x)=sin(x)/sqrt(y),y,x);

Lineales: y'=y/x+1, y(1)=3

ode2('diff(y,x)=y/x+1,y,x);

ic1(%,x=1,y=3);

Homogéneas: y'=(x+y)/(2x+y)

ode2('diff(y,x)=(x+y)/(2*x+y),y,x);

Factor integrante: y'=(4x^3-y)/(2x-x^4/y)

ode2('diff(y,x)=(4*x^3-y)/(2*x-x^4/y),y,x);

### EDO de segundo orden

Coeficientes constantes: y''+2y'+2y=sen x

ode2('diff(y,x,2)+2*'diff(y,x)+2*y=sin(x),y,x);

Condiciones iniciales: y''+2y'+2y=sen x, y(pi/2)=3, y'(pi/2)=5

ic2(%,x=%pi/2,y=3,'diff(y,x)=5);

Condiciones frontera: y''+2y'+2y=sen x, y(0)=-2/5, y(pi)=2/5

bc2(ode2('diff(y,x,2)+2*'diff(y,x)+2*y=sin(x),y,x),x=0,y=-2/5,x=%pi,y=2/5);

Euler: x^2·y''+x·y'+y=1, y(1)=3, y'(1)=-2

ode2(x^2*'diff(y,x,2)+x*'diff(y,x)+y=1,y,x);
ic2(%,x=1,y=3,'diff(y,x)=-2);

### Transformada de Laplace

Transformada: L(log(x))

laplace(log(x),x,z);

Heaviside: L(H_2)

laplace(unit_step(x-2),x,z);

Transformada inversa: L^-1((z-1)/(z^3+z))

ilt((z-1)/(z^3+z),z,x);

Ecuaciones lineales: y1'=y1+y2, y2'=y1-y2, y1(0)=2, y2(0)=3

desolve(['diff(y1(x),x)=y1(x)+y2(x),'diff(y2(x),x)=y1(x)-y2(x)],[y1(x),y2(x)]);

atvalue(y1(x),x=0,2)$
atvalue(y2(x),x=0,3)$
desolve(['diff(y1(x),x)=y1(x)+y2(x),'diff(y2(x),x)=y1(x)-y2(x)],[y1(x),y2(x)]);

## Cálculo numérico

### Aritmética

Coma flotante: sqrt(2)

sqrt(2);
sqrt(2.0);

float(sqrt(2));

Calculadora: sqrt(2)+pi/3+e^2

sqrt(2)+%pi/3+%e^2;
numer:true$
sqrt(2)+%pi/3+%e^2;

Precisión baja: sqrt(2), N=6 (OJO: doble precisión interna)

fpprintprec:6$
sqrt(2);
sqrt(2)^2-2;

Precisión alta: 1/3, N=50 (OJO: numer lo anula)

reset(fpprintprec)$
fpprec:50$
bfloat(1/3);
numer:false$
bfloat(1/3);

Tiempo de CPU: sum(n=1,25000) 1/n^2

showtime:true$
sum(1/n^2,n,1,25000);
showtime:false$

### Factorizaciones

LU: (0  1 2)
    (1  1 3)
    (-1 0 5)

lu_factor(matrix([0,1,2],[1,1,3],[-1,0,5]));
get_lu_factors(%);

Cholesky: (4   2 -2)
          (2   3 -3)
          (-2 -3 14)

cholesky(matrix([4,2,-2],[2,2,-3],[-2,-3,14]));

### Ecuaciones no lineales

Bisección: cos(x)=x, [a,b]=[0,1]

find_root(cos(x)=x,x,0,1);

Newton: cos(x)=x, x0=0

load(mnewton)$
mnewton(x-cos(x),x,0);

Newton: x+cos(y)=2,y+sin(x)=3, (x0,y0)=(0,0) (OJO: y depende de x)

mnewton([x+cos(y)-2,y+sin(x)-3],[x,y],[0,0]);

diff(y,x);
kill(y)$
diff(y,x);
mnewton([x+cos(y)-2,y+sin(x)-3],[x,y],[0,0]);

### Interpolación

Lagrange: p(0)=1, p(2)=-1, p(3)=0

load(interpol)$
lagrange([[0,1],[2,-1],[3,0]]);
expand(%);

### Cuadratura

Cuadratura: int[0,1] exp(-x^2/2)/sqrt(2·pi)

quad_qags(exp(-x^2/2)/sqrt(2*%pi),x,0,1);

### Ecuaciones diferenciales

Runge-Kutta: y1''=-y1/(y1^2+y2^2)^(3/2), y2''=-y2/(y1^2+y2^2)^(3/2), y1(0)=1, y2(0)=0, y1'(0)=0, y2'(0)=1, T=2*pi

rk([y3,y4,-y1/(y1^2+y2^2)^(3/2),-y2/(y1^2+y2^2)^(3/2)],
  [y1,y2,y3,y4],
  [1,0,0,1],
  [t,0,2*%pi,0.1]
);

last(%);

## Gráficas

### Funciones en 2-D

Una: y=sen x en [0,2·pi]

plot2d([sin(x),cos(x)],[x,0,2*%pi])$

Varias: y=sen x, y=cos x, y=tg x en [-2·pi,2·pi] (OJO: asíntotas)

plot2d([sin(x),cos(x),tan(x)],[x,-2*%pi,2*%pi])$

plot2d([sin(x),cos(x),tan(x)],[x,-2*%pi,2*%pi],[y,-2,2])$

Poligonal: y=sen x en [0,2·pi]

l:makelist(2*%pi*n/10,n,0,10);
plot2d([discrete,l,sin(l)])$

l:makelist(2*%pi*n/100,n,0,100);
plot2d([discrete,l,sin(l)])$

Puntos: y=sen x en [0,2·pi]

plot2d([discrete,l,sin(l)],
  [style,points]
)$

Paramétricas: Círculo (OJO: escalas distintas, picos)

plot2d([parametric,cos(t),sin(t)],[t,0,2*%pi])$

plot2d([parametric,cos(t),sin(t)],[t,0,2*%pi],
  [gnuplot_preamble, "set size ratio 1;"]
)$

plot2d([parametric,cos(t),sin(t)],[t,0,2*%pi],
  [nticks,100],
  [gnuplot_preamble, "set size ratio 1;"]
)$

### Superficies en 3-D

Superficie: z=2^(-x^2+y^2) en [-3,3]x[-2,2]

plot3d(2^(-x^2+y^2),[x,-3,3],[y,-2,2])$

Perspectiva: z=2^(-x^2+y^2) en [-3,3]x[-2,2]

plot3d(2^(-x^2+y^2),[x,-3,3],[y,-2,2],
  [azimuth,0],
  [elevation,90]
)$

Malla: z=2^(-x^2+y^2) en [-3,3]x[-2,2]

plot3d(2^(-x^2+y^2),[x,-3,3],[y,-2,2],
  [palette,false]
)$

Cilíndicas: Cono

plot3d(r,[r,0,1],[theta,0,2*%pi],
  [transform_xy,polar_to_xy]
)$

Esféricas: Casquete esférico

plot3d(1,[theta,0,2*%pi/3],[phi,0,2*%pi],
  [transform_xy,spherical_to_xyz]
)$

### Campos de direcciones

Ecuación: y'=(y^2-1) (OJO: ejes excesivos)

load(plotdf)$
plotdf(y^2-1)$

plotdf(y^2-1,
  [xradius,3],
  [yradius,3]
)$

Sistema: y1'=y2, y2'=-y1

plotdf([y,-x])$

Trayectoria: y''-(1-y^2)·y'+y=0

plotdf([y,(1-x^2)*y-x],
  [xradius,5],
  [yradius,5],
  [trajectory_at,2,0],
  [direction,forward]
)$

### Animaciones

EDO: y''+y=sen omega·x, y(0)=y'(0)=0 (OJO: omega definido)

ic2(ode2('diff(y,x,2)+y=sin(%omega*x),y,x),x=0,y=0,'diff(y,x)=0);

kill(%omega)$
ic2(ode2('diff(y,x,2)+y=sin(%omega*x),y,x),x=0,y=0,'diff(y,x)=0);

l:makelist(1+i/20,i,20,1,-1);
with_slider(
  %omega,
  l,
  (%omega*sin(x))/(%omega^2-1)-sin(%omega*x)/(%omega^2-1),
  [x,-10*%pi,10*%pi],
  [y,-15,15]
)$

Direcciones: y1'=y2, y2'=-a·y1

plotdf(
  [y,-2^b*x],
  [sliders,"b=-2:2"]
)$

### Gráficos avanzados

Implícitas: x^2+4y^2=4

load(draw)$
draw2d(implicit(x^2+4*y^2=4,x,-2.5,2.5,y,-1.5,1.5))$

h

debugmode(true);

Paramétricas: banda Moëbius (OJO: plano visible)

draw3d(
  xu_grid=50,
  yv_grid=25,
  parametric_surface(cos(a)*(3+b*cos(a/2)),sin(a)*(3+b*cos(a/2)),b*sin(a/2),
    a,-%pi,%pi,b,-1,1)
)$

draw3d(
  xu_grid=50,
  yv_grid=25,
  parametric_surface(cos(a)*(3+b*cos(a/2)),sin(a)*(3+b*cos(a/2)),b*sin(a/2),
    a,-%pi,%pi,b,-1,1),
  surface_hide=true
)$

Regiones: Circunferencia inscrita

draw2d(
  x_voxel=200,
  y_voxel=200,
  region(abs(x)+abs(y)<1 and x^2+y^2>0.49,x,-1.5,1.5,y,-1.5,1.5)
)$

Rellenos: f(x)=2-3·x^2, g(x)=sqrt(1-x^2/4)

draw2d(fill_color=grey,
  filled_func=2-3*x^2,
  explicit(sqrt(1-x^2/4),x,-0.59,0.59),
  filled_func=false,
  line_width=2,
  color=red,
  explicit(2-3*x^2,x,-sqrt(2/3),sqrt(2/3)),
  color=blue,
  explicit(sqrt(1-x^2/4),x,-1,1)
)$

Vectores: f(x)=x^2, x=1

draw2d(explicit(x^2,x,-2,2),
  line_width=2,
  color=red,
  explicit(1+2*(x-1),x,-2,2),
  color=black,
  point_size=3,
  point_type=filled_circle,
  points([[1,1]]),
  line_width=1,
  head_length=0.1,
  vector([1.5,0.5],[-0.4,0.4]),
  label(["punto de tangencia",1.5,0.3]),
  xtics={["x",1]},
  ytics_rotate=true,
  ytics={["f(x)",1]},
  grid=true,
  title="recta tangente"
)$

Contornos: e^(-x^2-y^3/100)

draw3d(color=gray,
  explicit(exp(-x^2-y^3/100),x,-1,1,y,-4,4),
  yv_grid=10,
  contour_levels=10,
  contour=both,
  surface_hide=true
)$

Múltiples: f(x)=1,x,x^2,x^3

multiplot_mode(screen)$
draw2d(
  title="Constante",
  allocation=[[0,1/2],[1/2,1/2]],
  xtics=none,
  ytics=none,
  explicit(1,x,-1,1)
)$
draw2d(
  title="Recta",
  allocation=[[1/2,1/2],[1/2,1/2]],
  xtics=none,
  ytics=none,
  explicit(x,x,-1,1)
)$
draw2d(
  title="Parabola",
  allocation=[[0,0],[1/2,1/2]],
  xtics=none,
  ytics=none,
  explicit(x^2,x,-1,1)
)$
draw2d(
  title="Cubica",
  allocation=[[1/2,0],[1/2,1/2]],
  xtics=none,
  ytics=none,
  explicit(x^3,x,-1,1)
)$
multiplot_mode(none)$

### Dibujo de curvas 3D en paramétricas

Elipse corte del plano x+y+z=1 y x^2+y^2=1

draw3d(
parametric(cos(u),sin(u),1-cos(u)-sin(u),u,0,2*%pi),
color= turquoise,
line_width = 2)$

Intersección de plano y cilindro

draw3d(
explicit(1-x-y,x,-1,1,y,-1,1),
color = royalblue,
surface_hide = true,
parametric_surface(cos(u),sin(u),v,
u, 0, 2*%pi, v, -2.5, 2.5),
color= pink,
surface_hide = true,
line_width = 10,
parametric(u,0,0,u,0,1),
color= red,
line_width = 10,
parametric(0,u,0,u,0,1),
color= red,
line_width = 10,
parametric(0,0,u,u,-2.5,2.5),
color= red,
line_width = 10,
parametric(cos(u),sin(u),1-cos(u)-sin(u),u,0,2*%pi),
color= red,
line_width = 10,
parametric(u,0,0,u,0,1),
color= red,
line_width = 10,
title = "Surface & curves" )$

Ejemplo del manual. Varios objetos juntos

draw3d(
explicit(exp(sin(x)+cos(x^2)),x,-3,3,y,-3,3),
color = royalblue,
parametric(cos(5*u)^2,sin(7*u),u-2,u,0,2),
color= turquoise,
line_width = 2,
parametric(t^2,sin(t),2+t,t,0,2),
surface_hide = false,
title = "Surface & curves" )$

Caracola de mar

draw3d(title
= "Sea shell",
xu_grid
= 100,
yv_grid
= 25,
view
= [100,20],
surface_hide
= true,
parametric_surface(0.5*u*cos(u)*(cos(v)+1),
0.5*u*sin(u)*(cos(v)+1),
u*sin(v) - ((u+3)/8*%pi)^2 - 20,
u, 0, 13*%pi, v, -%pi, %pi) )$

### Malla de puntos en 3D

load(draw)$
draw3d(
mesh([[1,1,3],
[7,3,1],[12,-2,4],[15,0,5]],
[[2,7,8],
[4,3,1],[10,5,8], [12,7,1]],
[[-2,11,10],[6,9,5],[6,15,1], [20,15,2]])) $

## Programación

### Bucles

For: 10! (OJO: bucle no sale por pantalla)

x:1;
for n:1 thru 10 do x:x*n;

x:1;
for n:1 thru 10 do (
  x:x*n,
  disp(x)
)$

While: 10!

x:1$
n:0$
while (n<10) do (
  n:n+1,
  x:x*n,
  disp(x)
)$

### Condicionales

If: Fibonnacci

f(n):=if (n=1 or n=2) then (1) else (f(n-2)+f(n-1))$
for n:1 thru 10 do disp(f(n))$

### Salida por pantalla

Evaluada: x

x;

disp(x)$

No evaluada: x

'x;
disp('x)$

Nominal: x

display(x)$

Compuesta: x

print("La variable ",'x," vale 10!=",x)$

### Funciones

Sencillas: f(x,y)=(x,sin(y),H(sin(x+y)) (OJO: (x,y) definido)

kill(x,y)$
f(x,y):=[x,sin(y),hessian(sin(x+2*y),[x,y])]$
f(x,y);

Procedimientos: sum(i=1...n) i, n! (OJO: variable globales)

f(n):=block(
  a:0,
  x:1,
  for i:1 thru n do print('sum(j,j,1,i),"=",a:a+i,", ",i,"!=",x:x*i)
)$

f(5)$

a;
x;
kill(a,x)$
a;
x;

f(n):=block([a,x],
  a:0,
  x:1,
  for i:1 thru n do print('sum(j,j,1,i),"=",a:a+i,", ",i,"!=",x:x*i)
)$

f(5)$

x;
a;

### Ejemplos prácticos

LU(A,[no_output]): Factorización LU de una matriz cuadrada A (con output paso a paso opcional)

LU(A,[no_output]):=block([n,L,U,l,u,e], /*Variables bloqueadas*/
local(l,u,e), /*Vectores locales*/
globalsolve:true, /*Parámetro de la función linsolve para asignar la solución*/
n:matrix_size(A)[1], /*Dimensión de la matriz A*/
for i:1 thru n do 
  (
  for j:1 thru n do (
    if i>j then u[i,j]:0 else if i=j then l[i,j]:1 else l[i,j]:0 /*Coeficientes de L y U conocidos*/
    )
  ),
L:genmatrix(l,n,n),U:genmatrix(u,n,n), /*L y U antes de empezar*/
if emptyp(no_output)=true then print("Sistema original: ",L,".",U,"=",A), /*Salida por pantalla (condicional)*/ 
for i:1 thru n do 
  (
  for j:1 thru n do
    (
    e[i,j]:sum(l[i,k]*u[k,j],k,1,min(i,j))=A[i][j], /*Cada ecuación*/
    linsolve(e[i,j],if i>j then l[i,j] else u[i,j]), /*Solución de cada ecuación*/
    if emptyp(no_output)=true then ( /*Salida por pantalla (condicional)*/ 
      print(""), /*Una línea en blanco*/
      print("Ecuación fila ",i," columna ",j,":",e[i,j],
", Solución:",if i>j then 'l[i,j]=l[i,j] else 'u[i,j]=u[i,j]), /*Ecuaciones por orden*/
      print("Sistema actualizado: ",L:genmatrix(l,n,n),".",U:genmatrix(u,n,n),"=",A) /* Sistema actualizado*/
     )
 )
  ),
globalsolve:false, /*Parámetro de la función linsolve para no asignar la solución*/
L:genmatrix(l,n,n),U:genmatrix(u,n,n), /*L y U definitivas*/
[L,U] /*Salida de la función*/
)$

LU(matrix([4,2,-2], [2,2,-3],[-2,-3,14]))$

construye_regla(c,[no_output]): Coeficientes y error de una regla de cuadratura (con output paso a paso opcional)

construye_regla(c,[no_output]):=block([A,X,E,B,b,C,p,n], /*Variables bloqueadas*/
local(b), /*Vectores locales*/
globalsolve:true, /*Parámetro de la función linsolve para asignar la solución*/
powerdisp:true, /*Parámetro para ordenar los subíndices de menor a mayor*/
n:length(c), /*Dimensión del sistema*/
A:transpose(vandermonde_matrix(c)), /*Matriz del sistema*/
X:makelist(b[i],i,1,n), /*Vector de coeficientes */
for i:1 thru n do e[i]:A[i].X=1/i, /*Cada ecuación del sistema*/
E:makelist(e[i],i,1,n), /*Vector de ecuaciones*/
B:linsolve(E,X), /*Solución del sistema*/
C:0,
for i:n while(is(equal(C,0))) do (
  p:i, /*Orden*/
  C:expand((1/(p+1)-sum(b[i]*c[i]^p,i,1,n))/p!) /*Constante de error*/
  ),
if emptyp(no_output)=true then ( /*Salida por pantalla (condicional)*/ 
  print("El sistema es ",transpose(E)),
  powerdisp:false, /*Parámetro para escribir (b-a) y no (-a+b)*/
  print("Los coeficientes son ",B," y ",makelist(a[i]=('b-'a)*b[i],i,n)),
  print("La regla es ",integrate('f('x),'x,'a,'b),"~",('b-'a)*sum(b[i]*'f('a+c[i]*('b-'a)),i,1,n)),
  print("El orden es p=",p," y la constante es ",'C[p],"=",C),
  print("El error es ",integrate('f('x),'x,'a,'b)-('b-'a)*sum(b[i]*'f('a+c[i]*('b-'a)),i,1,n),"~",C*('b-'a)^(p+1),at('diff('f('x),'x,p),'x='theta))
),
[makelist(b[i],i,1,n),p,C]
)$

construye_regla([0,1/2,1])$

Forma_Lagrange(x,f,t,[no_output]): Función que calcula en la forma de Lagrange el polinomio interpolador de la función f(t) basado en los nodos x (con output paso a paso opcional)

Lagrange(x,f,t,[no_output]):=block([l,L,n], /*Variables bloqueadas*/
n:length(x)-1, /*Grado del polinomio*/
for i:0 thru n do l[i]:product((t-x[j+1])/(x[i+1]-x[j+1]),j,0,i-1)*product((t-x[j+1])/(x[i+1]-x[j+1]),j,i+1,n), /*Base de Lagrange*/
L:sum(subst(t=x[i+1],f)*l[i],i,0,n), /*Polinomio interpolador*/
if emptyp(no_output)=true then ( /*Salida por pantalla (condicional)*/ 
  plot2d([f,L,[discrete,makelist([x[i],subst(t=x[i],f)],i,1,n+1)]],[t,lmin(x),lmax(x)],
    [style,lines,lines,points],[point_type, bullet],[legend,f,sconcat("L(",t,")"),""],
    [xlabel,""],[color,red,blue,magenta]),
  print("Los nodos son ",makelist('x[i]=x[i+1],i,0,n))
),
L
)$

n:4$
Lagrange(makelist(-5+10*i/n,i,0,n),1/(1+x^2),x)$

n:8$
Lagrange(makelist(-5+10*i/n,i,0,n),1/(1+x^2),x)$

n:16$
Lagrange(makelist(-5+10*i/n,i,0,n),1/(1+x^2),x)$

Chevyshev(n,I,f,t,[no_output]): Función que calcula los nodos de Chevyshev para grado n en el intervalo I y dibuja el polinomio interpolador de f(t) basado en ellos (con output paso a paso opcional)

Chevyshev(n,I,f,t,[no_output]):=block([a,b,x,L], /*Variables bloqueadas*/
[a,b]:I, /*Intervalo*/
x:float(makelist(a+(b-a)/2*(1+cos(%pi*(2*n-2*i+1)/(2*n+2))),i,0,n)), /*Nodos de Chevyshev (en coma flotante)*/
L:Lagrange(x,f,t,vacio),
if emptyp(no_output)=true then ( /*Salida por pantalla (condicional)*/ 
plot2d([f,L,[discrete,makelist([x[i+1],subst(t=x[i+1],f)],i,0,n)]],[t,lmin(x),lmax(x)],
    [style,lines,lines,points],[point_type, bullet],[legend,f,sconcat("L(",t,")"),""],
    [xlabel,""],[color,red,blue,magenta]), 
  print("Los nodos son ",makelist('x[i]=x[i+1],i,0,n))
),
L
)$

Chevyshev(4,[-5,5],1/(1+x^2),x)$

Chevyshev(8,[-5,5],1/(1+x^2),x)$

Chevyshev(16,[-5,5],1/(1+x^2),x)$

## LaTeX

### Salida

Expresiones: e^(sqrt(pi)) (OJO: entorno matemático)

tex(exp(%pi^(1/2)))$

tex1(exp(%pi^(1/2)));

Formales: int[0,pi]sen xdx

f(x):=sin(x)$
'integrate(f(x),x,0,%pi);
tex1(%);

Igualdades: int[0,pi]sen x·dx

'integrate(f(x),x,0,%pi)=integrate(f(x),x,0,%pi);
tex1(%);

Definición: epsilon_0

tex1(e0);
texput(e0,"\\varepsilon_{0}")$
tex1(e0);

Modificación: conjugado(z)

declare(z,complex)$
tex1(conjugate(z));

texconjugado(e):=block([a:first(e)],
  concat("\\overline\{",tex1(a),"\}")
)$
texput(conjugate,texconjugado)$
tex1(conjugate(z));

### Gráficas

Terminal (plot): sen x en eps

plot2d(sin(x),[x,0,2*%pi],
  [gnuplot_term,ps],
  [gnuplot_out_file,"plot.eps"]
)$

Terminal (draw): sen x en eps

draw2d(
  terminal='eps,
  file_name="draw",
  explicit(sin(x),x,0,2*%pi)
)$

Tipografía (plot): sen x en eps-LaTeX (OJO: incompleto, color, puntos)

plot2d(sin(x),[x,0,2*%pi],
  [gnuplot_term,"epslatex"],
  [gnuplot_out_file,"plotlatex.tex"]
)$

plot2d(sin(x),[x,0,2*%pi],
  [gnuplot_term,"epslatex standalone color solid"],
  [gnuplot_out_file,"plotlatex.tex"]
)$

Tipografía (draw): sen x en eps-LaTeX (OJO: 5.29.0, gplbacktext)

load(draw)$
draw2d(
  dimensions=[1600,900],
  terminal='epslatex_standalone,
  file_name="drawlatex",
  explicit(sin(x),x,0,2*%pi)
)$

### Interacción

Ejecución: latex-maxima.tex (OJO: requiere max-tex-tools.lisp, max-tex-init.mac y maxiplot.sty)

Reptición: latex-personalizado.tex

## Ejercicios

R:2^2$
f(x,y):=sqrt(1-x^2-y^2)$
normavector(v):=sqrt(v[1]^2+v[2]^2)$
v:[1,1]$
vv:[-1,1]$
A:[.5,.5]$
fv(punto,w,t):=f((punto+t*w/normavector(w))[1],(punto+t*w/normavector(w))[2])$
plot2d([fv(A,v,t),fv(A,vv,t)],[t,0,1])$

fv(0.5,0.5,t);

# Borrador

A:matrix([1,2,3,4],[2,3,4,1],[3,4,1,2],[4,1,2,3]);

determinant(A);

A:matrix([1,0],[0,0]);

B:matrix([0,1],[0,0]);

(A-B).(A+B);

A^^2-B^^2

A^^2-B^^2+A.B-B.A;

A^^2-B^^2-A.B+B.A;

A:matrix([1,-1,1],[-2,1,0],[1,1,2]);

determinant(A);

invert(A);

B:matrix([2,3,4],[3,1,2],[1,0,0]);

A.B
;

transpose(A.B);

transpose(A).transpose(B);

transpose(B).transpose(A);