# Interpolación polinomial con MAXIMA

## Interpolación lineal

La regla de interpolación más sencilla es la lineal. Todos sabemos (o deberíamos saber) que por dos puntos en el plano pasa una línea recta 
(según el Teorema de los Puntos Gordos, por tres puntos también para una recta, si estos son suficientemente gordos)

La ecuación de una recta que pasa por los puntos (x1,y1) y (x2,y2) es:

    f(x) = y2 + (y2-y1)/(x2-x1) * (x-x2) 

/* Así definimos la recta que une los puntos (x1,y1) con (x2,y2) */
recta(x,x1,y1,x2,y2) := y2 + (x-x2)*(y2-y1)/(x2-x1);

/* Por ejemplo la recta que pasa por los puntos (1,2) y (3,1) es */
recta(x,1,2,3,1);

/* Representamos esta recta */
wxplot2d(recta(x,1,2,3,1),[x,0,4]);

/* Por si alguien no se fía, veamos que pasa por los nodos, dibujando los puntos correspondientes con la opción discrete de wxplot2d */
nodos : [[1,2],[3,1]];
wxplot2d(
    [recta(x,1,2,3,1),[discrete,nodos]],
    [x,0,4],
    [style, [lines,2,1],[points,3,2,1]]
    );

 
Si en lugar de 2, tenemos varios puntos, entonces podemos realizar una interpolación lineal entre ellos, es decir, podemos construir la línea poligonal que los une a  todos. Esta función es continua a trozos y unirá dos puntos consecutivos mediante el correspondiente segmento. Para ello utilizamos el paquete "interpol"  (¡anda, como el cuerpo de policía internacional!) de MAXIMA que hay que cargar previamente con el comando 

load(interpol);

/* Aquí se carga en memoria la librería de funciones interpol */
load(interpol);

 
Para realizar Interpolación lineal, usaremos los siguientes comandos:

linearinterpol(lista de puntos)
charfun (prop)
charfun2(x,a,b)

/* Vamos a definir la función anterior usando la librería interpol */

/* Vamos a borrar la función recta que hemos creado antes (para que se vea que no se hacen trampas) */
remfunction(recta);
functions;

/* Usamos linearinterpol, en este caso, en la definición, los puntos no están dentro de la función y los puntos vienen dados como vector de vectores */

define( recta(x),linearinterpol([[1,2],[3,1]]));

/* Y la pintamos otra vez (con sus dos bolicas rojas) */
wxplot2d(
    [recta(x),[discrete,nodos]],
    [x,0,4],
    [style, [lines,2,1],[points,3,2,1]]
    );

/* Como hemos dicho, podemos usar más puntos */ 
define(g(x),linearinterpol([[2,2],[3,3],[4,1],[5,4]]));

/* y dibujamos para comprobar que es una línea poligonal (eso para las visitas, en confianza es una línea quebrada) */
wxplot2d(g(x),[x,2,5]);

/* Se me olvidaron las bolicas rojas */
nodos : [ [2,2], [3,3], [4,1], [5,4] ];
wxplot2d(
    [g(x),[discrete,nodos]],
    [x,2,5],
    [style, [lines,2,1],[points,3,2,1]]
    );

Se ha hecho uso de la función 

    charfun(prop)

que es una función definida en la librería interpol, que devuelve el valor 0 si el argumento prop es FALSO y devuelve 1 si el argumento prop es VERDADERO, 
por ejemplo:

charfun(3<1);

charfun(3>1);

Otra función incluida en la librería interpol es

    charfun2(x,a,b)

o función característica de intervalo. Devuelve el valor 1 si a<= x < b, y 0 en caso contrario.

charfun2(2,1,3);

charfun2(4,1,3);

charfun2(3,1,3);

charfun2(2,1,3);

Si x no está determinado, el resultado viene dado en términos de charfun

remvalue(x);
charfun2(x,1,10);

## Interpolación de Lagrange

 
El caso n=2 (2 puntos), se puede (y se debe) extender a un número arbitrario de puntos. La cosa es: dados (n+1) pares de puntos (x0,y0),....(xn,yn), buscamos el polinomio de grado <=n que pasa por todos los nodos.

