Resolución de las ecuaciones del Modelo de Lorenz

    Una vez simplificado el Modelo de Lorenz, se llega a las expresiones siguientes:

{1},{2} y {3}

    Donde Pr, r y b son los parámetros de las ecuaciones.

    Estas ecuaciones se pueden resolver numéricamente empleando diferentes métodos.

    El más simple es emplear el método de Euler y uno de los que quizás sean más efectivos sea el método de Runge - Kutta. Los dos consisten en realizar una aproximación para el cálculo de la derivada.

    El método de Euler consiste en:

    Si tenemos una ecuación del tipo y'[x_i] = F[x_i,y[x_i]], entonces dada una conción inicial, podremos calcular la evolución de y[x], tomando y_i+1[x_i+1]=y_i[x_i]+F[x_i,y[x_i]] Dx_i, siendo Dx_i =x_i+1-x_i y tan pequeño como podamos, por tanto cuanto más pequeño sea x_i, más preciso es el resultado. Evidentemente hay que dar la C.I. Este método es de primer orden.

    Bien, en VB para resolver las ecuaciones, el código que hay que escribir es el siguiente, en su expresión más sencilla:

    Para el Método de Euler:

For i = 1 To N

    dx = Pr * (y - x) * deltaT

    dy = (-x * z + r * x - y) * deltaT

    dz = (x * y - b * z) * deltaT

    x = x + dx

    y = y + dy

    z = z + dz

    Print #1, x, y, z

Next i

    Donde se han definido las variables todas como de doble precisión, salvo i y N que son de enteras de doble longitud.

    El contenido de estas variables, se corresponden fácilmente con sus nombres, así:

        x,y y z representan a las correspondientes variables del problema

        dx,dy y dz representan a los incrementos de las variables de x, y y z, en cada paso de integración.

      Pr,r y b representan a los parámetros del sistema de ecuaciones.

      deltaT Representa a al incremento de tiempo

      i Es un iterador o un contador y N es el Nº total de pasos de integración que se pretende realizar.

    Con la instrucción con la instrucción Print se guarda el resultado en un fichero.

    Antes de ejecutar estas instrucciones, hay que introducir el valor de todas las constantes (Pr, r, b, deltaT y N) y la condición inicial de x,y, y z.

    Evidentemente nos podemos extender en todos los detalles de introducción de la variables, en el programa de definición de las mismas, etc. Pero se puede considerar innecesario. También se puede pensar en otras formas de escribir el programa con la idea de que sea más bonito, etc. pero también se puede considerar que es innecesario, en cualquier caso se puede consultar el apéndice donde se incluye el listado completo del programa que se adjunta. Donde el método lo hemos implementado como una clase en Visual Basic, al estilo de lo que se hace en C++. Nos referiremos a los resultados obtenidos con el método de cálculo.

    En primer lugar como se vio en la teoría, uno de los valores de los parámetros más interesantes para resolver estas ecuaciones son:

        Que b varíe entre: 8/3 y 4.

        Que Pr varíe entre 10 y 16.

        Que r se libre, mayor que 0.

    Primero nos fijaremos en los valores b = 4 y Pr = 16. Al pedir el programa que resuelva el problema, nos encontramos con circunstancias muy similares pero con soluciones bien diferentes dependiendo de: Las condiciones iniciales y el valor de r.

    De esta forma se puede reflejar en una tabla algunos de los resultados obtenidos:

    Como se puede observar, dadas las condiciones iniciales, para un cierto rango de valores, obtenemos una solución convergente, mientras que para otro rango, no se da ninguna solución convergente, también podemos ver que a partir de cierto valor de r la solución es convergente a un ciclo límite, aunque si cambiamos el Dt empleado para resolver las ecuaciones, podemos encontrar que determinadas soluciones que parecían ser un ciclo límite no lo son. Esto no advierte de que hay que tener mucho cuidado a la hora de decir cual es la solución concreta, pues el sistema tiene cierta sensibilidad al método de cálculo.

    Veamos pues que por otra parte que sucede cuando se cambian las condiciones iniciales:

    Como se puede observar en la tabla, se obtienen dos posibles soluciones estacionarias, que se deben corresponder precisamente con la solución algebraica que se puede obtener fácilmente:

