Anterior

Procedimiento numérico

Resolver el sistema formado por una ecuación diferencial de segundo orden y una ecuación diferencial de primer orden por el procedimiento de Runge-Kutta

d 2 x d t 2 = − 1 4 π ε 0 q 2 4 m 1 x 2 + q 2 B 2 m 2 ( x 0 − x ) d y d t = q B m ( x 0 − x )

con las condiciones iniciales siguientes: en el instante t=0, x=x0, y=0, vx=0.

Con los datos: masa del protón m=1.6725·10-27 kg, carga q=1.6·10-19 C, campo magnético B en gauss (10-4 T), posición ( x, y) en mm (10-3 m) y tiempo t en ms (10-3 s) tenemos el siguiente sistema de ecuaciones diferenciales

d 2 x d t 2 =−34.43946 1 x 2 +91.51825· B 2 ( x 0 −x) dy dt =9.56652·B·( x 0 −x) 

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

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

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

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

		    k4=h*(vx+l3);
		    l4=h*f(x+k3, y+q3, vx+l3,  t+h);
		    q4=h*g(x+k3, y+q3, vx+l3,  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;
        t+=h;

//estado final
	    e.x=x;
	    e.y=y;
	    e.vx=vx;
	    e.t=t;
    }
    abstract public double f(double x, double y, double vx, double t);
    abstract public double g(double x, double y, double vx, double t);
}
public class Particula extends RungeKutta{
    double campo;
    double x0;
    public Particula(double campo, double x0,  double h){
      super(h);
      this.campo=campo;
      this.x0=x0;
    }
    public double f(double x, double y, double vx, double t){
         double z=-34.43946/(x*x)+91.51825114*campo*campo*(x0-x);
         return z;
    }
    public double g(double x, double y, double vx,double t){
        double z=9.56651719*campo*(x0-x);
         return z;
    }

}
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;
    }
}
public class MiCanvas extends Canvas {
//objeto particula
   Estado estado=new Estado(0.0, 2.0, 0.0,  0.0, 0.0);
   Particula particula=null;
//...
 void setNuevo(double campo, double x0){
    estado=new Estado(0.0, x0, 0, 0, 0);
    particula=new Particula(campo, x0, 0.005);
 }
 void mover(){
     particula.resolver(estado);
  if(estado.x<=0.0){
//se detiene
      estado.x=0.0;
   }
 }
Anterior