## Sistema de ecuaciones diferenciales de segundo orden

$\begin{array}{l}\frac{{d}^{2}{x}_{R}}{d{t}^{2}}-2·0.231\frac{d{y}_{R}}{dt}-{0.231}^{2}{x}_{R}=-11519.568\frac{\left({x}_{R}+0.731\right)}{{d}_{T}^{3}}-141.394\frac{\left({x}_{R}-59.552\right)}{{d}_{L}^{3}}\\ \frac{{d}^{2}{y}_{R}}{d{t}^{2}}+2·0.231\frac{d{x}_{R}}{dt}-{0.231}^{2}{y}_{R}=-11519.568\frac{{y}_{R}}{{d}_{T}^{3}}-141.394\frac{{y}_{R}}{{d}_{L}^{3}}\\ {d}_{T}=\sqrt{{\left({x}_{R}+0.731\right)}^{2}+{y}_{R}^{2}}\text{ }\text{ }\text{ }{d}_{L}=\sqrt{{\left({x}_{R}-59.552\right)}^{2}+{y}_{R}^{2}}\end{array}$

La constante J se escribe

$J=\frac{1}{2}\left({\stackrel{˙}{x}}_{R}^{2}+{\stackrel{˙}{y}}_{R}^{2}\right)-\frac{1}{2}{0.231}^{2}\left({x}_{R}^{2}+{y}_{R}^{2}\right)-\left(\frac{11519.568}{{d}_{T}}+\frac{141.394}{{d}_{L}}\right)$

En el instante t=0, la nave espacial parte de la posición

x0=-rT+r·cosθ
y0
= r·sinθ

 public class Sistema extends RungeKutta{ final double rT=0.73095; final double rL=59.55162; final double w=0.23072; // wTierra-Luna*un día final double cL=141.39403; final double cT=11519.56835; pulic Sistema(double h){ super(h); } public double f(double x, double y, double vx, double vy, double t){ double dL=Math.sqrt((x-rL)*(x-rL)+y*y); double dT=Math.sqrt((x+rT)*(x+rT)+y*y); double z=2*w*vy+w*w*x-cL*(x-rL)/(dL*dL*dL)-cT*(x+rT)/(dT*dT*dT); return z; } public double g(double x, double y, double vx, double vy, double t){ double dL=Math.sqrt((x-rL)*(x-rL)+y*y); double dT=Math.sqrt((x+rT)*(x+rT)+y*y); double z=-2*w*vx+w*w*y-cL*y/(dL*dL*dL)-cT*y/(dT*dT*dT); return z; } public double energia(double x, double y, double vx, double vy){ double dL=Math.sqrt((x-rL)*(x-rL)+y*y); double dT=Math.sqrt((x+rT)*(x+rT)+y*y); double z=(vx*vx+vy*vy)/2-(cL/dL+cT/dT)-w*w*(x*x+y*y)/2; return z; } } //Objeto de la clase Sistema void setNuevo(double v0, double ang){ double x0=-rT+rNave*Math.cos(ang); double y0=rNave*Math.sin(ang); double v0x=-(rNave*wNave+v0*3600*24/6.37e6)*Math.sin(ang); double v0y=(rNave*wNave+v0*3600*24/6.37e6)*Math.cos(ang); estado=new Estado(0.0, x0, y0, v0x, v0y); sis=new Sistema(calculaPaso()); eInicial=oscilador.energia(x0, y0, v0x, v0y); } double calculaPaso(){ double dL=Math.sqrt((estado.x-rL)*(estado.x-rL)+estado.y*estado.y); double dT=Math.sqrt((estado.x+rT)*(estado.x+rT)+estado.y*estado.y); double distMed=Math.sqrt(dT*dL); double dt=0.0005*Math.pow((distMed/(60.282*0.1)), 1.5); //para un día return dt; } void mover(){ sis.resolver(estado); sis.setPaso(calculaPaso()); double energia=oscilador.energia(estado.x, estado.y, estado.vx, estado.vy); error=Math.abs((energia-eInicial)*100/eInicial); if(error>1.0){ //se para ; } } }