{4},{5},{6} y {7}

    Donde se puede observar que para r < 1 no hay solución.

    De aquí se obtiene que:

x = y = 6

z = 9

    Lo que está de acuerdo perfectamente con el resultado numérico.

    Ahora, evidentemente, podemos construir un programa que nos de las soluciones estacionarias para los diferentes valores de r. Por en vez de eso vamos a fijarnos en otras cuestiones, ¿Que es lo que pasa cuando no encontramos con una solución que no converge o que parece un ciclo límite?. Pues bien, veamos:

    Vamos a emplear en esta ocasión los parámetros b = 8/3 y Pr = 10 y r variable, donde vamos a ver que es lo que sucede para valores de r diferentes y en especial que es lo que pasa en los puntos r=24,74 valor para el cual según la teoría se debe de dar una transformación de Hopf y a r = 28 donde aparece lo que se conoce como atractor de Lorenz:

    Ahora, como se ve nuevamente para ciertos valores, la solución par tiempos grandes, puede ser una solución estacionaria, y para otros valores no. En concreto para r = 28 no a aparece solución estacionaria y si se hace una representación de los datos obtenidos, se puede obtener la siguiente gráfica:

    Donde se observa, que el sistema, evoluciona de una forma un tanto peculiar entorno a las dos soluciones estacionarias predichas en las expresiones {4},{5},{6} y {7}. En esta situación se da lo que llamamos un comportamiento caótico. Esto se puede confirmar al comprobar que el espectro de frecuencias de alguno de los parámetros que conforman la trayectoria es continuo, o al menos no es discreto.

    Así, representando por ejemplo el parámetro x vs t, obtenemos la siguiente gráfica:

Si tomamos una muestra lo suficientemente amplia, obtendremos lo siguiente:

    Aunque aquí solo se visualiza una gráfica con 10⁴ datos, en realidad la muestra tomada es de:482192 datos. Bien, si a dicha muestra le realizamos la fast Fourier transform, y después de ordenarla y de representar su valor absoluto, se obtiene la siguiente gráfica:

    Ampliando sobre la zona de interés:

    Como podemos ver el espectro de frecuencias para la componente x es tal, que no es discreto, confirmandonos la idea de movimiento caótico para el sistema que estamos estudiando.

    Por otra parte las soluciones estacionarias predichas por las ecuaciones {4} a la {7}, siendo r=28, b= 8/3 son X=Y=±8.485281374 y Z = 27. Si tomamos estos parámetros como condiciones iniciales, ¿Que es lo que sucede? Pues sorprendentemente el sistema se queda en dicho punto, como se puede comprobar numéricamente, pero si tomamos un punto próximo al mismo, este no converge, sino que se queda oscilando entorno al mismo. De esa forma se ha obtenido la siguiente gráfica:

    Donde las condiciones iniciales han sido: X = Y = 8,48528, Z = 27,0001,y Dt = 0,0000001.

    Cuando pasa suficiente tiempo, el sistema evoluciona hasta dar una nueva figura como la anterior del atractor de Lorenz:

    Desde otra perspectiva:

    Evidentemente, no podemos asegurar este comportamiento para todas las condiciones iniciales entorno a las soluciones estacionarias ya que tendríamos que hacer un análisis exhaustivo del comportamiento en dicho entorno.

    Por otra parte si nos fijamos en lo que sucede para valores como r = 10, obtenemos la gráfica ya esperada según se ha comentado anteriormente, tiende a una situación estacionaria.

    Ahora, si nos fijamos en r = 100, hemos indicado en la tabla que el sistema tiende a un ciclo límite como se puede ver en la representación gráfica:

    Observese que aquí, se han tomado 110800 datos, de forma que si hacemos una representación gráfica de t vs x, se obtiene:

    Donde se puede observar a simple vista que, los valores de X son periódicos, haciendo el análisis de Furier correspondiente, se obtiene:

    Como podemos ver no tiene nada que ver con el encontrado en el caso de r = 28.

    Empleando el método de Runge - Kutta de cuarto orden, el algoritmo para encontrar la solución es mucho más complicado, para él, dada una C.I. la solución se obtiene tomando:

