Caída de un lápiz con su punta apoyada en un suelo rugoso

Procedimiento numérico

Sistema de ecuaciones diferenciales de segundo orden

Cuando no desliza

$\frac{d^{2} θ}{d t^{2}} = \frac{3 g}{2 L} sin θ$

Para obtener el ángulo θ que hace la varilla con el suelo en función del tiempo, se integra la ecuación diferencial mediante procedimientos numéricos con las condiciones iniciales siguientes: en el instante t=0, la varilla está inclinada un ángulo θ=θ₀ y parte del reposo, ω=dθ/dt=0

Cuando desliza

$\begin{array}{l} \frac{d^{2} θ}{d t^{2}} = \frac{(6 g - 3 L {(\frac{d θ}{d t})}^{2} \cos θ) (sin θ - μ \cos θ)}{L (1 + 3 sin θ (sin θ - μ \cos θ))} \\ m \frac{d^{2} x}{d t^{2}} = F = \pm μ N N = m g - m \frac{L}{2} {\cos θ {(\frac{d θ}{d t})}^{2} + sin θ \frac{d^{2} θ}{d t^{2}}} \end{array}$

El estado inicial corresponde al estado final de la etapa anterior

public class Estado {
	double t;
	double x;
	double y;
	double vx;
	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 abstract class RungeKutta {
	double h;
public RungeKutta(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;
//estado inicial
	double x=e.x;
	double y=e.y;
	double vx=e.vx;
	double vy=e.vy;
	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);

	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;
	t+=h;

//estado final
	e.x=x;
	e.y=y;
	e.vx=vx;
	e.vy=vy;
	e.t=t;
}
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);

abstract public double fuerza_X(double x, double vx);
abstract public double fuerza_Y(double x, double vx);
}

public class Sistema extends RungeKutta{
	double lonVarilla;
public Sistema(double lonVarilla, double h){
	super(h);
	this.lonVarilla=lonVarilla;
}
public double f(double x, double y, double vx, double vy, double t){
	return (14.7*Math.sin(x)/lonVarilla);
}
public double g(double x, double y, double vx, double vy, double t){
	double temp=(f(x, 0.0, vx, 0.0, 0.0)*Math.cos(x)-vx*vx*Math.sin(x))
*lonVarilla/2;
	return temp;
}

public double fuerza_X(double x, double vx){
	double temp=(f(x, 0.0, vx, 0.0, 0.0)*Math.cos(x)-vx*vx*Math.sin(x)
)*lonVarilla/2;
	return temp;
}
public double fuerza_Y(double x, double vx){
	double temp=9.8-(f(x, 0.0, vx, 0.0, 0.0)*Math.sin(x)+vx*vx*Math.cos(x))*
lonVarilla/2;
	return temp;
} 
}

public class Sistema_Roza extends Sistema {
	double mu;
public Sistema_Roza(double lonVarilla, double mu, double h){
	super(lonVarilla, h);
	this.mu=mu;
}

public double f(double x, double y, double vx, double vy, double t){
	double temp=(6*9.8-3*lonVarilla*vx*vx*Math.cos(x))*(Math.sin(x)-mu*Math.cos(x))
	/(lonVarilla*(1.0+3*Math.sin(x)*(Math.sin(x)-mu*Math.cos(x))));
	return temp;
}
public double g(double x, double y, double vx, double vy, double t){
	double temp=9.8-lonVarilla*(vx*vx*Math.cos(x)+f(x, y, vx, vy, t)*Math.sin(x))/2;
	return (mu*temp);
}
public double fuerza_X(double x, double vx){
	double temp=mu*fuerza_Y(x, vx);
	return temp;
}

}

//Creación de objetos

	estado=new Estado(0.0, angInicial, 0.0, 0.0, 0.0);
	sistema=new Sistema(lonVarilla, dt);

//movimiento, resuelve las ecuaciones diferenciales

 	sistema.resolver(estado);
	double vPunta=estado.vy-estado.vx*lonVarilla*Math.cos(estado.x)/2;
	if((mu<muMax)&& (vPunta>0) && !bMueveDcha){ //velocidad de la punta del lápiz
		bMueveDcha=true;
		bMueveIzq=false;
		sistema=new Sistema(lonVarilla, dt);
	}
	if(Math.abs(sistema.fuerza_X(estado.x, estado.vx))>
	mu*sistema.fuerza_Y(estado.x, estado.vx)) &&!bMueveIzq){
		if(sistema.fuerza_X(estado.x, estado.vx)<0.0)
			sistema=new Sistema_Roza(lonVarilla, -mu, dt);
		else
			sistema=new Sistema_Roza(lonVarilla, mu, dt);
		bMueveIzq=true;
	}