Siguiente Anterior

Sistema de ecuaciones diferenciales de primer orden

El procedimiento de Runge-Kutta es igualmente efectivo en la resolución de un sistema de ecuaciones diferenciales de primer orden.

dx dt =f(x,y,t) dy dt =g( x,y,t)

El procedimiento de aplicación del método de Runge-Kutta a cada una de las ecuaciones diferenciales, con las condiciones iniciales

x( t 0 )= x 0 y( t 0 )= y 0

se esquematiza en la tabla adjunta. Como vemos además de los cuatro números k1, k2, k3, k4 para la primera ecuación diferencial precisamos otros cuatro números l1, l2, l3, l4 para la segunda ecuación diferencial. A partir del valor de x en el instante t, se determina el valor de x en el instante t+h, y a partir del valor de y en el instante t se determina el valor de y en el instante t+h mediante las fórmulas de la última fila de la tabla.

dx dt =f(x,y,t) dy dt =g(x,y,t)
k 1 =h·f(x,y,t) k 2 =h·f( x+ 1 2 k 1 ,y+ 1 2 l 1 ,t+ 1 2 h ) k 3 =h·f( x+ 1 2 k 2 ,y+ 1 2 l 2 ,t+ 1 2 h ) k 4 =h·f( x+ k 3 ,y+ l 3 ,t+h ) l 1 =h·g(x,y,t) l 2 =h·g( x+ 1 2 k 1 ,y+ 1 2 l 1 ,t+ 1 2 h ) l 3 =h·g( x+ 1 2 k 2 ,y+ 1 2 l 2 ,t+ 1 2 h ) l 4 =h·g( x+ k 3 ,y+ l 3 ,t+h )
x(t+h)=x(t)+ 1 6 ( k 1 +2 k 2 +2 k 3 + k 4 ) y(t+h)=y(t)+ 1 6 ( l 1 +2 l 2 +2 l 3 + l 4 )

La jerarquía de clases

De modo similar, definimos la clase base abstracta denominada RungeKutta cuya función miembro resolver nos calcule la solución aproximada del sistema de ecuaciones diferenciales de primer orden (x, y) en el instante final tf, cuando le pasamos el estado inicial, es decir, el valor x0 e y0 en el instante inicial t0

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;
//estado inicial en el instante t
        double x=e.x;
        double y=e.y;
        double t=e.t;

           k1=h*f(x, y, t);
            l1=h*g(x, y, t);
            k2=h*f(x+k1/2, y+l1/2, t+h/2);
            l2=h*g(x+k1/2, y+l1/2, t+h/2);
            k3=h*f(x+k2/2, y+l2/2, t+h/2);
            l3=h*g(x+k2/2, y+l2/2, t+h/2);
            k4=h*f(x+k3, y+l3, t+h);
            l4=h*g(x+k3, y+l3, t+h);
            x+=(k1+2*k2+2*k3+k4)/6;
            y+=(l1+2*l2+2*l3+l4)/6;
            t+=h;
//estado final en el instante t+h
        e.t=t;
        e.x=x;
        e.y=y;
    }
    abstract public double f(double x, double y, double t);
    abstract public double g(double x, double y, double t);
}

La función resolver recibe en su segundo parámetro el estado inicial e del sistema. Ya que los objetos que se pasan a una función se pueden modificar en el curso de su llamada. La función resolver modifica dicho objeto, devolviendo el estado final cuando se le pasa el estado inicial.

La definición de la clase Estado cuyos objetos guardan el estado x, y del sistema en el instante t

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

Para obtener el estado e (objeto de la clase Estado) accedemos a sus miembros dato públicos t, x, y

	System.out.println("t " +e.t);
	System.out.println("x " +e.x);
	System.out.println("y " +e.y);

Para crear el estado e a partir de los valores de t, x, e y, llamamos al constructor de la clase  Estado.

        double t=0;
        double x=1000;
        double y=0;
	Estado e=new Estado(t, x, y);

Ejemplo

Continuando con el ejemplo de la desintegración radioactiva consideremos ahora, una serie radioactiva de tres elementos en la que, una sustancia radiactiva A se desintegra y se transforma en otra sustancia radiactiva B, que a su vez se desintegra y se transforma en una sustancia C estable. Las ecuaciones diferenciales que gobiernan el proceso y sus soluciones analíticas son, respectivamente,

dx dt =axx= x 0 exp(at) dy dt =axbyy= a ba x 0 ( exp(at)exp(bt) )

La solución analítica que aparece a la derecha, se ha obtenido con las condiciones iniciales t=0, x=x0 e y=0. La segunda solución se obtiene siempre que a sea distinto de b. En el caso de que a sea igual a b, la solución analítica para y es

y= x 0 aexp(at)

La interpretación del sistema de ecuaciones diferenciales no es complicada. En la unidad de tiempo, desaparecen ax núcleos de la sustancia A al desintegrarse (primera ecuación). En la unidad de tiempo, se producen ax núcleos de la sustancia B, y a su vez desaparecen bx núcleos de la sustancia B que al desintegrarse se transforman en núcleos de la sustancia C estable (segunda ecuación).

Para codificar este ejemplo, primero se han de definir las funciones f, g en la clase derivada Funcion de la clase abstracta RungeKutta.

public class Sistema extends RungeKutta{
    public Sistema(double h){
      super(h);
    }
    public double f(double x, double y, double t){
        return (-0.1*x);
    }
    public double g(double x, double y, double t){
        return (0.1*x-0.2*y);
    }
}

Definimos el estado inicial en el instante t0=0, el número de núcleos de la sustancia radiactiva A, suponemos que no hay inicialmente núcleos de sustancia B. Establecemos el instante t en el que deseamos conocer el número de núcleos que hay de la sustancia A, y el número de núcleos de la sustancia B, llamando a la función resolver desde un objeto de la clase Sistema.

public class Aplicacion {
  public static void main(String[] args) {
      double h=0.5;       //paso
      double x0=5000;
      Estado e=new Estado(0.0, 5000.0, 0);
      double t=20.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("valor aprox. de x "+(int)e.x);
    System.out.println("valor aprox. de y "+(int)e.y);

  }
}

Como en el apartado anterior fijamos el número de núcleos de la sustancia A en 5000, y el de la B ningún núcleo. El valor del paso h lo tomamos como 0.5. Los resultados aproximados para el instante t=20 están en concordancia con los obtenidos a partir del resultado exacto, como queda reflejado en la tabla adjunta

t=20 Aproximado Exacto
A 676 676
B 585 585

 

 

 

Siguiente Anterior