Sistema de ecuaciones diferenciales de segundo orden

Algunas páginas del Curso Interactivo de Física describen situaciones físicas en las que se aplica el procedimiento numérico de Runge-Kutta para resolver un sistema de dos ecuaciones diferenciales de segundo orden, que puede generalizarse a un sistema de n ecuaciones diferenciales de segundo orden.

Descripción

Sea el sistema de dos ecuaciones diferenciales de segundo orden

$\frac{d^{2} x}{d t^{2}} = f (x, y, v_{x}, v_{y}, t) \frac{d^{2} y}{d t^{2}} = g (x, y, v_{x}, v_{y}, t)$

con las condiciones iniciales

$x (t_{0}) = x_{0} {(\frac{d x}{d t})}_{t_{0}} = v_{x 0} y (t_{0}) = y_{0} {(\frac{d y}{d t})}_{t_{0}} = v_{y 0}$

Este sistema, se puede transformar en un sistema equivalente formado por cuatro ecuaciones diferenciales de primer orden. Aplicando dos veces el esquema descrito para una ecuación diferencial de segundo orden, obtenemos el esquema descrito en las siguientes tablas

$\frac{d x}{d t} = v_{x}$	$\frac{d v_{x}}{d t} = f (x, y, v_{x}, v_{y}, t)$
$\begin{array}{l} k_{1} = h v_{x} \\ k_{2} = h (v_{x} + \frac{1}{2} l_{1}) \\ k_{3} = h (v_{x} + \frac{1}{2} l_{2}) \\ k_{4} = h (v_{x} + l_{3}) \end{array}$	$\begin{array}{l} l_{1} = h \cdot f (x, y, v_{x}, v_{y}, t) \\ l_{2} = h \cdot f (x + \frac{1}{2} k_{1}, y + \frac{1}{2} q_{1}, v_{x} + \frac{1}{2} l_{1}, v_{y} + \frac{1}{2} m_{1}, t + \frac{1}{2} h) \\ l_{3} = h \cdot f (x + \frac{1}{2} k_{2}, y + \frac{1}{2} q_{2}, v_{x} + \frac{1}{2} l_{2}, v_{y} + \frac{1}{2} m_{2}, t + \frac{1}{2} h) \\ l_{4} = h \cdot f (x + k_{3}, y + q_{3}, v_{x} + l_{3}, v_{y} + m_{3}, t + h) \end{array}$
$x (t + h) = x (t) + \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$	$v_{x} (t + h) = v_{x} (t) + \frac{1}{6} (l_{1} + 2 l_{2} + 2 l_{3} + l_{4})$

$\frac{d y}{d t} = v_{y}$	$\frac{d v_{y}}{d t} = g (x, y, v_{x}, v_{y}, t)$
$\begin{array}{l} q_{1} = h v_{y} \\ q_{2} = h (v_{y} + \frac{1}{2} m_{1}) \\ q_{3} = h (v_{y} + \frac{1}{2} m_{2}) \\ q_{4} = h (v_{y} + m_{3}) \end{array}$	$\begin{array}{l} m_{1} = h \cdot g (x, y, v_{x}, v_{y}, t) \\ m_{2} = h \cdot g (x + \frac{1}{2} k_{1}, y + \frac{1}{2} q_{1}, v_{x} + \frac{1}{2} l_{1}, v_{y} + \frac{1}{2} m_{1}, t + \frac{1}{2} h) \\ m_{3} = h \cdot g (x + \frac{1}{2} k_{2}, y + \frac{1}{2} q_{2}, v_{x} + \frac{1}{2} l_{2}, v_{y} + \frac{1}{2} m_{2}, t + \frac{1}{2} h) \\ m_{4} = h \cdot g (x + k_{3}, y + q_{3}, v_{x} + l_{3}, v_{y} + m_{3}, t + h) \end{array}$
$y (t + h) = y (t) + \frac{1}{6} (q_{1} + 2 q_{2} + 2 q_{3} + q_{4})$	$v_{y} (t + h) = v_{y} (t) + \frac{1}{6} (m_{1} + 2 m_{2} + 2 m_{3} + m_{4})$

La codificación de la función que resuelve el sistema de dos ecuaciones diferenciales de segundo orden, repite dos veces, el código de la función resolver correspondiente a una ecuación diferencial de segundo orden, ya que hay que calcular cuatro pares de parámetros más, los correspondientes a la segunda ecuación diferencial de segundo orden.

La jerarquía de clases

De modo análogo a los apartados anteriores, estableceremos una jerarquía de clases, la clase base abstracta define el procedimiento numérico y la clase derivada define el sistema físico particular.

