Anterior

Procedimiento numérico

Ecuación diferencial de segundo orden y sistema de dos ecuaciones diferenciales de segundo orden

( I c +m L 2 4 −m L 2 2 μ sinθcos⁡θ ) d 2 θ d t 2 =( m L 2 2 μ cosθcos⁡θ ) ( dθ dt ) 2 +mg L 2 sinθ−mg L 2 μcos⁡θ

Se resuelve la ecuación diferencial por procedimientos numéricos con las siguientes condiciones iniciales: En el instante t=0, dθ/dt=0, θ=θ0.

d 2 x d t 2 =g− L 2 { cos⁡θ ( dθ dt ) 2 +sinθ d 2 θ d t 2 } ( I c +m L 2 4  sinθ(sinθ−μcos⁡θ) ) d 2 θ d t 2 =( −m L 2 4 cosθ(sinθ−μcos⁡θ) ) ( dθ dt ) 2 +mg L 2 (sinθ−μcos⁡θ)

Se resuelve este sistema de dos ecuaciones diferenciales con las siguientes condiciones iniciales: en el instante t1, el ángulo que forma la escalera con la vertical es θ1, la velocidad angular de rotación es (dθ/dt)1 y la velocidad horizontal del centro de masas es

( dx dt ) 1 = L 2 cos⁡ θ 1 ( dθ dt ) 1 

public abstract class State {
	double t;
	double x;
	double vx;
public State(double t, double x, double vx) {
	this.t=t;
	this.x=x;
	this.vx=vx;
}
}
public class Estado extends State{
public Estado(double t, double x, double vx) {
	super(t, x, vx);
}
}
public class Estado1 extends State{
	double y;
	double vy;
public Estado1(double t, double x, double y, double vx, double vy) {
	super(t, x, vx);
	this.y=y;
	this.vy=vy;
}
}
public interface RK {
abstract public void resolver(State e);
}
      
public abstract class RungeKutta implements RK{
	double h;
RungeKutta(double h){
	this.h=h;
}
public void resolver(State 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 v=e.vx;
	double t=e.t;

	k1=h*v;
	l1=h*f(x, v, t);
	k2=h*(v+l1/2);
	l2=h*f(x+k1/2, v+l1/2, t+h/2);
	k3=h*(v+l2/2);
	l3=h*f(x+k2/2, v+l2/2, t+h/2);
	k4=h*(v+l3);
	l4=h*f(x+k3, v+l3, t+h);
//nuevo estado del sistema
	x+=(k1+2*k2+2*k3+k4)/6;
	v+=(l1+2*l2+2*l3+l4)/6;
//cambia el estado de la partícula
	e.x=x;
	e.vx=v;
	e.t=t+h;
}
abstract public double f(double x, double v, double t);
}
public abstract class RungeKutta1 implements RK{
	double h;
RungeKutta1(double h){
	this.h=h;
}
public void resolver(State 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=((Estado1)e).y;
	double vx=e.vx;
	double vy=((Estado1)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;
	((Estado1)e).y=y;
	e.vx=vx;
	((Estado1)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);
}
      
public class Sistema extends RungeKutta{
	double mInercia;
	final double m=10.0;
	final double lonVarilla=1.0;
	double mu;
Sistema(double mu, double h){
	super(h);
	this.mu=mu;
	this.mInercia=m*lonVarilla*lonVarilla/12;
}
public double f(double x, double v, double t){ 
// x es ángulo, v es velocidad angular
	double seno=Math.sin(x);
	double coseno=Math.cos(x);
	double A=mInercia+m*lonVarilla*lonVarilla/4
-m*lonVarilla*lonVarilla*mu*seno*coseno/2;
	double B=m*lonVarilla*lonVarilla*mu*coseno*coseno/2;
	double C=m*9.8*lonVarilla*seno/2-m*9.8*mu*lonVarilla*coseno;
	return ((B*v*v+C)/A);
}
public double fuerza_X(double x, double v){
	double Fx=mu*9.8+(-mu*Math.sin(x)+Math.cos(x))*f(x, v, 0.0)*lonVarilla/2-
		(mu*Math.cos(x)+Math.sin(x))*v*v*lonVarilla/2;
	return (m*Fx);
}
public double fuerza_Y(double x, double vx){
	double temp=9.8-(f(x, vx, 0.0)*Math.sin(x)+vx*vx*Math.cos(x))*lonVarilla/2;
	return (m*temp);
}
}
      
public class Sistema1 extends RungeKutta1{
double mInercia;
	final double m=10.0;
	final double lonVarilla=1.0;
	double mu;
Sistema1(double mu, double h){
	super(h);
	this.mu=mu;
	this.mInercia=m*lonVarilla*lonVarilla/12;
}
public double g(double x, double y, double vx, double vy, double t){
	return (-mu*fuerza_Y(x, vx)/m);
}
public double f(double x, double y, double vx, double vy, double t){
	double seno=Math.sin(x);
	double coseno=Math.cos(x);
	double A=mInercia+m*lonVarilla*lonVarilla*seno*(seno-mu*coseno)/4;
	double B=-m*lonVarilla*lonVarilla*coseno*(seno-mu*coseno)/4;
	double C=m*9.8*lonVarilla*(seno-mu*coseno)/2;
	return ((B*vx*vx+C)/A);
}

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 (m*temp);
}
}

public class MiCanvas extends Canvas {
 //objeto sistema
	RK sistema;
	double angInicial=Math.PI/6;
	State estado=new Estado(0.0, angInicial, 0.0);
//otros miembros dato...
//funciones miembro
 void setNuevo(int ang, double muEst, double muDin){
	this.angInicial=ang*Math.PI/180;
	this.mu=muDin;
	tipo=1;
	estado=new Estado(0.0, angInicial, 0.0);
	if(angInicial>=Math.atan(2*muEst)){
		sistema=new Sistema(mu, dt);
		tipo=1;
	}else{
	tipo=3;
	}
	xExtremo=lonVarilla*Math.sin(angInicial);
	repaint();
}

void mover(){
	switch(tipo){
		case 1:
			sistema.resolver(estado);
		   xExtremo=lonVarilla*Math.sin(estado.x);
			double fX=((Sistema)sistema).fuerza_X(estado.x, estado.vx);
			if(fX<0){
				double fY=((Sistema)sistema).fuerza_Y(estado.x, estado.vx);
				sistema=new Sistema1(mu, dt);
				double x0=lonVarilla*Math.sin(estado.x)/2; 
				double vx0=lonVarilla*Math.cos(estado.x)*estado.vx/2;
				estado=new Estado1(estado.t, estado.x, x0, estado.vx, vx0);
				fY=((Sistema1)sistema).fuerza_Y(estado.x, estado.vx);
				tipo=2;
			}
			break;
		case 2:
			sistema.resolver((Estado1)estado);
			xExtremo=((Estado1)estado).y+lonVarilla*Math.sin(estado.x)/2;
			double vPunta=((Estado1)estado).vy+estado.vx*lonVarill
a*Math.cos(estado.x)/2;
			double dist=((Estado1)estado).y-lonVarilla*Math.sin(estado.x)/2;
			if(dist<0){
				parent.hilo.putMsg(Hilo.PAUSE);
				System.out.println("vuelve a tocar");
			}
			if(estado.x>Math.PI/2){
				estado.x=Math.PI/2;
				parent.hilo.putMsg(Hilo.PAUSE);
			}
			if(vPunta<0){
				parent.hilo.putMsg(Hilo.PAUSE);
				System.out.println("se para");
			}
			break;
		case 3:
			break;
		default:
		break;
	}
	repaint();
}
Anterior