y_i+1[x_i+1] = y_i[x_i]+1/6(Dy₀+2Dy₁+2Dy₃+Dy₃)

    de tal forma que:

Dy₀ = F[x_i,y_i]Dx

Dy₁ = F[x_i+Dx/2, y_i+1/2Dy₀/2] Dx

Dy₂ = F[x_i+Dx/2, y_i+Dy_1/2] Dx

Dy₃ = F[x_i+1, y_i+Dy₂]Dx

    Según se indica por ejemplo en el libro: Introducción a los métodos numéricos. A.A.Samarski. Ed. Mir Moscú .

    Esto lo hemos implementado al igual que en caso anterior como una clase en V.B. el código que hay que escribir es un poco más complicado en este caso. En el módulo de clase se incluyen todas las funciones y propiedades para la solución del problema.

'Class Name:LRK

Option Explicit

Private m_x As Double

Private m_y As Double

Private m_z As Double

Private m_dx As Double

Private m_dy As Double

Private m_dz As Double

Private m_Pr As Double

Private m_b As Double

Private m_r As Double

Private m_iteracion As Long

Private m_t As Double

Private m_Dt As Double

Public Sub LRK_inicialize()

    m_x = 10#

    m_y = 10#

    m_z = 10#

    m_dx = 0.001

    m_dy = 0.001

    m_dz = 0.001

    m_Pr = 16

    m_b = 8 / 3

    m_r = 28

    m_iteracion = 0#

    m_Dt = 0.001

    m_t = 0#

End Sub

Public Property Get t() As Double

    t = m_t

End Property

Public Property Let t(ByVal vNewValue As Double)

    If vNewValue < 0# Then vNewValue = 0#

    MsgBox ("Se ha cambiado el valor de t = " & m_t & " por t = 0")

    m_t = vNewValue

End Property

Public Property Get dt() As Double

    dt = m_Dt

End Property

Public Property Let dt(ByVal vNewValue As Double)

    m_Dt = vNewValue

End Property

Public Property Get iteracion() As Long

    iteracion = m_iteracion

End Property

Public Property Let iteracion(ByVal vNewValue As Long)

    m_iteracion = vNewValue

End Property

Public Property Get X() As Double

    X = m_x

End Property

Public Property Let X(ByVal vNewValue As Double)

    m_x = vNewValue

End Property

Public Property Get Y() As Double

    Y = m_y

End Property

Public Property Let Y(ByVal vNewValue As Double)

    m_y = vNewValue

End Property

Public Property Get Z() As Double

    Z = m_z

End Property

Public Property Let Z(ByVal vNewValue As Double)

    m_z = vNewValue

End Property

Public Property Get dX() As Double

    dX = m_dx

End Property

Public Property Let dX(ByVal vNewValue As Double)

    m_dx = vNewValue

End Property

Public Property Get dY() As Double

    dY = m_dy

End Property

Public Property Let dY(ByVal vNewValue As Double)

    m_dy = vNewValue

End Property

Public Property Get dZ() As Double

    dZ = m_dz

End Property

Public Property Let dZ(ByVal vNewValue As Double)

    m_dz = vNewValue

End Property

Public Property Get r() As Double

    r = m_r

End Property

Public Property Let r(ByVal vNewValue As Double)

    m_r = vNewValue

End Property

Public Property Get b() As Double

    b = m_b

End Property

Public Property Let b(ByVal vNewValue As Double)

    m_b = vNewValue

End Property

Public Property Get Pr() As Double

    Pr = m_Pr

End Property

Public Property Let Pr(ByVal vNewValue As Double)

    m_Pr = vNewValue

End Property

Public Function putx(x_int As Double) As Integer

    m_x = x_int

putx = 1

End Function

Private Function Dfx(X As Double, Y As Double) As Double

    Dfx = m_Pr * (Y - X)

End Function

Private Function Dfy(X As Double, Y As Double, Z As Double) As Double

    Dfy = -X * Z + m_r * X - Y

End Function

