# RESOLVIENDO EDOS CON MAXIMA

/*  RESOLVIENDO EDOS CON MAXIMA*/

## Resolviendo EDOs analíticamente

/* Resolviendo EDOs analiticamente */

edo2:'diff(z,x)=x-z;

ode2(edo2,z,x)$
expand(%);

expand(ic1(%,x=1,z=2));

define(s1(x),second(%));

wxplot2d(s1(x),[x,1,5])$

edo:'diff(y,x)+x*y=2*x;

ode2(edo,y,x)$
expand(%);

expand(ic1(%,x=0,y=1));

define(s1(x),second(%));

 wxplot2d(s1(x),[x,-5,5])$

edo3:'diff(w,x)=x/w;
   ode2(edo3,w,x)$
   sol:expand(%);
   expand(ic1(sol,x=0,w=3));

'diff(w,x)=1+w^2;
        ode2('diff(w,x)=1+w^2,w,x)$
        sol:expand(%);
        expand(ic1(sol,x=0,w=0));

sol:solve(%,w);

part(%,1);

 define(solw(x),second(%));

wxplot2d(solw(x),[x,0,%pi/2-.1])$

## Resolviendo EDOs numéricamente 

/* Resolviendo EDOs numericamente */

ode2('diff(z,x)=x-sin(z),z,x);

load(diffeq);

f(x,y):=1+y; h:1/20;

solnum:runge1(f,0,1,h,1);

 wxplot2d([discrete,solnum[1],solnum[2]])$

ode2('diff(w,x)=1+w,w,x)$
       sol:expand(%);
       expand(ic1(sol,x=0,w=1));

define(solw(x),second(%));

 wxplot2d([[discrete,solnum[1],solnum[2]],solw(x)],[x,0,1])$

load(dynamics)$

h:1/20;kill(x,y);
       numsolrk:rk(x-sin(y),y,1,[x,0,1,h])$

numsolrk;

wxplot2d([discrete,numsolrk],
         [xlabel, "x"], [ylabel, "y"],[color,blue])$

/* No siempre las cosas funcionan */
kill(all);

load(diffeq);

kill(f)$ /* Probar con 1.2, 1.4 y 1.6 */
x0:0$ x1:1.1$ y0:0$ h:1/20$
f(x,y):= 1+y^2;

solnum:runge1(f,x0,x1,h,y0)$

solnum[1];

wxplot2d([discrete,solnum[1],solnum[2]])$

/* Qué está ocurriendo */
'diff(w,x)=1+w^2;
        ode2('diff(w,x)=1+w^2,w,x)$
        sol:expand(%);
        expand(ic1(sol,x=0,w=0));

sol:solve(%,w);part(%,1);define(solw(x),second(%));

wxplot2d( [[discrete,solnum[1],solnum[2]],solw(x) ],[x,x0,x1]);

## Aplicaciones

/* APLICACIONES */
kill(all);

/* Ejemplo 1 geometría */;

kill(y)$
edo2:'diff(y,x,1) = 2*y/x;

ode2(edo2,y,x)$
expand(%);

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

define(soledo2(x),second(%));

wxplot2d(soledo2(x),[x,0,10], [xlabel, "x"], [ylabel, "y(x)"]);

/* Ejemplo 2 geometría */;

kill(w,x)$
edo3:'diff(w,x,1) = x+w;
ode2(edo3,w,x);

ic1(%,x=0,w=0);

define(s1(x),second(%));

wxplot2d(s1(x),[x,0,2], [xlabel, "x"], [ylabel, "y(x)"]);

/* Ejemplo 1 Circuito RL */;

kill(i,t,L,R,U)$
lru:L*'diff(i,t)+R*i-U=0;

ode2(lru,i,t);

ic1(%,t=0,i=0);

expand(%);

define(i(L,R,U,t),second(%));

plot2d( i(1,1,3,t) ,[t,0,10], [xlabel, "t"], [ylabel, "i(t)"]);

kill(i,t,L,R,U)$
lru:L*'diff(i,t)+R*i-U*sin(w*t)=0;

ode2(lru,i,t);

ic1(%,t=0,i=0);expand(%);

define(i(L,R,U,w,t),second(%));

float(i(1,1,3,2*%pi*50,t));

wxplot2d( float(i(1,1,3,2*%pi,t))  ,[t,0,10], [xlabel, "t"], [ylabel, "i(t)"]);

/* Ejemplo 4 ecuación barométrica */;

kill(p,l,h)$
pres:'diff(p,h)+l*p=0;

ode2(pres,p,h);
ic1(%,h=0,p=1);

define(pre(h,l),second(%));

float(lambda:1/8000);

float(pre(1654,lambda));

float(pre(3478,lambda));

float(pre(3710,lambda));

/* Ejemplo 5 Conducitividad del calor */;

kill(T,r,Q,a)$
Qq:a*'diff(T,r)*r^2+Q=0;

sol:ode2(Qq,T,r);

eq1:ic1(sol,r=r1,T=T1);

define(T(r,Q),second(%));

solve(eq1,Q);

define(Q(r,T),second(first(%)));

Q(r2,T2);

T(r,Q(r2,T2));

define(Te(r,r1,r2,T1,T2),ratsimp(%));

wxplot2d(Te(r,20,100,100,30),[r,34,100], [xlabel, "r"], [ylabel, "T(r)"]);

/* Ejemplo 6 reacciones químicas */

kill(x,t,k,a,b,sol)$
rq:'diff(x,t) = k*(a-x)*(b-x);

sol:ode2(rq,x,t);

ic1(%,t=0,x=0);

logcontract(%);

solve([%], [x]);

define(xx(t,k,a,b),factor(second(first(%))));

xx(t,1,10,7);
wxplot2d(xx(t,1,10,7),[t,0,2], [xlabel, "t"], [ylabel, "x(t)"]);

/* Ejemplo 8 velocidad de escape */

kill(vel,v,r,g,R,sol)$
vel:v*'diff(v,r) = - g*R^2/r^2 ;

sol:ode2(vel,v,r);

ic1(%,r=R,v=v0);

solve([%],v^2);

define(vq(r,g,R,v0),expand(second(first(%))));

limit(vq(r,g,R,v0),r,inf);

kill(all);

/* Ejemplo 9 Caida de un cuerpo en un medio viscoso */
kill(vel,v,t,r,k,sol)$
vel:'diff(v,t) = g- k*v^r ;

sol:ode2(vel,v,t);

ic1(%,t=0,v=v0);

assume(vf*w>1,g>0,k>0,w>0);

(1/g)*integrate(1/(1-w^2*v^2),v,v0,vf);

logcontract(%);

dede:%;

num(dede)-2*g*w*t;

solve(%,vf);

define(solcaso2(v0,w,g,t),second(first(%)));

limit(solcaso2(v0,w,g,t) ,t,inf);

/* Ejemplo 9 Caida de un cuerpo en un medio viscoso */
kill(vel,v,t,r,k,sol)$
vel:'diff(v,t) = g- k*v^2 ;

sol:ode2(vel,v,t);

ic1(%,t=0,v=v0);

solve(%,v);

assume(vf*w>1,g>0,k>0,w>0);

(1/g)*integrate(1/(1-w^2*v^2),v,v0,vf);

logcontract(%);

dede:%;

num(dede)-2*g*w*t;

solve(%,vf);

define(solcaso2(v0,w,g,t),second(first(%)));

limit(solcaso2(v0,w,g,t) ,t,inf);

kill(vel,v,t,r,k,sol)$
assume(vf*w>1,g>0,k>0,w>0);
vel:'diff(v,t) = g - k*v^3 ;

(1/g)*integrate(1/(1-w^2*v^3),v,v0,vf);

logcontract(%);

kill(vel,v,t,r,k,sol)$
assume (g>0, k>0,k>0,w>0)$
vel:'diff(v,t) = g*(1- w^2*v^3) ;

sol:ode2(vel,v,t);
ic1(%,t=0,v=v0);

logcontract(%);

kill(vel,v,t,r,k,sol)$
assume (g>0, k>0,k>0,w>0)$
vel:'diff(v,t) = g*(1- w^2*v^1.7) ;

sol:ode2(vel,v,t);
ic1(%,t=0,v=v0);

/* Resolucion numerica con Runge-Kutta */

/* Runge-Kuta 1. sintaxis runge(f, x0, x1, h, y0) donde f es la funcion y'=f(x,y)
x0 y x1 inicio y final del intervalo de integracion, h el paso, yo, ondición inicial */

load(diffeq);

kill(k,g,w)$
k:1$ g:9.8$ w:sqrt(k/g)$
f(x,y):= g*(1- w^2*y^1.7); h:1/20;

solnum:runge1(f,0,3,h,1)$

wxplot2d([discrete,solnum[1],solnum[2]])$

kill(k,g,w)$
k:1$ g:9.8$ w:sqrt(k/g)$ x0:0$ x1:3$
f(x,y):= g*(1- w^2*y^2); h:1/20;

solnum:runge1(f,0,3,h,1)$

solcaso2(1,w,g,t);

wxplot2d([[discrete,solnum[1],solnum[2]],solcaso2(1,w,g,t),1/w],[t,0,3]);

Nu:(x1-x0)/h;

makelist(x0+(x1-x0)*(k-1)/Nu,k,1,Nu+1);

solcaso2(1,w,g,%);

%-solnum[2];

/* Otra implementación del método de Runge-Kutta */;

/* Runge-Kuta 1. sintaxis rk(f,y,y0,[x,x0, x1, h]) donde f es la funcion y'=f(x,y)
x0 y x1 inicio y final del intervalo de integracion, h el paso, yo, ondición inicial */

load(dynamics)$

kill(k,g,w)$
k:1$ g:9.8$ w:sqrt(k/g)$ x0:0$ x1:3$
f(t,v):= g*(1- w^2*v^2); h:1/20;

numsolrk:rk(f(t,v),v,1,[t,0,3,h])$

numsolrk;

wxplot2d( [ [discrete,numsolrk], 1/w],[x,x0,x1] ,[xlabel, "t"], [ylabel, "v"],
[legend,"solución","1/w"], [color,blue,red]);

kill(k,g,w)$
k:1$ g:9.8$ w:sqrt(k/g)$ x0:0$ x1:3$
f(t,v):= g*(1- w^2*v^1.7); h:1/20;

numsolrk2:rk(f(t,v),v,1,[t,0,3,h])$

wxplot2d( [ [discrete,numsolrk2] ],[x,x0,x1] ,[xlabel, "t"], [ylabel, "v"],
[color,blue,red]);