Python中含可变边界的定积分计算问题求助

我正在计算含可变边界的定积分，需要获取数值输出但一直失败，以下是我的两种尝试方案，请问如何让其中一种正常运行？

第一种计算方式 #

import numpy as np
from sympy import *

Nt = 17                                         
alpha = .99
t = np.linspace(0, .85, Nt)     
s = np.linspace(0, .85, Nt) 

for k in range(Nt): 
  for i in range(k): 
    Int = integrate( (t[k] - s) ** - int(alpha), (s, t[i], t[i + 1] )) 
    print(Int) 

报错信息 #

ValueError: Invalid limits given: ((array([0.      , 0.053125, 0.10625 , 0.159375, 0.2125  , 0.265625,
       0.31875 , 0.371875, 0.425   , 0.478125, 0.53125 , 0.584375,
       0.6375  , 0.690625, 0.74375 , 0.796875, 0.85    ]), 0.0, 0.053125),)

第二种计算方式 #

import numpy as np
from sympy import *
from numba import jit, prange

Nt = 17                                         
alpha = .99
t = np.linspace(0, .85, Nt)     
s = np.linspace(0, .85, Nt) 

@jit(nopython=True)
def Int(alpha):
  Int = 0
  for k in prange(Nt): 
    for i in prange(k): 
      Int = Int + integrate( (t[k] - s) ** - int(alpha), (s, t[i], t[i + 1] )) 
      print(Int)   
  return Int

此方法无任何输出。

编辑说明：修改后的第一种方法 #

我对第一种方法做了小幅修改，得到了一些数值结果，但不确定结果是否正确，欢迎各位提供意见。

修改后的代码 #

Nt = 17                                         
alpha = .99
t = np.linspace(0, .85, Nt)     
s = symbols('s')

for k in range(Nt): 
  for i in range(k): 
    Int = integrate((t[k] - s) ** - int(alpha), (s, t[i], t[i + 1] )) 
    print(Int) 

输出结果 #

0.0531250000000000
0.0531250000000000
0.0531250000000000
0.0531250000000000...

问题分析与解决方案 #

1. 原第一种方法报错原因 #

你把积分变量s定义成了numpy数组，但sympy的integrate函数要求积分变量是符号变量，而非数值数组，这就导致积分边界和变量冲突，触发了无效边界的错误。

2. 原第二种方法无输出原因 #

numba的nopython=True模式只支持纯数值计算的Python代码，完全不兼容sympy的符号计算函数（比如integrate）。同时prange用于并行数值循环，在这里也无法发挥作用，导致函数内部代码根本没有执行，所以没有输出。

3. 修改后代码的结果问题 #

修改后你把s改成了符号变量，这部分是对的，但犯了一个关键错误：-int(alpha)。因为alpha=0.99，int(alpha)会直接截断为0，被积函数变成(t[k]-s)^0=1，积分结果自然就是积分区间的长度（t[i+1]-t[i] = 0.85/(17-1)=0.053125），这显然不是你想要的结果——你应该用-alpha而非-int(alpha)。

正确的实现方式 #

方式一：用sympy符号积分后转数值

import numpy as np
from sympy import symbols, integrate, N

Nt = 17                                         
alpha = 0.99
t = np.linspace(0, 0.85, Nt)     
s = symbols('s')

for k in range(Nt): 
    for i in range(k):
        # 正确使用alpha而非int(alpha)
        integral_expr = integrate((t[k] - s) ** (-alpha), (s, t[i], t[i + 1]))
        # 转换为数值输出
        numeric_result = N(integral_expr)
        print(numeric_result)

方式二：用scipy数值积分（更适合大规模计算）

如果不需要符号解析解，直接用scipy的数值积分函数quad效率更高：

import numpy as np
from scipy.integrate import quad

Nt = 17                                         
alpha = 0.99
t = np.linspace(0, 0.85, Nt)     

def integrand(s, t_k, alpha):
    return (t_k - s) ** (-alpha)

for k in range(Nt): 
    for i in range(k):
        result, _ = quad(integrand, t[i], t[i+1], args=(t[k], alpha))
        print(result)

这两种方式都能得到正确的数值结果，其中方式二更适合大规模循环计算，速度更快。

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