我试图创建一个软件，通过将白矮星分成一堆层并对其进行粗略的物理模拟来计算白矮星内部的密度分布。我写了一些函数，试图进行矢量化，但我总是得到错误TypeError: return arrays must be of ArrayType。
从我在这个网站上的其他问题中收集到的信息来看，罪魁祸首可能是在不同类型的实体上运行操作。所以我尝试将frompyfunc()语句重写为vectorize语句，并尝试将类型设置为np.double，"double"和"float64"，但我最终得到了“无效类型”。是的，在执行此操作时，我记得要去掉ninput和noutput参数。
我还尝试将所有数组的类型设置为double，并澄清在我的原始函数中哪些变量是全局变量，但这些都没有帮助。
这就是

def selfBelow(m,r,x):   #The gravitational self-force of a layer below the border
    if x>r:
        return 3*G*m**2*x*(5*r**3 + 6*r**2*x + 3*r*x**2 + x*3) / (5*(r**2 + r*x + x**2)**3)
    else:
        return 1/0 #Return an error in this case because we need x to be greater than r

这是我对它的矢量化尝试：FSB = np.frompyfunc(selfBelow,2,1)
这是我尝试调用的向量化函数：forces[0]-=FSB(masses[0], 0, radii[0])
当我运行它时，我得到这个错误：

Traceback (most recent call last):
  File "C:\...\White Dwarf.py", line 78, in <module>
    forces[0]-=FSB(masses[0], 0, radii[0])
TypeError: return arrays must be of ArrayType

下面是完整的代码，以防你需要：

import numpy as np
import matplotlib as plt

h=1.054571817e-34   #Reduced Planck constant
G=6.67430e-11   #Gravitational constant
m_e=9.10938370153e-31   #Mass of an electron in kg
u=1.66053907e-27    #Atomic mass unit in kg
m_Sun=1.9885e30
r_Earth=6.3781e6
n=2     #Number of baryons per electron
C=100 #Number of layers

def volume(r, R):
    return 4*np.pi/3*(R**3-r**3)

def selfBelow(m,r,x):   #The gravitational self-force of a layer below the border
    if x>r:
        return 3*G*m**2*x*(5*r**3 + 6*r**2*x + 3*r*x**2 + x*3) / (5*(r**2 + r*x + x**2)**3)
    else:
        return 1/0

def selfAbove(m,x,R):   #The gravitational self-force of a layer abpve the border
    return 9*G*m**2*x**2*(R**2 + 3*R*x + x**2)/(10*(R**2+R*x+x**2)**3)

def elseBelow(M,m,r,x):                             #The gravitational force energy from the cumulative masses on the layer below
    if x>r:
        return 3*G*M*m*x*(2*r+x)/(2*(r**2+r*x+x**2)**2)
    else:
        return 1/0

def elseAbove(M,m,x,R):                             #The gravitational force from the cumulative masses on the layer below
    if x>r:
        return 3*G*M*m*x*(2*R+x)/(2*(R**2+R*x+x**2)**2)
    else:
        return 1/0

def fermiBelow(n,r,x):    #The fermi pressure from the layer below
    if x>r:
        return -(9*pow(3/2,1/3)*pow(np.pi,2/3)*h**2*x**2)/(10*m_e)*pow(n/(x**3-r**3),5/3)
    else:
        return 1/0

def fermiAbove(n,x,R):
    if R>x:
        return -(9*pow(3/2,1/3)*pow(np.pi,2/3)*h**2*x**2)/(10*m_e)*pow(n/(x**3-r**3),5/3)
    else:
        return 1/0


V = np.frompyfunc(volume,2,1)
FSB = np.frompyfunc(selfBelow,2,1)
FSA = np.frompyfunc(selfAbove,2,1)
FEB = np.frompyfunc(elseBelow,2,1)
FEA = np.frompyfunc(elseAbove,2,1)
FFB = np.frompyfunc(fermiBelow,2,1)
FFA = np.frompyfunc(fermiAbove,2,1)

sunMasses=1

M=m_Sun*sunMasses
Particles=M/(n*u)

startRadius = h**2/(G*M**2*m_e)*pow(Particles,5/3)*pow(9*np.pi/4,2/3)
radii = startRadius*np.arange(1,C+1)/C
masses = np.ndarray(C, np.double)
masses[0:C] = M*V(np.arange(0,C),np.arange(1,C+1))/volume(0,C)
cumuMasses  = np.ndarray(C, np.double)
cumuMasses[0:C] = M*V(0,np.arange(1,C+1))/volume(0,C) #Cumulative mass of all layers below the border
N = np.ndarray(C,np.double)
N[0:C] = Particles*V(np.arange(0,C),np.arange(1,C+1))/volume(0,C)
forces      =np.zeros(C, np.double)
print(radii.dtype, masses.dtype, cumuMasses.dtype, N.dtype, forces.dtype)

for i in range(C):
    print(i)
    forces.fill(0)
    forces[0]-=FSB(masses[0], 0, radii[0]) #Special case because radii does not include 0
    forces[1:C]-=FSB(masses[1:C], radii[0:-1], radii[1:C]) #The 0:-1 index means up to, but not including, the last element
                                                        #0 is not included for reason above    
    forces[0:-1]-=FSA(masses[1:C], radii[0:-1], radii[1:C])

    forces[1:C]-=FEB(cumuMasses[0:-1], masses[1:C], radii[0:-1], radii[1,C])
    forces[0:-1]-=FEA(cumuMasses[0:-1], masses[1:C], radii[0:-1], radii[1,C])
    
    forces[0]-=FFB(N[0], 0, radii[0]) #Special case because radii does not include 0
    forces[1:C]-=FFA(N[1:], radii[0:-1], radii[1:C])

#print(forces)

我怀疑罪魁祸首可能是我在开始时声明的变量，这些变量可能都是Python的普通对象，但我不知道如何修复。