Método de Runge-Kutta orden 4 para sistemas y Ecuaciones de Lorenz

El método de Runge-Kutta orden 4 para sistemas es una extensión del método numérico usado para resolver una ecuación de la forma \frac{dy}{dx}=f(x,y) al caso de un sistema de n de tales ecuaciones, es decir permite resolver:

     \begin{align*} \frac{du_1}{dt} &= f_1(t,u_1,\cdots,u_n) \\ \frac{du_2}{dt} &= f_2(t,u_1,\cdots,u_n) \\  \vdots \quad  &= \quad \vdots \\ \frac{du_n}{dt}  &= f_n(t,u_1,\cdots,u_n) \end{align*}

Si \mathbf{u}=(u_1,u_2,\cdots,u_n), entonces el sistema anterior se puede representar vectorialmente como \mathbf{u}' = \mathbf{f}(t,\mathbf{u}), de forma que si se aproxima \mathbf{u}(t_j) con \mathbf{w}_j  (con t_0 \leq t_j \leq t_n),  entonces podemos escribir el método como:

 \mathbf{w}_{j+1} = \mathbf{w}_j+\frac{h}{6}(\mathbf{k}_1+2 \mathbf{k}_2 + 2 \mathbf{k}_3 + \mathbf{k}_4)

donde

  •  \mathbf{k}_1 = \mathbf{f}(t_j, \mathbf{w}_j)
  •  \mathbf{k}_2 = \mathbf{f}(t_j+\frac{h}{2}, \mathbf{w}_j+\frac{h}{2} \mathbf{k}_1)
  •  \mathbf{k}_3 = \mathbf{f}(t_j+\frac{h}{2}, \mathbf{w}_j+\frac{h}{2} \mathbf{k}_2)
  •  \mathbf{k}_4 = \mathbf{f}(t_j+h, \mathbf{w}_j+h \mathbf{k}_3)

para t_j = t_0 + j h y una condición inicial \mathbf{w}(t_0)=\mathbf{w}_0. En el archivo ordinarias.rkt, se implementa en Scheme/DrRacket (versión 5.3)  este método.

Ejemplo: veamos el caso de las ecuaciones clásicas de Lorenz (ver atractor de Lorenz):

     \begin{align*} \frac{dx}{dt} &= \sigma (y-x) \\ \frac{dy}{dt} &= \rho x - y -xz \\  \frac{dz}{dt} &=xy-\beta z \end{align*}

Estas ecuaciones con los valores \sigma =10, \beta = \frac{8}{3} y \rho = 28, generan la trayectoria de solución indicada en la figura, cuando se inicia en el punto (10,7,7).

Le invitamos a explorar el código en lorenz.rkt (y archivos auxiliares: ordinarias.rkt, derivaciones.rkt, vectoriales.rkt, basicas.rkt, para Scheme/DrRacket [versión 5.3]). Las instrucciones básicas utilizadas son:

(define sigma 10)
(define beta 8/3)
(define rho 28)

(define fs (list
            (lambda(t x y z)
              (* 10 (- y x)))
            (lambda(t x y z)
              (- (* 28 x) y (* x z)))
            (lambda(t x y z)
              (- (* x y) (* 8/3 z)))))
(define sols
  (srk4-h fs 0 30 '(10 7 7) 0.005))

Note que (en este caso) se incluye en la lista de puntos solución (j,t_j,x(t_j),y(t_j),z(t_j)) donde t_j=t_0 + j h para h=0.005

Observación: si se desea resolver una ecuación diferencial de orden superior, es suficiente utilizar cambios de variable y convertir dicha ecuación a un sistema equivalente de ecuaciones diferenciales de primer orden. Por ejemplo la ecuación

\frac{d^2x}{dt^2}=-\frac{k}{m}x

asociada a un resorte, se puede convertir a un sistema de dos ED de primer orden, introduciendo la variable u=\frac{dx}{dt}, con lo cual \frac{du}{dt}=\frac{d^2x}{dt^2} y entonces el sistema equivalente es:

     \begin{align*} \frac{du}{dt} &= -\frac{k}{m}x \\ \frac{dx}{dt} &= u \end{align*}

Nota: se corrigió en mordinarias.rkt la función srk4-h (y en consecuencia la gráfica de Lorenz), favor de actualizar a versión 0.4; gracias y una disculpa [2011.12.23]