No es una poligonal como antes, sino un polinomio que pasa por todos los puntos.

/* Vamos a empezar desde el principio */

/* Nos cargamos todas las cosas, para que no aparezcan cosas raras (como las matemáticas) */
kill(all)$;

/* Empecemos despacio, caso 2 nodos (sí, sé que se ha visto, pero de otro modo) */

/* Aquí los nodos (la variable independiente, la "x" de toda la vida) */
nodos:[1,2];

/* Aquí los valores de la función buscada en los nodos (la variable dependiente, vamos la variable "y", la compañera de la "x" en estos negocios) */
valor:[3,7];

/* Resumen: buscamos una función de grado 1 (recordemos que para n+1 nodos el grado es n, por tanto si hay 2 nodos el grado es .....tachán....¡1! */

/* 
    Buscamos una función de la forma 
    
    f(x) = a*x +b 
    
    que pase por los puntos (1,3) y (2,7), o sea, que cumple
    
    f(1) = 3 ==> a*1 + b = 3
    f(2) = 7 ==> a*2 + b = 7
    
    Usemos MAXIMA:
*/  

/* Definimos la función */
define(f(x),a*x+b);

/* Planteamos las 2 ecuaciones */
ec1: f(nodos[1])=valor[1];
ec2: f(nodos[2])=valor[2];

/* Resolvemos el sistema en los coeficientes a y b (que son los que no conocemos, claro)*/
solve([ec1,ec2],[a,b]);

/* y ya tenemos la función */;
define(f(x),4*x-1);

/* Pintamos (por si las moscas) */
nodos : [ [1,3], [2,7] ];
wxplot2d(
    [f(x),[discrete,nodos]],
    [x,1,2],
    [style, [lines,2,1],[points,3,2,1]]
    );

Vemos un ejemplo con 3 nodos y por tanto un polinomio de grado 2

/* Aquí las "x" y las "y" */
nodos:  [1,2,3];
valores: [3,7,1];

/* La función que buscamos: polinomio de grado 2 (ya sabes 3 nodos, grado 2, 1001 nodos, grado 1000 )*/
define(p(x),a*x^2+b*x+c);

/* Las ecuaciones (que obviamente será 3, tantos como nodos) */
ec1: p(1) = 3;
ec2: p(2) = 7;
ec3: p(3) = 1;

/* Resolvemos el sistema lineal en los coeficientes*/
solve([ec1,ec2,ec3],[a,b,c]);

/* Et voilà!, ya tenemos la función */
define(p(x),-5*x^2+19*x-11);

/* que pintamos, por si quedaran dudas, con los puntos correspondientes  */
puntos : [[1,3], [2,7], [3,1] ]$
wxplot2d(
    [p(x),[discrete,puntos] ],
    [x,1,3],
    [style, [lines,2,1],[points,3,2,1]]
    );

 
Obviamente si n es grande, este procedimiento parece demasiado largo (¡toma 20 nodos y veras qué risa!), 
así que vamos a hacerlo usando el paquete interpol

/* Borramos todo */
kill(all)$

/* Cargamos el paquete */
load(interpol);

/* Definimos los puntos por donde pasará el polinomio, notar que ahora todos están juntitos y por parejas */
puntos : [[1,3],[2,7],[3,1]];

/* y usamos la función lagrange que lleva interpol */
define(p(x),lagrange(puntos));

/*  Vamos a darle pinta de polinomio */
define(p(x),expand(p(x)));

/* Dibujando, que es gerundio */
puntos : [[1,3], [2,7], [3,1] ]$
wxplot2d(
    [p(x),[discrete,puntos] ],
    [x,1,3],
    [style, [lines,2,1],[points,3,2,1]]
    );

Un ejemplo con muchos, muchos nodos (bueno, solo 5)

/* Lo hacemos todo de golpe, que ya tenemos experiencia en estas cosas */
puntos : [ [1,3], [2,7], [3,1], [4,-2], [5,3] ]$
define(p(x),expand(lagrange(puntos)));
wxplot2d(
    [p(x),[discrete,puntos] ],
    [x,1,5],
    [style, [lines,2,1],[points,3,2,1]]
    );

