Scipy修正贝塞尔函数k0运行报错：复数输入不支持问题求助

问题根源 #

scipy.special.k0和k1仅支持实数输入，而你的计算中beta*rw是复数数组，触发了类型不兼容错误。

解决方案 #

使用scipy.special.kv函数替代，它是修正贝塞尔函数的通用实现，支持复数参数：

• 阶数为0时，kv(0, z)等价于k0(z)（但支持复数）
• 阶数为1时，kv(1, z)等价于k1(z)（但支持复数）

同时修正代码中未定义的phase函数，改用numpy官方的np.angle计算相位。

修改后的完整代码 #

import numpy as np
from scipy.special import kv  # 替换原k0、k1的导入

freqM2 = 1.932274 /86400 #cps
rw = 10*1e-3 
rc = rw #[m]
omegaM2 = 2*np.pi*freqM2

leakage = np.arange(1e-11,1e-3, 1e-6)
BKu = 40*1000000000

eo = 2.5*1e-7 * 0.045

def leaky_model(leakage, S, T, omega, BKu, rc, rw, eo, rho=1000, g=9.81):
    
    # specific leakage = Kprime/bprime
    i= 1j
    beta = ((leakage/T) + i*omega*S/T )**(1/2)
    
    # 用kv(0, ...)替代k0，kv(1, ...)替代k1
    upsilon = 1 + (rc/rw)**2 * (i*omega*rw / (2*T*beta)) * kv(0, beta*rw)/kv(1, beta*rw)
    
    hwo = i*omega*S/((i*omega*S + leakage)*upsilon) * (BKu*eo / (rho*g))  # 修正括号确保运算优先级正确
    
    A = np.abs(hwo/((BKu*eo / (rho*g))))
    eta = np.angle(hwo/(BKu*eo / (rho*g)))  # 用np.angle替代未定义的phase
    
    return A, eta


A1, eta1 = leaky_model(leakage, S=1e-2, T =1e-4, omega=omegaM2, BKu=BKu, rc=rc, rw=rw, eo=eo, rho=1000, g=9.81)

额外说明 #

1. 修正了BKu*eo / rho*g的括号，确保运算优先级符合物理公式逻辑（原代码会先计算eo/rho再乘g，实际应为eo/(rho*g)）
2. kv函数支持复数数组的向量化运算，完全兼容原代码的数组计算逻辑

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