Anterior

Procedimiento numérico

Resolver el sistema de ecuaciones diferenciales

d θ 1 d τ = − ( γ − 1 ) θ 1 ξ d ξ d τ + δ ( γ − 1 ) θ 1 ξ ( d ξ d τ ) 2 − π 8 ( γ − 1 ) δ 2 1 ξ ( d ξ d τ ) 3 d θ 2 d τ = ( γ − 1 ) θ 2 1 − ξ d ξ d τ + δ ( γ − 1 ) θ 2 1 − ξ ( d ξ d τ ) 2 + π 8 ( γ − 1 ) δ 2 1 1 − ξ ( d ξ d τ ) 3 d 2 ξ d τ 2 = θ 1 ξ − θ 2 1 − ξ − δ ( θ 1 ξ + θ 2 1 − ξ ) d ξ d τ + π 8 δ 2 ( 1 ξ − 1 1 − ξ ) ( d ξ d τ ) 2

con las siguientes condiciones iniciales, en el instante τ=0, ξ=ξ0, (dξ/dτ)0=0,θ1=θ10, θ2=1-θ10

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

public abstract class RungeKutta {
    double h;
    RungeKutta(double h){
      this.h=h;
    }
    public void resolver(Estado e){
//variables auxiliares
	    double k1, k2, k3, k4;
	    double l1, l2, l3, l4;
	    double m1, m2, m3, m4;
	    double q1, q2, q3, q4;

//condiciones iniciales
	    double x=e.x;
	    double v=e.v;
	    double t=e.t;
	    double y=e.y;
	    double z=e.z;

	    //for(double t=t0; t<tf; t+=h){
		    k1=h*v;
		    l1=h*f(x, v, t, y, z);
        m1=h*g(x, v, t, y, z);
        q1=h*c(x, v, t, y, z);

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

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

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

//nuevo estado del sistema
		    x+=(k1+2*k2+2*k3+k4)/6;
		    v+=(l1+2*l2+2*l3+l4)/6;
        y+=(m1+2*m2+2*m3+m4)/6;
        z+=(q1+2*q2+2*q3+q4)/6;

	   // }
//cambia el estado 
	    e.x=x;
	    e.v=v;
	    e.y=y;
	    e.z=z;
	    e.t=t+h;
    }
    abstract public double f(double x, double v, double t, double y, double z);
    abstract public double g(double x, double v, double t, double y, double z);
    abstract public double c(double x, double v, double t, double y, double z);

}
public class Sistema extends RungeKutta{
    final double gamma=5.0/3;
    double delta;
    Sistema(double delta, double h){
      super(h);
      this.delta=delta;
    }
    public double f(double x, double v, double t, double y, double z){
         double temp=y/x-z/(1-x)-delta*v*(Math.sqrt(y)/x+
Math.sqrt(z)/(1-x))+(Math.PI/8)*delta*delta*v*v*(y/x-z/(1-x));
         return temp;
    }
    public double g(double x, double v, double t, double y, double z){
         double temp=-(gamma-1)*y*v/x+delta*(gamma-1)*v*v*Math.sqrt(y)/x
-(Math.PI/8)*(gamma-1)*delta*delta*v*v*v/x;
         return temp;
    }
    public double c(double x, double v, double t, double y, double z){
         double temp=(gamma-1)*z*v/(1-x)+delta*(gamma-1)*v*v*Math.sqrt(z)/(1-x)
+(Math.PI/8)*(gamma-1)*delta*delta*v*v*v/(1-x);
         return temp;
    }

}
//objetos de la clase Sistema
    estado=new Estado(0.0, x0, 0.0, temp1, 1.0-temp1);
    sistema=new Sistema(delta,  0.02);
Anterior