Demostración de la conservación del momento lineal

Procedimientos numéricos

Suistema de ecuaciones diferenciales de segundo orden

Resolver el sistema de ecuaciones diferenciales de segundo orden

$\begin{array}{l} m \frac{d^{2} x_{1}}{d t^{2}} = - k {(x_{1} - x_{2})}^{r} \\ m \frac{d^{2} x_{2}}{d t^{2}} = k {(x_{1} - x_{2})}^{r} - k {(x_{2} - x_{3})}^{r} \\ m \frac{d^{2} x_{i}}{d t^{2}} = k {(x_{i - 1} - x_{i})}^{r} - k {(x_{i} - x_{i + 1})}^{r} \\ m \frac{d^{2} x_{n}}{d t^{2}} = k {(x_{n - 1} - x_{n})}^{r} \end{array}$

Tenemos un sistema de n ecuaciones diferenciales acopladas que se resuelve por procedimientos numéricos con las siguientes condiciones iniciales

En el instante t=0,

Todas las partículas están en sus posiciones de equilibrio. x_i=0, i=1...n
Todas las partículas están en reposo, excepto la partícula incidente que lleva velocidad v en el instante t=0 en el que choca con la segunda partícula. v₁=v, v_i=0, i=2…n

La clase base que describe el procedimiento numérico

public abstract class RungeKutta {
  double h;
  int n;
  public RungeKutta(int n, double h) {
     this.h=h;
     this.n=n;
  }

  public void setPaso(double h){
        this.h=h;
  }
  public void resolver(Estado e){
     double[] k1=new double[n+1];
     double[] k2=new double[n+1];
     double[] k3=new double[n+1];
     double[] k4=new double[n+1];
     double[] l1=new double[n+1];
     double[] l2=new double[n+1];
     double[] l3=new double[n+1];
     double[] l4=new double[n+1];
//condiciones iniciales
      double[] x=e.x;
	    double[] v=e.v;
	    double t=e.t;
          
//variables auxiliares
     k1[0]=h*v[0];
     l1[0]=h*f(x[0], x[0], x[1]);
     k2[0]=h*(v[0]+l1[0]/2);
     l2[0]=h*f(x[0]+k1[0]/2, x[0]+k1[0]/2, x[1]+k1[1]/2);
     k3[0]=h*(v[0]+l2[0]/2);
     l3[0]=h*f(x[0]+k2[0]/2, x[0]+k2[0]/2, x[1]+k2[1]/2);
     k4[0]=h*(v[0]+l3[0]);
     l4[0]=h*f(x[0]+k3[0], x[0]+k3[0], x[1]+k3[1]);
     x[0]+=(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;
     v[0]+=(l1[0]+2*l2[0]+2*l3[0]+l4[0])/6;

     for(int i=1; i<n; i++){
     k1[i]=h*v[i];
     l1[i]=h*f(x[i-1], x[i], x[i+1]);
     k2[i]=h*(v[i]+l1[i]/2);
     l2[i]=h*f(x[i-1]+k1[i-1]/2, x[i]+k1[i]/2, x[i+1]+k1[i+1]/2);
     k3[i]=h*(v[i]+l2[i]/2);
     l3[i]=h*f(x[i-1]+k2[i-1]/2, x[i]+k2[i]/2, x[i+1]+k2[i+1]/2);
     k4[i]=h*(v[i]+l3[i]);
     l4[i]=h*f(x[i-1]+k3[i-1], x[i]+k3[i], x[i+1]+k3[i+1]);
     x[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
     v[i]+=(l1[i]+2*l2[i]+2*l3[i]+l4[i])/6;
  }
     k1[n]=h*v[n];
     l1[n]=h*f(x[n-1], x[n], x[n]);
     k2[n]=h*(v[n]+l1[n]/2);
     l2[n]=h*f(x[n-1]+k1[n-1]/2, x[n]+k1[n]/2, x[n]+k1[n]/2);
     k3[n]=h*(v[n]+l2[n]/2);
     l3[n]=h*f(x[n-1]+k2[n-1]/2, x[n]+k2[n]/2, x[n]+k2[n]/2);
     k4[n]=h*(v[n]+l3[n]);
     l4[n]=h*f(x[n-1]+k3[n-1], x[n]+k3[n], x[n]+k3[n]);
     x[n]+=(k1[n]+2*k2[n]+2*k3[n]+k4[n])/6;
     v[n]+=(l1[n]+2*l2[n]+2*l3[n]+l4[n])/6;

     t+=h;
//cambia el estado de la partícula
     for(int i=1; i<n+1; i++){
	    e.x[i]=x[i];
	    e.v[i]=v[i];
    }
	    e.t=t;
  }
  public abstract double f(double x, double y, double z);
}

La clase derivada que describe el sistema físico particular

public class Sistema extends RungeKutta{
    double exponente;
    double cte;
    public Sistema(int n, double exponente, double cte, double h){
        super(n, h);
        this.exponente=exponente;
        this.cte=cte;
    }
    public double f(double x0, double x, double x1){
        if((x0-x)<0)  x0=x;
        if((x-x1)<0)  x1=x;
        double temp=cte*Math.pow((x0-x), exponente)-
cte*Math.pow((x-x1), exponente);
        return temp;
    }
}

La clase cuyos objetos guardan el estado del sistema

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

Se establece el estado incial

    final double cte=1.638e10/0.512;
    double[] x=new double[n];
    double[] v=new double[n];
    for(int i=0; i<n; i++){
        x[i]=0.0;
        v[i]=0.0;
    }
    v[0]=Math.sqrt(2*9.8*0.01);
    Estado estado=new Estado(0.0, x, v);

Se crea un objeto de la clase derivada

 Sistema sis=new Sistema(n-1, 1.5, cte, 1.0e-6);

Se llama a la función resolver que determina el estado del sistema en el instante t+h conocido el estado en el instante t

sis.resolver(estado);