Anterior

Procedimiento numérico

Resolver el sistema de ecuaciones diferenciales

d 2 V 1 d t 2 = S 2 m ( k 1 ( V 1 ) γ 1 k 2 ( V T V 1 ) γ 2 ) e m d V 1 d t d k 1 d t = α e R c v 1 S 2 ( d V 1 d t ) 2 ( V 1 ) γ 1 1 d k 2 d t = ( 1 α ) e R c v 2 S 2 ( d V 1 d t ) 2 ( V T V 1 ) γ 2 1

por procedimientos numéricos, con las siguientes condiciones iniciales, en el instante t=0.

V 1 = V 10 ( d V 1 d t ) = 0 k 1 = p 10 V 10 γ 1 k 2 = p 20 V 20 γ 2

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 class Gas {
  public double presion;
  public double volumen;
  public double nMoles;
  public double  cV;
  public double temperatura;
  double gamma;
  public Gas(double presion, double volumen, double nMoles, boolean bMono) {
    this.presion=presion;
    this.volumen=volumen;
    this.nMoles=nMoles;
    temperatura=presion*volumen*100/(nMoles*8.3143);
    this.cV=(bMono)?(3*8.3143/2):(5*8.3143/2);
    this.gamma=(bMono)?(5.0/3):(7.0/5);
  }
}

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{
    Gas[] gas;
    double roza;
    double masa;
    double kA, kB;
    final double fraccion=0.5;
    public Sistema(Gas[] gas, double roza, double masa, double h){
      super(h);
      this.gas=gas;
      this.roza=roza;
      this.masa=masa;
      kA=8.3143*fraccion*roza/(4*gas[0].cV);
      kB=8.3143*(1.0-fraccion)*roza/(4*gas[1].cV);
    }
    //volumen en litros, presión en 1.0e5 Pa
    public double f(double x, double v, double t, double y, double z){
         double temp=4*(y/Math.pow(x, gas[0].gamma)-z/
Math.pow((gas[0].volumen+gas[1].volumen-x), gas[1].gamma))-roza*v; return (temp/masa); } public double g(double x, double v, double t, double y, double z){ return (kA*v*v*Math.pow(x, (gas[0].gamma-1))); } public double c(double x, double v, double t, double y, double z){ return (kB*v*v*Math.pow((gas[0].volumen+gas[1].volumen-x), (gas[1].gamma-1))); } } public class MiCanvas extends Canvas{ //... Gas[] gas=new Gas[2]; //objeto sistema Sistema sistema; Estado estado; public MiCanvas(ProcesoApplet7 p) { gas[0]=new Gas(9.0, 1.0, 0.08, true); gas[1]=new Gas(2.0, 2.0, 0.1, false); estado=new Estado(0.0, 1.0, 0.0, 0.0, 0.0); } void setNuevo(Gas[] gas, int nGrafica){ this.gas=gas; double cte_0=gas[0].presion*Math.pow(gas[0].volumen, gas[0].gamma); double cte_1=gas[1].presion*Math.pow(gas[1].volumen, gas[1].gamma); estado=new Estado(0.0, gas[0].volumen, 0.0, cte_0, cte_1); } void mover(){ sistema.resolver(estado); gas[0].volumen=estado.x; gas[1].volumen=volTotal-estado.x; gas[0].presion=estado.y/Math.pow(gas[0].volumen, gas[0].gamma); gas[1].presion=estado.z/Math.pow(gas[1].volumen, gas[1].gamma); gas[0].temperatura=100*gas[0].presion*gas[0].volumen/(8.3143*gas[0].nMoles); gas[1].temperatura=100*gas[1].presion*gas[1].volumen/(8.3143*gas[1].nMoles); } //... }
Anterior