# Ecuaciones diferenciales ordinarias (EDO)

## EDO de primer orden usando "ode2"

La estrategia consiste en hallar la solución general (ode2)
y determinar el parámetro mediante la condición inicial (ic1),
más, si fuera necesario, despejar la x.

print("Ponemos nombre a la ecuación")$
E1 : 'diff(x,t)=t*(x+1);

print("Calculamos la solución general y le ponemos nombre")$
SGE1 : ode2(E1,x,t),fullratsimp;

print("Hacemos intervenir la condición inicial para fijar el parámetro")$
ic1( SGE1,t=1,x=3);

kill(all)$
print("Ponemos nombre a la ecuación")$
E1 : 'diff(x,t) = (exp(t)+1)/x;

print("Calculamos la solución general y le ponemos nombre")$
SGE1 : ode2(E1,x,t),fullratsimp;

print("Hacemos intervenir la condición inicial para fijar el parámetro")$
ic1( SGE1,t=1,x=3);

print("Despejamos x")$
solve(ic1( SGE1,t=1,x=3), x);

kill(all)$
print("Ponemos nombre a la ecuación")$
E1 : t^2 * 'diff(x,t) + 3*t*x = sin(t)/t;

print("Calculamos la solución general y le ponemos nombre")$
SGE1 : ode2(E1,x,t),fullratsimp;

print("Hacemos intervenir la condición inicial para fijar el parámetro")$
ic1( SGE1,t=%pi,x=0);

## EDO de segundo orden usando "ode2"

La estrategia consiste en hallar la solución general (ode2)
y determinar los 2 parámetros mediante las condiciones iniciales (ic2),
o de valores frontera (bc2)
más, si fuera necesario, despejar la x.

kill(all)$
print("Ponemos nombre a la ecuación")$
E1 : 'diff(x,t,2) -10* 'diff(x,t) + 9*x = 5*t;

print("Calculamos la solución general y le ponemos nombre")$
SGE1 : ode2(E1,x,t),fullratsimp;

