Siguiente Anterior

Integrales elípticas

Algunas páginas del Curso Interactivo de Física describen situaciones físicas en las que es necesario calcular una integral elíptica de primera o de segunda especie.

Descripción

Procedimiento de Carlson para hallar las integrales elípticas de primera y segunda especie.

F(k,θ)= ∫ 0 θ dϕ 1− k 2 sin⁡ 2 ϕ   E(k,θ)= ∫ 0 θ 1− k 2 sin⁡ 2 ϕ dϕ 

public class RD {
  static final double ERRTOL=0.05;
  static final double TINY=1.0e-25;
  static final double BIG=4.5e21;
  static final double C1=3.0/14.0;
  static final double C2=1.0/6.0;
  static final double C3=9.0/22.0;
  static final double C4=3.0/26.0;
  static final double C5=0.25*C3;
  static final double C6=1.5*C4;

  public static double  rd(double x, double y, double z){
      double alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,
sqrty, sqrtz,sum,xt,yt,zt;
      if (Math.min(x,y) < 0.0 || Math.min(x+y,z) < TINY 
|| Math.max(Math.max(x,y),z) > BIG)
          System.out.println("invalid arguments in rd");
      xt=x;
      yt=y;
      zt=z;
      sum=0.0;
      fac=1.0;
      do {
          sqrtx=Math.sqrt(xt);
          sqrty=Math.sqrt(yt);
          sqrtz=Math.sqrt(zt);
          alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
          sum += fac/(sqrtz*(zt+alamb));
          fac=0.25*fac;
          xt=0.25*(xt+alamb);
          yt=0.25*(yt+alamb);
          zt=0.25*(zt+alamb);
          ave=0.2*(xt+yt+3.0*zt);
          delx=(ave-xt)/ave;
          dely=(ave-yt)/ave;
          delz=(ave-zt)/ave;
      } while (Math.max(Math.max(Math.abs(delx),Math.abs(dely)),
Math.abs(delz)) > ERRTOL);
      ea=delx*dely;
      eb=delz*delz;
      ec=ea-eb;
      ed=ea-6.0*eb;
      ee=ed+ec+ec;
      return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee)+
delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*Math.sqrt(ave));
  }
}

public class RF {
    static final double ERRTOL=0.08;
    static final double TINY=1.5e-38;
    static final double BIG=3.0e37;
    static final double THIRD=(1.0/3.0);
    static final double C1=(1.0/24.0);
    static final double C2=0.1;
    static final double C3=(3.0/44.0);
    static final double C4=(1.0/14.0);

public  static double rf(double x, double y, double z) {
    double alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt;
    if (Math.min(Math.min(x,y),z)<0.0 || Math.min(Math.min(x+y,x+z),y+z)<TINY 
|| Math.max(Math.max(x,y),z)>BIG)
        System.out.println("invalid arguments in rf");
        xt=x;
        yt=y;
        zt=z;
        do {
            sqrtx=Math.sqrt(xt);
            sqrty=Math.sqrt(yt);
            sqrtz=Math.sqrt(zt);
            alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz;
            xt=0.25*(xt+alamb);
            yt=0.25*(yt+alamb);
            zt=0.25*(zt+alamb);
            ave=THIRD*(xt+yt+zt);
            delx=(ave-xt)/ave;
            dely=(ave-yt)/ave;
            delz=(ave-zt)/ave;
        } while (Math.max(Math.max(Math.abs(delx),Math.abs(dely)),
Math.abs(delz))> ERRTOL);
        e2=delx*dely-delz*delz;
        e3=delx*dely*delz;
        return (1.0+(C1*e2-C2-C3*e3)*e2+C4*e3)/Math.sqrt(ave);
    }
}

Las integrales elípticas de primera y segunda especie.

public class Eliptica {
  public static double segunda(double phi, double ak) {
    double cc,q,s;
    s=Math.sin(phi);
    cc=Math.cos(phi)*Math.cos(phi);
    q=(1.0-s*ak)*(1.0+s*ak);
    return s*(RF.rf(cc,q,1.0)-(s*ak)*(s*ak)*RD.rd(cc,q,1.0)/3.0);
  }
  public static double primera(double phi, double ak) {
        double s=Math.sin(phi);
        return (s*RF.rf((Math.cos(phi)*Math.cos(phi)),(1.0-s*ak)*(1.0+s*ak),1.0));
  }
}

Referencia

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition. Special functions. Sección 6.11 Elliptic Integrals and Jacobian Elliptic Functions. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java

Siguiente Anterior