public abstract class RungeKutta {
      double h;
  public RungeKutta(double h) {
     this.h=h;
  }
  public void setPaso(double h){
        this.h=h;
  }
   public void resolver(Estado e){
//variables auxiliares
	    double k1, k2, k3, k4;
	    double l1, l2, l3, l4;
	    double q1, q2, q3, q4;
	    double m1, m2, m3, m4;
//estado en el instante t
	    double x=e.x;
	    double y=e.y;
	    double vx=e.vx;
	    double vy=e.vy;
	    double t=e.t;

		    k1=h*vx;
		    l1=h*f(x, y, vx, vy, t);
		    q1=h*vy;
		    m1=h*g(x, y, vx, vy, t);

		    k2=h*(vx+l1/2);
		    l2=h*f(x+k1/2, y+q1/2, vx+l1/2, vy+m1/2, t+h/2);
		    q2=h*(vy+m1/2);
		    m2=h*g(x+k1/2, y+q1/2, vx+l1/2, vy+m1/2, t+h/2);

		    k3=h*(vx+l2/2);
		    l3=h*f(x+k2/2, y+q2/2, vx+l2/2, vy+m2/2, t+h/2);
		    q3=h*(vy+m2/2);
		    m3=h*g(x+k2/2, y+q2/2, vx+l2/2, vy+m2/2, t+h/2);

		    k4=h*(vx+l3);
		    l4=h*f(x+k3, y+q3, vx+l3, vy+m3, t+h);
		    q4=h*(vy+m3);
		    m4=h*g(x+k3, y+q3, vx+l3, vy+m3, t+h);
    

		    x+=(k1+2*k2+2*k3+k4)/6;
		    vx+=(l1+2*l2+2*l3+l4)/6;
		    y+=(q1+2*q2+2*q3+q4)/6;
		    vy+=(m1+2*m2+2*m3+m4)/6;
        t+=h;
//estado en el instante t+h
	    e.x=x;
	    e.y=y;
	    e.vx=vx;
	    e.vy=vy;
	    e.t=t;
    }
    abstract public double f(double x, double y, double vx, double vy, double t);
    abstract public double g(double x, double y, double vx, double vy, double t);
}

El estado del sistema vendrá dado por cuatro números x, y, v_x, v_y y t, miembros dato de la clase que hemos denominado Estado.

public class Estado {
    public double t;
    public double x;
    public double y;
    public double vx;
    public double vy;
    public Estado(double t, double x, double y, double vx, double vy) {
        this.t=t;
        this.x=x;
        this.y=y;
        this.vx=vx;
        this.vy=vy;
    }
}

Ejemplo

Uno de los ejemplos más interesantes de resolución de un sistema de ecuaciones diferenciales de segundo orden es la descripción del movimiento de los planetas, el cual tiene una solución analítica sencilla en coordenadas polares. La trayectoria seguida por un planeta es una cónica, una elipse en particular, en uno de cuyos focos está el centro fijo de fuerzas, el Sol. En la figura, se muestra la fuerza que ejerce el Sol sobre un planeta, inversamente proporcional al cuadrado de las distancias que separan sus centros, y también se muestran sus componentes rectangulares

Teniendo en cuenta que la fuerza que ejerce el sol sobre un planeta viene descrita por la ley de la Gravitación Universal

$F = G \frac{M m}{r^{2}}$

donde M es la masa del Sol, m la masa del planeta y r la distancia entre el centro del Sol y del planeta. Las componentes de la aceleración del planeta serán

$\begin{array}{l} a_{x} = \frac{- F_{x}}{m} = \frac{- F \cos θ}{m} = - G \frac{M}{r^{3}} x \\ a_{y} = \frac{- F_{y}}{m} = \frac{- F \sin θ}{m} = - G \frac{M}{r^{3}} y \end{array}$

Uno de los problemas del tratamiento numérico con ordenador, es la de reducir el problema a números simples e inteligibles por el usuario de un vistazo. Las masa de los planetas y del Sol son números muy grandes: la masa de la Tierra es 5.98 10²⁴ kg., y 1.98 10³⁰ kg. la del Sol. La distancia media entre la Tierra y el Sol es también muy grande 1.49 10¹¹ m. y la constante G es muy pequeña 6.67 10^-11 en el Sistema Internacional de Unidades. Podemos simplificar el problema numérico, refiriéndonos a un hipotético Sol cuya masa sea tal que el producto GM=1 o bien, que se ha cambiado la escala de los tiempos de modo que se cumpla esta condición. Teniendo en cuenta que la aceleración es la derivada segunda de la posición, el sistema de dos ecuaciones diferenciales de segundo orden quedará

$\frac{d^{2} x}{d t^{2}} = - \frac{x}{r^{3}} \frac{d^{2} y}{d t^{2}} = - \frac{y}{r^{3}}$

En la clase derivada Sistema redefinimos las funciones f y g.

public class Sistema extends RungeKutta{
    public Sistema(double h){
      super(h);
    }
    public double f(double x, double y, double vx, double vy, double t){
        double r=Math.sqrt(x*x+y*y);
        return (-x/(r*r*r));
    }
    public double g(double x, double y, double vx, double vy, double t){
        double r=Math.sqrt(x*x+y*y);
        return (-y/(r*r*r));
    }
}

Los pasos para usar la clase Sistema son los siguientes

Establecer la magnitud del paso h

Introducir el estado inicial de la partícula (su posición y su velocidad inicial)

Establecer el tiempo final t, en el que deseamos calcular el nuevo estado del planeta.

Crear un objeto de la clase Sistema y llamar a la función miembro resolver.

Mostrar la posición y velocidad final del planeta

public class Aplicacion {
  public static void main(String[] args) {
      double h=0.1;       //paso
//Situación inicial
        double x0=1.0;
        double vy0=1.2;
        Estado e=new Estado(0.0, x0, 0.0, 0.0, vy0);

        double t=10.0;      //resolver la e. d. hasta este instante
        Sistema sis=new Sistema(h);
        for(int i=0; i<t/h; i++){
            sis.resolver(e);
        }
        System.out.println("posición aproximada x "+e.x);
        System.out.println("posición aproximada y "+e.y);
        System.out.println("velocidad aproximada vx "+e.vx);
        System.out.println("velocidad aproximada vy "+e.vy);
  }
}