print("Hacemos intervenir las condiciones iniciales 
para fijar los dos parámetros")$
ic2( SGE1,t=0,x=-1, 'diff(x,t)=0);

fullratsimp(%);

kill(all)$
print("Ponemos nombre a la ecuación")$
E1 : 'diff(x,t,2) -10* 'diff(x,t) + 9*x = 5*t;

print("Calculamos la solución general y le ponemos nombre")$
SGE1 : ode2(E1,x,t),fullratsimp;

print("Hacemos intervenir los valores frontera 
para fijar los dos parámetros")$
bc2( SGE1,t=0,x=-1, t=1,x=2);

## EDO usando "contrib_ode"

load(contrib_ode)$

E2 : t*(t+1)*'diff(x,t,2)+(t+5)*'diff(x,t,1)+(-4)*x;

print("Resolución con ode2:",  ode2(E2,x,t))$

contrib_ode(E2,x,t);
bc2( %[1], t=1,x=2, t=2,x=3);

E1 : 'diff(y,x,2)+(a*x+b)*y;
S1 : [y = bessel_y(1/3,2*(a*x+b)^(3/2)/(3*a))*%k2*sqrt(a*x+b)
         +bessel_j(1/3,2*(a*x+b)^(3/2)/(3*a))*%k1*sqrt(a*x+b)];

print("Veamos si S1 es solución de la ecuación E1, lo cual corresponde 
a que la respuesta sea 0. En este caso la respuesta es: ",
ode_check(E1,S1[1]))$

## La transformada de Laplace

La transforma de Laplace de la función f(t) es la función F(s)
definida mediante la fórmula siguiente,
supuesto que la integral impropia exista. 

A continuación hay algunos ejemplos:

laplace(1,t,s);

laplace(t^2,t,s);
    
assume(n+1>0)$ laplace(t^n,t,s);
    
laplace(cos(a*t),t,s);

print("La transformada de Laplace es un operador lineal")$
'laplace(a*f(t)+b*g(t),t,s) = laplace(a*f(t)+b*g(t),t,s);

print("Conportamiento respecto a la derivada")$
laplace('diff(f(t),t),t,s);
laplace('diff(f(t),t,2),t,s);

Para muchas funciones es posible recuperar la función original
a partir de su transformada de Laplace; 
el comandi "ilt" realiza la inversa de la transformada de Laplace

f(t):= t^2;
print("Su transforma de Laplace, que denotamos con F, es ", F : laplace(f(t),t,s))$
print("Aplicando   ilt se recupera la función f a partir de la F")$
ilt(F,s,t);

## EDO de primer y segundo orden usando "laplace"

La transformada de Laplace, en algunos casos, permite transformar un problema de EDO 
con valores iniciales, en un problema algebráico de solución más sencilla.
Resuelto el cual, mediante la inversa de la transformada de Laplace
se consigue obtener la solución del problema original

 Resuelta anteriormente con "ode2"
 ahora la resolvemos con "laplace"

Lo primero a tener en cuenta es que la incógnita x,
a diferencia de lo que ocurre con  "ode2" 
hay que escribirla como   x(t)  no simplemente como x

kill(all)$
print("Ponemos nombre a la ecuación")$
E1 : 'diff(x(t),t,2) -10*'diff(x(t),t)+ 9* x(t) = 5*t;

print("Hallamos la transformada de Laplace y le ponemos nombre")$
E2 : laplace(E1,t,s);


print("Aplicamos a E2 las condiciones iniciales y le ponemos nombre")$
E3 : ev(E2,x(0)=0,at('diff(x(t),t),t=0)=1);

print("Resolvemos la ecuación algebráica E3 tomando como incógnita  
laplace(x(t),t,s) ")$
E4 : solve(E3, laplace(x(t),t,s));

print("Convertimos la lista en ecuación ordinaria E5")$
E5 : E4[1];

print("Finalmente aplicamos a E5 la transformada de Laplace inversa;
(nótese el cambio de orden entre  t y s al utilizar ilt)
y así obtenemos la solucion del problema original")$
ilt(E5,s,t);

fullratsimp(%);

## EDO de segundo orden y sistemas lineales usando "desolve"

El comando  "desolve" permite resolver ecuaciones diferenciales
ordinarias de orden 2 con condiciones iniciales en punto 0

 Funciona para condiciones iniciales en el punto 0
 

kill(all)$
EC1: 'diff(x(t),t,2)+'diff(x(t),t)  +2* x(t) = sin(t);

print("Establecemos condiciones iniciales en t=0 
antes de utilizar desolve")$
atvalue(x(t), t=0, 0)$ atvalue( 'diff(x(t),t), t=0, 0)$

print("Obtenemos la solución con desolve")$
desolve( EC1, x(t));

OTRO PROCEDIMIENTO USANDO   "ode2"
La ecuación diferencial ha de ser escrita de un modo diferente

kill(all)$
EC1: 'diff(x,t,2)+'diff(x,t)  +2* x = sin(t);

SGE1 : ode2( EC1,x,t);

ic2( SGE1,t=0, x=0, 'diff(x,t)=0);

kill(all)$
print("Escribimos las ecuaciones y les ponemos nombre")$
EC1:'diff(f(x),x)=sin(x)+'diff(g(x),x);  
EC2:'diff(f(x),x)+x^2-f(x)=2*'diff(g(x),x,2);

print("Establecemos condiciones iniciales en t=0 
antes de utilizar desolve")$
atvalue(f(x),x=0,1);  atvalue('diff(g(x),x),x=0,a);

SOL : desolve([EC1,EC2],[f(x),g(x)])$
print("Separamos las dos soluciones f(x) y g(x)")$
SOL[1]; 
SOL[2];

kill(all)$
print("Escribimos las ecuaciones y les ponemos nombre")$
EC1:'diff(x(t),t)=2*x(t)-3*y(t);  
EC2:'diff(y(t),t)=-4*x(t)+5*y(t);

print("Establecemos condiciones iniciales en t=0 
antes de utilizar desolve")$
atvalue(x(t),t=0,1);  atvalue(y(t),t=0,1);

SOL : desolve([EC1,EC2],[x(t),y(t)])$
print("Separamos las dos soluciones x(t), y(t)")$
SOL[1]; 
SOL[2];

load(contrib_ode);
/* comentar su existencia y referir al manual de Maxima*/;

EC : 'diff(x,t)=t-x^2;

contrib_ode(EC,x,t);

## Métodos numéricos

Si los métodos analíticos fracasan, métodos numéricos, como el de Runge-Kutta, permiten
valores aproximados y hacen posible el dibujo de la solución.
Requiere cargar el paquete  dynamics

### Paquete "dynamics" 

load(dynamics)$

EC : 'diff(x,t)=t-x^2;

ode2(EC,x,t);

Maxima no sabe resolver la ecuación diferencial.
Vamos a tratar, al menos, de visualizar la curva x(t)
solución de la ecuación

puntos:rk(t-x^2, x,1, [t,0,20,0.1])$ 
wxplot2d([discrete,puntos]);

a:7$  /* cambiar a 10, 20, 50, 100,... */
vdp1 : rk([y,-x-(x^2-1)*y],[x,y],[1,1],[t,0,a,0.1])$
vdpp1: makelist([vdp1[k][2],vdp1[k][3]],k,1,10*a)$
wxplot2d([discrete,vdpp1])$

### Paquete "plotdf"

El paquete  "plotdf" permite crear gráficos del campo de direcciones de una 
ODE unidimensional y de campos de vectores en el plano, dando la posibilidad 
adicional de elegir órbitas de forma interactiva

load(plotdf)$

Una ecuación diferencial, tal como  y'=x-y^2, corresponde al problema de encontrar una función
y(x) que verifique la ecuación diferencial anterior. La información que proporciona
la ecuación es que para cada (x,y) del plano puedo conocer la dirección 
de la curva solución en dicho punto. El comando "plotdf" dibuja ese campo de direcciones. 

print("Campo de direcciones para y'=x-y^2, 
pinche en un punto para obtener la correspondiente
órbita que pasa po él.")$
plotdf(x-y^2,[x,y],[x,-1,5],[y,-5,5])$

print("Ecuación diferencial y'= exp(-x) +y 
con órbita pasando por un punto concreto;
pinche en otros puntos para obtener otras órbitas")$
plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$

La ecuación de Duffing   m*x"+c*x'+k*x+b*x^3 = 0   es de segundo orden, 
donde los coeficientes son parámetros.

Para representar el campo de direcciones de la misma podemos hacer el cambio
de variable   x'=y  que la transforme en un sistema lineal de dos ecuaciones
diferenciales 
    x'=y
    y'=-(c*y+k*x+b*x^3)/m

plotdf([y,-(k*x + c*y + b*x^3)/m],
              [parameters,"k=-1,m=1.0,c=0,b=1"],
              [sliders,"k=-2:2,m=-1:1"],[tstep,0.1])$

plotdf([e*(y-(x^3/3-x)),a-x],[x,-5,5],[y,-5,5],
   [parameters,"a=0,e=1"],[sliders,"e=1:20,a:0:2"],[nsteps,500],[xfun,"x^3/3-x"])$

plotdf([e*(y-(x^3/3-x)),a-x],[x,-3,3],[y,-3,3],
   [parameters,"a=0.983,e=8"],[sliders,"a=0.982:0.984"],[nsteps,500],[xfun,"x^3/3-x"])$