Movimiento de un disco que gira y desliza sobre una superficie plana y horizontal

Procedimientos numéricos

Integral definida por el procedimiento de Simpson
Integral doble por el procedimiento de Simpson
Procedimiento de Runge-Kutta para resolver un sistema de dos ecuaciones diferenciales de segundo orden.

Aro

$\begin{array}{l} \frac{d^{2} x}{d t^{2}} = - \frac{μ g}{2 π} \int_{0}^{2 π} \frac{v - ω R sin θ}{\sqrt{v^{2} + ω^{2} R^{2} - 2 v ω R sin θ}} d θ \\ \frac{d^{2} ϕ}{d t^{2}} = - \frac{μ g}{2 π R^{2}} \int_{0}^{2 π} (\frac{R (ω R - v sin θ)}{\sqrt{v^{2} + ω^{2} R^{2} - 2 v ω R sin θ}}) d θ \end{array}$

Disco

$\begin{array}{l} \frac{d^{2} x}{d t^{2}} = - \frac{μ g}{π R^{2}} \int_{0}^{2 π} d θ \int_{0}^{R} \frac{v - ω r sin θ}{\sqrt{v^{2} + ω^{2} r^{2} - 2 v ω r sin θ}} r d r \\ \frac{d^{2} ϕ}{d t^{2}} = - \frac{2 μ g}{π R^{4}} \int_{0}^{2 π} d θ \int_{0}^{R} (\frac{r^{2} (ω r - v sin θ)}{\sqrt{v^{2} + ω^{2} r^{2} - 2 v ω r sin θ}}) d r \end{array}$

public abstract class Simpson {
	int n;
	final double a=0.0;
	final double b=2*Math.PI;

public Simpson(int n){
	if(n%2==1) n++; //hace que el número de divisiones sea par
	this.n=n;
}
public double integral(){
	double h=(b-a)/n;
	double suma=f(a)+f(b);
	for(int i=1; i<n; i+=2){
		suma+=4*f(a+i*h);
	}
	for(int i=2; i<n; i+=2){
		suma+=2*f(a+i*h);
	}
	return (suma*h/3);
}
abstract public double f(double x);
}

public class Aro_Vcm extends Simpson{
	double v, w;
public Aro_Vcm(int n, double v, double w){
	super(n);
	this.v=v;
	this.w=w;
}
void setVelocidades(double v, double w){
	this.v=v;
	this.w=w;
}
public double f(double ang){ //radio R=unidad
	double num=(v-w*Math.sin(ang));
	double den=Math.sqrt(v*v+w*w-2*v*w*Math.sin(ang));
	return (num/den);
}
}

public class Aro_Vangular extends Simpson{
	double v, w;
public Aro_Vangular(int n, double v, double w){
	super(n);
	this.v=v;
	this.w=w;
}
void setVelocidades(double v, double w){
	this.v=v;
	this.w=w;
}
public double f(double ang){ //radio R=unidad
	double num=(w-v*Math.sin(ang));
	double den=Math.sqrt(v*v+w*w-2*v*w*Math.sin(ang));
	return (num/den);
}
}

public class Aro extends RungeKutta2{
	double mu;
	Aro_Vcm cm=null;
	Aro_Vangular angular=null;
Aro(double mu, double v0, double w0, double h){
	super(h);
	this.mu=mu;
	cm=new Aro_Vcm(90, v0, w0);
	angular=new Aro_Vangular(90, v0, w0);
}
public double f(double x, double y, double v, double w, double t){
	cm.setVelocidades(v, w);
	double temp=-mu*9.8*cm.integral()/(2*Math.PI);
	return temp;
}
public double g(double x, double y, double v, double w, double t){
	angular.setVelocidades(v, w);
	double temp=-mu*9.8*angular.integral()/(2*Math.PI);
	return temp;
}
}

