// Ejemplito de metodos numericos para ecuaciones diferenciales ordinarias
//
//
//     Atractor de Lorenz
//
//     El atractor de Lorenz es la representación de un sistema de ecuaciones diferenciales
//     autónomo (no depende explícitamente del tiempo), un sistema dinámico, que en principio
//     su autor lo propuso para tratar de comprender los fenómenos meteorológicos. Se
//     encontró con que a pesar de que los valores generados nunca se repiten y de que las
//     condiciones iniciales pueden hacer variar completamente los valores generados.(de ahí su
//     nombre de atractor extraño o caótico, y también que los que predicen el tiempo se
//     equivoquen tanto), el atractor toma una forma única y parece conservar cierto orden, sus
//     infinitas trayectorias nunca se cortan (pues eso implicaría que entraríamos en un ciclo
//     periódico en un sistema autónomo), y también podemos comprobar cómo este objeto
//     presenta la autosimilitud característica de los fractales. Esta extrema sensibilidad a las
//     condiciones iniciales fue lo que originó la frase del llamado Efecto Mariposa: "...el aleteo
//     de una mariposa en Australia sería suficiente para provocar un tornado en Tejas..."
//     Conviene señalar que las condiciones necesarias para que exista caos en un sistema de
//     ecuaciones diferenciales autónomo son: deben haber al menos tres ecuaciones
//     diferenciales y al menos tres variables y al menos alguna no linealidad.
//     El modelo atmosférico que utilizó Lorenz consiste en una atmósfera bidimensional
//     rectangular, cuyo extremo inferior está a una temperatura mayor que el superior. De esta
//     manera el aire caliente subirá y el aire frío bajará creándose corrientes que harán un
//     intercambio de calor por convección.
//
//     Las ecuaciones que describen este proceso son:
//     dx/dt = s(y-x)
//     dy/dt = rx - y - xz
//     dz/dt = xy - bz
//     En donde las variables, que únicamente dependen del tiempo son:
//
//     "x" representa el flujo convectivo
//     "y" es la distribución de temperaturas horizontal
//     "z" es la distribución de temperaturas vertical
//     Además tenemos de los tres parámetros que intervienen en las ecuaciones:
//
//     "s" cociente entre la viscosidad y la conductividad térmica
//     "r" la diferencia de temperaturas entre la capas inferior y superior
//     "b" el cociente entre la altura y el ancho de nuestro rectángulo
//
//
// EQUATIONS
//     x'=s*(y-x)
//     y'=r*x-y-x*z
//     z'=x*y-b*z
//
//     s=10
//     r=28
//     b=2.666667

 

function xf1=f1(t,x,y,z)

  // Parametros
  s=10;
  r=28;
  b=2.666667;

  // Ecuacion
  xf1=s*(y-x);

endfunction


function yf2=f2(t,x,y,z)

  // Parametros
  s=10;
  r=28;
  b=2.666667;

  // Ecuacion
  yf2=r*x-y-x*z;

endfunction


function zf3=f3(t,x,y,z)

  // Parametros
  s=10;
  r=28;
  b=2.666667;

  // Ecuacion
  zf3=x*y-b*z;

endfunction



N  = 20000; // Particion 
t0 = 0.0; // Tiempo inicial
t1 = 60.0; // Tiempo final
cix = 3.0; // Condicion inicial al tiempo t0
ciy = 3.0;
ciz = 3.0;

h = (t1 - t0) / N; // incremento en el tiempo
t = t0; // Asigna el tiempo inicial

// Asigna la condicion inicial 
e1 = cix;
e2 = ciy;
e3 = ciz;

// Solo para graficar
xgrid
X = zeros(N+1,1);
Y = zeros(N+1,1);
Z = zeros(N+1,1);
X(1) = cix;
Y(1) = ciy;
Z(1) = ciz;

// Calculo de las soluciones
for i=1:N 
   tiempo = t+h; // Solo para la visualizacion
 
   X(i+1) = e1 + h * f1(t,e1,e2,e3);  // Euler
   Y(i+1) = e2 + h * f2(t,e1,e2,e3); 
   Z(i+1) = e3 + h * f3(t,e1,e2,e3); 

   e1 = X(i+1); // evolucion para graficar
   e2 = Y(i+1); 
   e3 = Z(i+1);

   t = t0 + i * h; // Calcula en paso en el tiempo
end

// Grafica la solucion
plot3d(X,Y,Z)  // Calculada