# Resolver ecuaciones

Planteado con total generalidad, la resolución de ecuaciones
es un tema difícil: sólo para determinadas tipologías es posible
encontrar una solución exacta, para la mayor parte únicamente 
se consiguen obtener una solución aproximada.
 

## Soluciones exactas con "solve"

 Ejemplos de soluciones usando "solve"

print("1. Para este polinomio Maxima encuentra las cuatro soluciones
reales y complejas")$
solve(-4*x^3 - 12*x^2 + 12, x);
   

print("2. Pero para éste no es capaz de hacerlo")$ 
solve(x^5-2*x^2-1);
   
print("3. Nuevos ejemplos con funciones más complicadas")$
solve(2*log(x)=-1,x);
solve (asin (cos (3*x))*(f(x) - 1), x);
print("Maxima advierte de que pueden faltar soluciones")$
        /* Un reto para reflexionar */

print("4. Los sistemas lineales son más fáciles, también para Maxima")$
    /* Observe que no es necesario poner  =0 en las dos ecuaciones */
solve([2*x-7*y-2,x-y+3], [x,y]); 

print("5. Los polinomios pueden tener raíces múltiples")$
solve(x^3-x^2-5*x-3,x); multiplicities;
x^3-x^2-5*x-3= factor(x^3-x^2-5*x-3);

Vamos a resolver a continuación un problema de interpolación de Lagrange
consistente en determinar el polinomio p(x) de cuarto grado que pasa por
los 5 puntos siguientes [7,2],[8,2],[1,5],[3,2],[6,7].

A tal fin definimos un polinomio de cuarto grado con coeficientes a determinar
en base a los datos disponibles

print("Definimos un polinomio genérico de 4º grado")$
kill(all)$
p(x) := a*x^4+b*x^3+c*x^2+d*x+e;

print("Creamos una ecuación para cada punto y le ponemos nombre")$
e1 : ev(p(x),[x=7])=2; 
e2 : ev(p(x),[x=8])=2; 
e3 : ev(p(x),[x=1])=5; 
e4 : ev(p(x),[x=3])=2; 
e5 : ev(p(x),[x=6])=7;

print("Ahora resolvemos el sistema de ecuaciones")$
Sol : solve([e1,e2,e3,e4,e5],[a,b,c,d,e]);

### La representación gráfica es útil: reciclar salidas 

Para visualizar gráficamente la solución necesitamos sustituir esos valores
de a, b,... en la fórmula de p(x). Eso se puede hacer copiando y pegando a mano.

Pero sería lindo que ese trabajo mecánico lo hiciese Maxima;
puesto que ya sabe que a=??, b= ?? 

a;

Pues no lo sabe. 
¡¡ Claro, recordemos que para asignar valores Maxima utiliza : y no = !!

Hay que extraer, de esa lista anidada las partes de la derecha    rhs 
de las diferentes igualadades de la lista. Lo hacemos por etapas y
practicamos la gestión del comando   rhs 

Sol;
    /* la lista principal sólo tiene un elemento que recuperamos */
SOL: Sol[1];  

    /* ahora asignamos valoares a los parámetros desconocidos
       usando las partes de la derecha en los = de los elementos de SOL */
[a,b,c,d,e] : [ rhs(SOL[1]), rhs(SOL[2]), rhs(SOL[3]), rhs(SOL[4]), rhs(SOL[5]) ];
  
    /* ahora substituimos a,b,... en p(x) por los valores concretos
       que acabamos de asignar usando los  : 
subst([a=a, b=b, c=c,d=d, e=e], p(x));
P: %;

kill(a,b,c,d,e)$

Hemos sabido comunicarle lo que queríamos, con un lenguaje que Maxima ha entendido.

Veamos las gráficas y hacemos también prácticas con draw2d.
Si ha configurado su sistema para que Maxima cargue draw cuando lo necesite,
no es necesario que ejecute la celda siguiente; 
si no ha hecho la configuración ha de ejecutarla

load(draw)$

wxdraw2d( 
explicit(P, x, 0,9),
    color=red, point_type=circle,
points([[7,2],[8,2],[1,5],[3,2],[6,7]])
)$

## Soluciones aproximadas: realroots, allroots, algsys

f(x):= x- x^2 + 1+log (1+x)$
solve(f(x),x);

Lo anterior, realmente, no es una solución: a la derecha del signo =
sigue apareciendo la x.

Vamos a dibujar (eso es fácil) para hacernos una idea de la situación

wxdraw2d( 
explicit( f(x), x, -0.9, 5),
    color=red, 
explicit(0, x, -0.9, 5) )$

La gráfica nos da la pista: 
    hay una solución entre -0.9 y 1
    otra solución entre 1 y 3

A partir de ahí podemos encontrar soluciones aproximadas

s1 : find_root( f(x)=0, x,-0.9,1);

s2 : find_root( f(x), x, 1,3);

print("Comprobamos el valor de f en esos puntos")$
f(s1);
f(s2);

Es fácil argumentar que esas son las únicas soluciones (aproximadas) de f(x)=0.

f(x) := x^5-2*x^2-1;
solve( f(x),x);

wxdraw2d( 
explicit(f(x), x, -3,3),
    color=red,
explicit (0, x, -3,3)   
     )$

En el intervalo dibujado claramente hay, al menos, una cero del polinomio.
Podemos dibujar en un intervalo más pequeño, por ver si la situación se clarifica

wxdraw2d( 
explicit(f(x), x, -2,2),
    color=red,
explicit (0, x, -3,3));

find_root( f(x),x, -1, 2);
f(%);

Al ser f(x) un polinomio podríamos haber prescindido de las consideraciones
anteriores para elegir el intervalo (y visualizar la situación)
y haber usado dos comandos que sólo son aplicables a polinomios

realroots(f(x));
%,numer;

allroots(f(x));

Esta técnica corrobora que ese es el único cero real (aproximado) de f.

print("Otro comando para buscar ceros de polinomios")$
algsys([f(x)],[x]);

### Sistemas polinómicos de ecuaciones: "algsys"

Abordamos aquí la solución de un sistema de ecuaciones polinómicas
con varias variables, cuya solución, como se muestra, puede depender de parámetros 

poli1 : 2*x*(1-a)$
poli2 : b-a$
poli3 : a*(-y - x^2 + 1)$
poli4 : b*(y - (x - 1)^2)$

algsys( [poli1,poli2,poli3,poli4],[x,y,a,b]);

Otro ejemplo que determina los puntos de corte entre una circunferencia y una elipse

wxdraw2d( proportional_axes=xy, xrange=[-2.5,2.5], yrange=[-1.5,1.5],
          implicit (x^2+y^2=1,x,-1,1,y,-1,1),
          color=red,
          implicit(x^2/4+4*y^2=1,x,-2,2,y,-1,1))$      
algsys([x^2+y^2=1,(x-1)^2+y^2=1/4],[x,y]);

### Comparando métodos

En un polinomio cuyos ceros elegimos a priori vamos a comparar los tres métodos

kill(f)$
f(x):=(x-0.1)*(x+3)*(x^2+x+1);

print("Con solve", solve(f(x) ))$

kill(f)$
f(x):=(x-0.1)*(x+3)*(x^2+x+1);

print("Con solve", solve(f(x)) )$

print("Con realroots", realroots(f(x)))$
%,numer;

print("Con allroots", allroots(f(x)))$
%,numer;
            
print("Con algsys", algsys([f(x)],[x]))$

## En ocasiones, colaborar con Maxima permite obtener 
soluciones exactas

EC : log(a)-log(a+1)=log(x-1)-log(x);
solve(EC, x);

Maxima fracasa y no encuentra la solución. 
Le decimos que "contraiga" los logaritmos

EC2 : logcontract(EC);
print("Escrita así es capaz de resolverla")$
solve(EC2, x);

La reescritura de la ecuación que Maxima ha realizado con nuestra ayuda
sugiere que otra estrategia de ayuda podría ser decirle que reescriba
la ecuación hallando su imagen mediante la función   exp  (exponencial)

EC3 : map(exp, EC);
solve(EC3,x);

EC : 1/2=exp(k*500);
solve(EC,k);

De nuevo le podemos ayudar diciendo que transforme la ecuación tomando logaritmos

map(log, EC);
solve(%,k);