public abstract class SimpsonDoble {        
	public double integral(double a, double b, int n, int m){
        double h=(b-a)/(2*n);
        double J1=0.0, J2=0.0, J3=0.0;
        double x, y;
        double HX, K1, K2, K3, Q, L;
        for(int i=0; i<=2*n; i++){
            x=a+i*h;
            HX=(d(x)-c(x))/(2*m);
            K1=f(x, c(x))+f(x, d(x));
            K2=0.0;
            K3=0.0;
            for(int j=1; j<=2*m-1; j++){
                y=c(x)+j*HX;
                Q=f(x, y);
                if(j%2==0) K2+=Q;
                else K3+=Q;
            }
            L=(K1+2*K2+4*K3)*HX/3;
            if(i==0 || i==2*n)  J1+=L;
            else{
              if(i%2==0) J2+=L;
              else J3+=L;
           }
        }
        double J=h*(J1+2*J2+4*J3)/3;
        return J;
    }
    abstract public double f(double x, double y);
    abstract public double c(double x);
    abstract public double d(double x);
}

public class Disco_Vcm extends SimpsonDoble{
	double v, w;
public Disco_Vcm(double v, double w){ //d/a separación dividido radio
		this.v=v;
	this.w=w;
}
void setVelocidades(double v, double w){
	this.v=v;
	this.w=w;
}
public double f(double ang, double r){
	double num=r*(v-w*r*Math.sin(ang));
	double den=Math.sqrt(v*v+w*r*w*r-2*v*w*r*Math.sin(ang));
	return (num/den);
}
    public double c(double x){
      return(0.0);
  }
  public double d(double x){
      return(1.0);       //radio
  }
}

public class Disco_Vangular extends SimpsonDoble{
	double v, w;
 public Disco_Vangular(idouble v, double w){ 
	this.v=v;
	this.w=w;
 }
 void setVelocidades(double v, double w){
	this.v=v;
	this.w=w;
 }
 public double f(double ang, double r){
	double num=r*r*(w*r-v*Math.sin(ang));
	double den=Math.sqrt(v*v+w*r*w*r-2*v*w*r*Math.sin(ang));
	return (num/den);
 }
 public double c(double x){
      return(0.0);
  }
  public double d(double x){
      return(1.0);       //radio
  }
}

public abstract class RungeKutta2 {
	double h=0.01;
 public RungeKutta2(double h){
	this.h=h;
}
public void setPaso(double h){
	this.h=h;
}
public void resolver(Estado e){
//variables auxiliares
	double k1, k2, k3, k4;
	double l1, l2, l3, l4;
	double q1, q2, q3, q4;
	double m1, m2, m3, m4;
//condiciones iniciales
	double x=e.x;
	double y=e.y;
	double vx=e.v;
	double vy=e.w;
	double t=e.t;

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

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

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

	k4=h*(vx+l3);
	l4=h*f(x+k3, y+q3, vx+l3, vy+m3, t+h);
	q4=h*(vy+m3);
	m4=h*g(x+k3, y+q3, vx+l3, vy+m3, t+h);
//nuevo estado del sistema
	x+=(k1+2*k2+2*k3+k4)/6;
	vx+=(l1+2*l2+2*l3+l4)/6;
	y+=(q1+2*q2+2*q3+q4)/6;
	vy+=(m1+2*m2+2*m3+m4)/6;

//cambia el estado de la partícula
	e.x=x;
	e.y=y;
	e.v=vx;
	e.w=vy;
	e.t=t+h;
}
abstract public double f(double x, double y, double vx, double vy, double t);
abstract public double g(double x, double y, double vx, double vy, double t);
}

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

public class Disco extends RungeKutta2{
    double mu;
    Disco_Vcm cm=null;
    Disco_Vangular angular=null;
    Disco(double mu, double v0, double w0, double h){
        super(h);
        this.mu=mu;
        cm=new Disco_Vcm(v0, w0);
        angular=new Disco_Vangular(v0, w0);
    }
    public double f(double x, double y, double v, double w, double t){
        cm.setVelocidades(v, w);
        double temp=-mu*9.8*cm.integral(0.0, 2*Math.PI, 20, 20)/Math.PI;
        return temp;
    }
    public double g(double x, double y, double v, double w, double t){
        angular.setVelocidades(v, w);
        double temp=-2*mu*9.8*angular.integral(0.0, 2*Math.PI, 20, 20)/Math.PI;
        return temp;
    }
}

//creación de objetos
//v0 y w0 son las velocidades iniciales

Estado estado=new Estado(0.0, 0.0, 0.0, v0, w0);
RungeKutta2 objeto=new Disco(mu, v0, w0, dt);    //polimorfismo
if(bAro){
	objeto=new Aro(mu, v0, w0, dt);
}

//calcula la posición y la velocidad del disco o del aro en función del tiempo
objeto.resolver(estado);