Private Function Dfz(X As Double, Y As Double, Z As Double) As Double

    Dfz = X * Y - m_b * Z

End Function

Private Function RKx() As Double

    Dim Dx0 As Double, Dx1 As Double, Dx2 As Double, Dx3 As Double

    Dx0 = Dfx(m_x, m_y) * m_Dt

    Dx1 = Dfx(m_x + 0.5 * Dx0, m_y) * m_Dt

    Dx2 = Dfx(m_x + 0.5 * Dx1, m_y) * m_Dt

    Dx3 = Dfx(m_x + Dx2, m_y) * m_Dt

    RKx = m_x + 1 / 6 * (Dx0 + 2 * (Dx1 + Dx2) + Dx3)

End Function

Private Function RKy() As Double

    Dim Dy0 As Double, Dy1 As Double, Dy2 As Double, Dy3 As Double

    Dy0 = Dfy(m_x, m_y, m_z) * m_Dt

    Dy1 = Dfy(m_x, m_y + 0.5 * Dy0, m_z) * m_Dt

    Dy2 = Dfy(m_x, m_y + 0.5 * Dy1, m_z) * m_Dt

    Dy3 = Dfy(m_x, m_y + Dy2, m_z) * m_Dt

    RKy = m_y + 1 / 6 * (Dy0 + 2 * (Dy1 + Dy2) + Dy3)

End Function

Private Function RKz() As Double

    Dim Dz0 As Double, Dz1 As Double, Dz2 As Double, Dz3 As Double

    Dz0 = Dfz(m_x, m_y, m_z) * m_Dt

    Dz1 = Dfz(m_x, m_y, m_z + 0.5 * Dz0) * m_Dt

    Dz2 = Dfz(m_x, m_y, m_z + 0.5 * Dz1) * m_Dt

    Dz3 = Dfz(m_x, m_y, m_z + Dz2) * m_Dt

    RKz = m_z + 1 / 6 * (Dz0 + 2 * (Dz1 + Dz2) + Dz3)

End Function

Public Function paso(ByRef x_paso As Double, ByRef y_paso As Double, ByRef z_paso As Double) As Long

    x_paso = RKx()

    y_paso = RKy()

    z_paso = RKz()

    m_x = x_paso

    m_y = y_paso

    m_z = z_paso

    m_iteracion = m_iteracion + 1

    m_t = m_iteracion * m_Dt

    paso = m_iteracion

End Function

    Este algoritmo también lo podemos encontrar, el las librerías del compilador de FORTRAN: Fortran PowerStation 4.0, donde se implementan las subrutinas de Runge-Kutta de 4 orden, y también las de 6º orden.

Ejemplo de subrutina de la IMSL:

SUBROUTINE rk4(y,dydx,n,x,h,yout,derivs)

    INTEGER n,NMAX

    REAL h,x,dydx(n),y(n),yout(n)

    EXTERNAL derivs

    PARAMETER (NMAX=50)

    INTEGER i

    REAL h6,hh,xh,dym(NMAX),dyt(NMAX),yt(NMAX)

    hh=h*0.5

    h6=h/6.

    xh=x+hh

    do 11 i=1,n

        yt(i)=y(i)+hh*dydx(i)

11 continue

    call derivs(xh,yt,dyt)

    do 12 i=1,n

        yt(i)=y(i)+hh*dyt(i)

12 continue

    call derivs(xh,yt,dym)

    do 13 i=1,n

        yt(i)=y(i)+h*dym(i)

        dym(i)=dyt(i)+dym(i)

13   continue

    call derivs(x+h,yt,dyt)

    do 14 i=1,n

        yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))

14  continue

    return

END

SUBROUTINE derivs(x,y,dydx)

    REAL x,y(*),dydx(*)

    dydx(1)=-.013*y(1)-1000.*y(1)*y(3)

    dydx(2)=-2500.*y(2)*y(3)

    dydx(3)=-.013*y(1)-1000.*y(1)*y(3)-2500.*y(2)*y(3)

return

