椭圆积分Python代码报错‘cannot create mpf from array’求修复方案

错误原因 #

你遇到的cannot create mpf from array错误，核心原因是mpmath库的椭圆积分函数（ellipe、ellipk、ellippi）仅支持标量输入，但代码中k**2、n都是从numpy数组u生成的数组类型，mpmath无法将numpy数组转换为它的mpf（多精度浮点数）类型，因此报错。

修复方案 #

有两种简便的修复方式，根据需求选择：

方案1：用numpy向量化包装mpmath函数 #

通过np.vectorize()将mpmath的标量函数转换为支持数组的函数，无需改动太多代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
from mpmath import *

# 向量化mpmath的椭圆积分函数
vec_ellipe = np.vectorize(ellipe)
vec_ellipk = np.vectorize(ellipk)
vec_ellippi = np.vectorize(ellippi)

G = 6.6743*10**(-11) #引力常数
M = 10*(1.988*10**30) #黑洞质量(kg)
m = 1*(1.988*10**30) #伴星质量(kg)
Mt = M + m #总质量(kg)
q = m/M #质量比
c = 2.99792458*10**8 #光速(m/s)
Period = 10 #轨道周期(天)
P = Period*86400 #轨道周期(秒)
phi = 0.01*(np.pi/(180)) #视线倾角
t0 = Period/2 #脉冲中点
a = ((P**2*G*(Mt))/(4*np.pi**2))**(1/3) #半长轴
omega = ((G*(m + M))/a**3)**0.5 #角速度
vtran = a*omega #横向速度
Rs = (2*G*M)/c**2 #史瓦西半径
Rein = (2*Rs*a)**0.5 #爱因斯坦半径
te = (Rein/vtran)/86400 #爱因斯坦时间
u0 = (a/Rein)*phi #最近角距
pstar = ((1*(6.957*10**8))/Rein)
t = np.linspace(0,P,2500000) #时间向量
u = ((u0**2)+((((t/86400)-t0)/te))**2)**0.5 #角距离

n = ((4*u*pstar)/(u + pstar)**2)
k = ((4*n)/(4 + (u - pstar)**2))**0.5

# 使用向量化后的函数计算
Amp = (1/(2*(np.pi)))*(((((u+pstar)/(pstar**2))*np.sqrt(4+(u-pstar)**2)*vec_ellipe(k**2))-(((u-pstar)/(pstar**2))*(((8+u**2-pstar**2)/(np.sqrt(4+(u-pstar)**2)))*vec_ellipk(k**2)))+(((4*(u-pstar)**2)/(pstar**2*(u+pstar)))*(((1+pstar**2)/(np.sqrt(4+(u-pstar)**2)))*vec_ellippi(n,k**2)))))

方案2：改用Scipy的椭圆积分函数 #

Scipy的scipy.special模块提供了原生支持numpy数组的椭圆积分函数，替换mpmath的函数即可，性能比向量化更好：

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import ellipe, ellipk, ellippi  # 替换mpmath的函数

G = 6.6743*10**(-11) #引力常数
M = 10*(1.988*10**30) #黑洞质量(kg)
m = 1*(1.988*10**30) #伴星质量(kg)
Mt = M + m #总质量(kg)
q = m/M #质量比
c = 2.99792458*10**8 #光速(m/s)
Period = 10 #轨道周期(天)
P = Period*86400 #轨道周期(秒)
phi = 0.01*(np.pi/(180)) #视线倾角
t0 = Period/2 #脉冲中点
a = ((P**2*G*(Mt))/(4*np.pi**2))**(1/3) #半长轴
omega = ((G*(m + M))/a**3)**0.5 #角速度
vtran = a*omega #横向速度
Rs = (2*G*M)/c**2 #史瓦西半径
Rein = (2*Rs*a)**0.5 #爱因斯坦半径
te = (Rein/vtran)/86400 #爱因斯坦时间
u0 = (a/Rein)*phi #最近角距
pstar = ((1*(6.957*10**8))/Rein)
t = np.linspace(0,P,2500000) #时间向量
u = ((u0**2)+((((t/86400)-t0)/te))**2)**0.5 #角距离

n = ((4*u*pstar)/(u + pstar)**2)
k = ((4*n)/(4 + (u - pstar)**2))**0.5

# 直接使用scipy的函数，无需修改计算逻辑
Amp = (1/(2*(np.pi)))*(((((u+pstar)/(pstar**2))*np.sqrt(4+(u-pstar)**2)*ellipe(k**2))-(((u-pstar)/(pstar**2))*(((8+u**2-pstar**2)/(np.sqrt(4+(u-pstar)**2)))*ellipk(k**2)))+(((4*(u-pstar)**2)/(pstar**2*(u+pstar)))*(((1+pstar**2)/(np.sqrt(4+(u-pstar)**2)))*ellippi(n,k**2)))))

注意事项 #

• 方案2的性能远优于方案1，因为np.vectorize()本质是循环遍历数组元素，适合小批量数据；而Scipy的函数是底层优化过的，适合你代码中250万元素的大数组。
• 确认Scipy版本，旧版本的ellippi参数顺序可能与mpmath不同：Scipy的ellippi(n, m)对应mpmath的ellippi(n, m)，参数顺序一致，无需调整。

内容的提问来源于stack exchange，提问作者sorabella91