Sistema de ecuaciones diferenciales de primer orden

prev.gif (997 bytes)chapter.gif (1105 bytes)home.gif (1232 bytes)next.gif (1211 bytes)

Ecuaciones diferenciales

Sistema de ecuaciones diferenciales de primer orden

Código fuente


Sistema de ecuaciones diferenciales de primer orden

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

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

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.

 

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 {
    public void resolver(double tf, Estado e, double h){
//variables auxiliares
        double k1, k2, k3, k4;
        double l1, l2, l3, l4;
//estado inicial
        double x=e.x;
        double y=e.y;
        double t0=e.t;

        for(double t=t0; t<tf; t+=h){
            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;
        }
//estado final
        e.t=tf;
        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);
}

Una función solamente, puede devolver un valor, sin embargo, el estado del sistema en cualquier instante t viene dado por dos números x e y. Para ello, definimos la clase denominada Estado con tres miembros públicos x e y y el tiempo t.

La función resolver recibe en su segundo parámetro el estado inicial 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.

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,

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

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 Funcion extends RungeKutta{
    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 Funcion.

  double h=0.5;       //paso

  double x0=5000;
  Estado es=new Estado(0, x0, 0);

  double t=20.0;      //resolver la e. d. hasta este instante
  new Funcion().resolver(t, es, h);
  System.out.println("valor aprox. de x "+(int)es.x);
  System.out.println("valor aprox. de y "+(int)es.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

 

 

 

 

El código fuente

disco.gif (1035 bytes) RungeKutta.java, Funcion.java, Estado.java, RkApp2java