END

    También, en la librería llamada IMSL encontramos varias subrutinas, escritas para el método de R-K entre ellas se encuentran las subrutinas IVPRK/DIVPRK para cuarto y sexto orden con simple y doble precisión. Y la subrutina, DVERK de 5 y 6º orden.

    Bien, si decidimos emplear la subrutina DVERK, esta, se emplea utilizando la llamada:

    Call DVERK (n,fnc,x,y,xend,tol,ind,c,nw,w,ier)

    Donde los datos fundamentales que debemos de conocer son:

        n Número de ecuaciones diferenciales a resolver.

        Fnc es el nombre de la subrutina que implementa las ecuaciones diferenciales a resolver y debe tener la forma:

        Subrutine fnc(n,x,y,yprime)

    Donde: n es el número de ecuaciones, x la variable independiente, y la dependiente y yprime la derivada.

    X es la variable independiente.

    Y es la dependiente.

    Xend es el valor de x para el que la subrutina devuelve la y.

    Tol es el valor de la precisión con la que se pretende obtener la solución.

    Ind y C son indicadores.

    Nw y W son un indice y una matriz para trabajar con ella.

    Ier da valores sobre posibles errores en el cálculo.

    Con esta información se puede generar un programa simple, que resuelva las ecuaciones para unos determinados valores de los parámetros, hemos intentado hacerlo lo más simple posible, escribiendo el siguiente código:

    Program atractor

    integer N,Ind,NW,IER,K

    double precision Y(3),C(24),W(3,9),X,Tol,Xend

    external fnc1

    NW =3

    N=3

    X=0.0

    xend=0.0

    Y(1)=10.0

    Y(2)=10.0

    y(3)=10.0

    Tol=0.00001

    IND = 1

    write(1,*)"datos=["

    do 10 k = 1, 100000

        xend=xend+0.0001

        call dverk(n,fnc1,x,y,xend,tol,ind,c,nw,w,ier)

        if (ind.lt.0.or.ier.gt.0) go to 20

        print*, x,y(1),y(2),y(3)

        write(1,*) y(1),y(2),y(3)

10 continue

    write(1,*)"];"

    close(1)

    read(*,*)

    stop

20 continue

    write(1,*)"];"

    close(1)

    read(*,*)

    stop

    end

    subroutine fnc1(n,x,y,yprime)

        integer n

        double precision y(n), yprime(n),x

        double precision b

        b=8.0/3.0

                    yprime(1) = 10.0*(Y(1)-y(2))

                    yprime(2) = -Y(1)*Y(3)+28*Y(1)-y(2)

                    yprime(3) = y(1)*y(2)-b*y(3)

        x=x

        return

    end

    Como se puede observar al ejecutar el programa resultante, se obtiene la solución para el atractor de Lorenz para el que hemos hecho el análisis en frecuencias, anteriormente empleando el otro método.

    En este caso, la solución que se obtiene no se parece demasiado a lo obtenido con Euler, como se puede observar:

    En cambio si empleamos la clase creada para la resolución de este problema en VB, el resultado es bastante similar:

    Y donde se espera que para las condiciones iniciales, y los parámetros dados, el resultado sea mejor.

    Este fallo en el cálculo de RK para 5 y 6ºorden puede deberse a muchos motivos, entre otros como es lógico a no haber implementado correctamente las ecuaciones en el programa FORTAN.

    Lógicamente la forma correcta de comparar estos resultados es por ejemplo representando en una misma gráfica los dos resultados obtenidos:

    Como se puede observar en la gráfica, los dos métodos dan un resultado identico para aproximadamente los primeros mil valores, mientras que la solución deja de parecerse a partir de dichos datos. Esto se puede visualizar mejor en la siguiente gráfica.

    Esto lógicamente, indica lo que ya sabíamos sobre estas ecuaciones y es que cuando la solución no es convergente a unos valores estacionarios, la evolución es muy sensible a las condiciones iniciales, en este caso se ha puesto de manifiesto a partir del método de cálculo, hemos de suponer que de los dos métodos el que está pintado en Azul, correspondiente al de RK, debería de dar una mejor solución, pero dada las características de estas ecuaciones, realmente no se puede afirmar nada.

nº:

        Juan Ignacio Lagares González.

Julio de 2001.

Regresar