También es posible usar los puntos dentro de una matriz de dos columnas

puntos : matrix([1,3], [2,7], [3,1], [4,-2], [5,3])$
define(p(x),expand(lagrange(puntos)));

## Polinomio de Taylor

La función

    taylor (f(x),x,x0,n)

Devuelve el desarrollo de Taylor de grado "n" de la función "f(x)", centrada en "x0"

f(x0) + f'(x0)*(x-x0)/1! + f''(x0)*(x-x0)^2/2! + f'''(x0)*(x-x0)^3/3! +......

/* Empezamos por lo fácil, la función coseno */
f(x):=cos(x);

/* Su polinomio de taylor de orden 2, centrado en 0, será: */
taylor(f(x),x,0,2);

/* Vemos los puntos suspensivos que nos indican que se puede seguir (hasta el infinito y más allá)*/
/* No obstante estos puntos no molestan y podemos definir funciones */
define(t2(x),taylor(f(x),x,0,2));

/* Que podemos pintar, para compararlas con la función original */
wxplot2d([f(x),t2(x)],[x,0,2]);

/* Se puede aumentar el orden */
taylor(cos(x),x,0,5);

/* O cambiar de función y de punto */
taylor(log(x),x,1,7);

/* También podemos expresar el polinomio sin puntos suspensivos y en orden decreciente de potencias */
expand(taylor(log(x),x,1,7));

Aumentar el grado conduce a una mejor aproximación cerca del punto donde se centra la serie

/* Probemos con uno de orden 8*/
f(x):=cos(x);
define( 
    p8(x),  
    expand(taylor(f(x),x,0,8))
);
wxplot2d([f(x),p8(x)],[x,-4,4],[y,-2,2]);

/* Parece que todo ha ido perfecto, pero si aumentamos la escala, veremos que hay cierta discrepancia */
wxplot2d([f(x),p8(x)],[x,-8,8],[y,-2,2]);

/* Eso es debido a que la aproximación del polinomio a la función es de tipo local, todos sabemos que el coseno está acotado entre -1 y 1 y un polinomio ....pues no.
El problema se "resuelve" aumentando el grado el polinomio */

define( 
    p14(x),  
    expand(taylor(f(x),x,0,14))
);
wxplot2d([f(x),p14(x)],[x,-8,8],[y,-2,2]);

/* Claro que esto no les pasa a todas las funciones. Vemos lo que le ocurre a la función arco tangente */

/* Empezamos por orden 2 */
define(p2(x),taylor(atan(x),x,0,2));
wxplot2d([atan(x),p2(x)],[x,-8,8],[y,-2,2]);

/* Seguimos por orden 4 */
define(p4(x),taylor(atan(x),x,0,4));
wxplot2d([atan(x),p4(x)],[x,-8,8],[y,-2,2]);

/* Aumentamos a grado 8 */
define(p8(x),taylor(atan(x),x,0,8));
wxplot2d([atan(x),p8(x)],[x,-8,8],[y,-2,2]);

/* Y nos arriesgamos con grado 64 */
define(p64(x),taylor(atan(x),x,0,64));
wxplot2d([atan(x),p64(x)],[x,-8,8],[y,-2,2]);

/* Parece que la cosa funciona bien sólo en el intervalo [-1,1] */

Ya hemos visto que si nos molestan los puntos suspensivos, podemos eliminarlos con expand, aunque tenemos otras opciones

    ratdisrep(expresion): devuelve el polinomio de Taylor convertido en expresión racional, sin puntos suspensivos, es más o menos, como expand

    taytorat(expresion): devuelve el polinomio de taylor sin puntos suspensivos y con denominador común

Para manejar los polinomios de Taylor como funciones en MAXIMA, es preferible usar una de estas opciones.

p: taylor(cos(x),x,0,8);
q: ratdisrep(taylor(cos(x),x,0,8));
r : taytorat(taylor(cos(x),x,0,8));

