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