0

Estoy intentando pasar éste código de Matlab a Python, es la solución numérica de la integral elíptica completa generalizada de tercer tipo, la cual es la siguiente:

introducir la descripción de la imagen aquí

Éste es el código en Matlab que tengo:

function res = integral_eliptica(kc,p,c,s)
if kc==0
    res=NaN;
    quit;
end
errtol=0.000001;
k=abs(kc);
pp=p;
cc=c;
ss=s;
em=1;
if p>0
    pp=sqrt(p);
    ss=s/pp;
else
    f=kc^2;
    q=1-f;
    g=1-pp;
    f=f-pp;
    q=q*(ss-c*pp);
    pp=sqrt(f/g);
    cc=(c-ss)/g;
    ss=-q/(g^2*pp)+cc*pp;
end
f=cc;
cc=cc+ss/pp;
g=k/pp;
ss=2*(ss+f*g);
pp=g+pp;
g=em;
em=k+em;
kk=k;
while abs(g-k)>g*errtol
    k=2*sqrt(kk);
    kk=k*em;
    f=cc;
    cc=cc+ss/pp;
    g=kk/pp;
    ss=2*(ss+f*g);
    pp=g+pp;
    g=em;
    em=k+em;
end
res=(pi/2)*(ss+cc*em)/(em*(em+pp));
end

Y esto es lo que he intentado hacer en Python:

import numpy as np
from numpy import npymath
def cel(kc,p,c,s):
  if kc==0:
      return cel = NPY_NAN
  errtol=0.000001
  k=ABS(kc)
  pp=p
  cc=c
  ss=s
  em=1.0
  
  if (p>0)
    pp=(p)**(1/2)
    ss=s/pp
  else
    f=kc*kc
    q=1.0-f
    g=1.0-pp
    f=f-pp
    q=q*(ss-c*pp)
    pp=(f/g)**(1/2)
    cc=(c-ss)/g
    ss=-q/(g*g*pp)+cc*pp
  
  f=cc
  cc=cc+ss/pp
  g=k/pp
  ss=2*(ss+f*g)
  pp=g+pp
  g=em
  em=k+em
  kk=k

  while (abs(g-k)>g*errtol)
      k=2*(kk)**(1/2)
      kk=k*em
      f=cc
    cc=cc+ss/pp
    g=kk/pp
    ss=2*(ss+f*g)
    pp=g+pp
    g=em
    em=k+em
  
  return cel=(pi/2)*(ss+cc*em)/(em*(em+pp))

Sin embargo, no logro hacer que funcione, no estoy muy familiarizado con la sintaxis de Python para lo que requiero. Otra opción que he probado para resolver la integral es usando su equivalente con R_F y R_J de scipy.special:

introducir la descripción de la imagen aquí

Pero me aparece el siguiente error:

AttributeError: module 'scipy.special' has no attribute 'elliprf'

6
  • Buen día, ¿Podrías proporcionar dos conjuntos de valores de entrada y sus respectivos resultados? Commented el 27 jun. 2022 a las 0:11
  • Basado en la fórmula que adjuntas, algo como esto debería funcionar, pero debes proporcionar datos para comprobar las funciones. Tal y cómo está, es difícil reproducir tus salidas.
    – Jaime
    Commented el 27 jun. 2022 a las 0:30
  • @Jaime Menéndez Hice un código para graficar las líneas de campo magnético de una espira, utilizo sin problemas scipy.special, pero ahora estoy haciendo lo mismo para un solenoide completo de N espiras y parte de las componentes requieren la solución de la integral eliptica completa, esto es lo que llevo ahora Commented el 27 jun. 2022 a las 1:32
  • 2
    Buen día @EstebanGonzález, la pregunta de Jaime y la mía es básicamente la misma, si nos proporcionas ejemplos de los datos de entrada y de salida podemos verificar si lo que hemos intentado hacer para resolver tu problema funciona o no. En mi caso, honestamente no quiero verificar todo el proceso a mano y hacer los cálculos, por eso si ya tienes valores de entrada y lo que debería salir de la función entonces podemos verificar si funciona correctamente o no y poder redactar una respuesta. Commented el 27 jun. 2022 a las 1:38
  • 2
    Las preguntas que buscan ayuda con la depuración deben incluir el comportamiento esperado, un problema específico o error, y el código más corto necesario para reproducirlo en la propia pregunta. Las preguntas sin un planteamiento claro del problema no son útiles para otros lectores. Véase: Cómo crear un Ejemplo mínimo, completo y verificable. Commented el 27 jun. 2022 a las 1:39

0

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

Examina otras preguntas con la etiqueta o formula tu propia pregunta.