En los polinomios de Taylor es preferible usar las potencias en orden creciente, esto lo podemos conseguir cambiándole el valor a la variable powerdisp, que por defecto tiene el valor false, como podemos comprobar

powerdisp;

Si cambiamos su valor a true

powerdisp : true;

Ordenamos según potencias crecientes

ratdisrep(taylor(cos(x),x,0,8));

Vemos con un ejemplo, la utilidad de "tantas" funciones para el polinomio de Taylor

/* Potencias crecientes */
powerdisp: true$ 

/* Definimos un par de funciones */
f(x):= sin(x); 
g(x):= tan(x);

/* Calculamos el polinomio de Taylor de orden 5 de cada función, centrado en 0 */
taylor(f(x),x,0,5);  
taylor(g(x),x,0,5);  

/* Calculamos la diferencia entre estos desarrollos, que vuelve a ser otro polinomio de Taylor */
taylor(g(x),x,0,5)-taylor(f(x),x,0,5);

/* Si hacemos la diferencia con la función */
f(x)- taylor(f(x),x,0,5);

/* Pero si hacemos la diferencia usando expand o trunc, obtendremos la verdadera "diferencia" entre el polinomio y la función */
(f(x)-trunc(taylor(f(x),x,0,5)));
(f(x)-expand(taylor(f(x),x,0,5)));

Otras funciones relacionadas con el polinomio de Taylor son: 

La función

    taylorp(expresion)

indica si "expresion" es un polinomio de Taylor o un polinomio "normal", la diferencia está en que un polinomio de Taylor incluye los puntos suspensivos al final.

taylorp(p);
taylorp(q);

La función:

    coeff(expr,x,n): Devuelve el coeficiente de la potencia x^n

coeff(x^2+x+2,x,2);

coeff(x^2+x+2,x^2);

coeff(taylor(cos(x),x,0,8),x^4);

Comprobamos en esta animación, cómo evoluciona el ajuste del polinomio con la función en función del  a la función con el tiempo

f(x):= sin(x)+cos(x);
with_slider(
n,
makelist(i,i,1,20),
[f(x),taylor(f(x),x,0,n)],[x,-10,10],[y,-3,3]
);

### Resto de Lagrange

Podemos definir el resto de Lagrange usando su definición:
                   D_(n+1) f(c)
R_n  =      -------------- (x-x0)^(n+1)
                         (n+1)!

donde D_n indica la derivada (n+1)-ésima de la función f, evaluada en un punto c, que está entre x y x0.

R(n,x0):=diff(subst(c,x,f(x)),c,n+1)*(x-x0)^(n+1) /(n+1)!$
f(x):= log(1+x); taylor(f(x),x,0,4); 
R(4,0); 
f(x) = expand( trunc(taylor(f(x),x,0,4)) ) + R(4,0);

## Polinomio de Taylor en varias variables

Usamos la misma función, pero indicando todas las variables

F(x,y):=exp(x+y);
define(G(x,y),expand(taylor(F(x,y),[x,y],[0,0],[3,3])));

expand(%);

wxplot3d( 
    [ F(x,y),G(x,y),
    [x,-1,1],[y,-1,1]]);

/* Plano tangente = Polinomio de Taylor de orden 1 en los puntos (1,2) y (0,0)*/
F(x,y):= 14-x^2-y^2;
define(G0(x,y),expand(taylor(F(x,y),[x,y],[1,2],[1,1])));
define(G1(x,y),expand(taylor(F(x,y),[x,y],[0,0],[1,1])));

/* Polinomio de Taylor de orden 2 en (0,0) */
F(x,y):= (x^2-3*x)*exp(y^2);
define(G(x,y),expand(taylor(F(x,y),[x,y],[0,0],2)));

/* Polinomio de Taylor de orden 2 en (1,1,0) */
F(x,y,z):= x/y - z*exp(x);
define(G(x,y,z),expand(taylor(F(x,y,z),[x,y,z],[1,1,0],2)));

/* Polinomio de Taylor de orden 3 en (1,1) */
F(x,y):= x^y;
define(G(x,y),expand(taylor(F(x,y),[x,y,z],[1,1],3)));