看起来你现在的问题是把固定上限的积分和随t变化的积分上限搞混了——你的代码当前计算的是从0到Tmax的单一积分值，而你需要的是每个时刻t∈[0,Tmax]对应的积分结果（也就是A(t)数组，每个元素对应一个t点的积分值）。

问题分析 #

你原代码里的逻辑是计算：
$$A = \int_0^{Tmax} e^{-\gamma(Tmax - s)} B(s) ds$$
这是一个固定上限Tmax的积分，所以返回的是标量。但你实际需要的应该是对每个t_i，计算：
$$A(t_i) = \int_0^{t_i} e^{-\gamma(t_i - s)} B(s) ds$$
这是一个随t_i变化的积分，每个t_i对应一个积分结果，最终要得到长度为N的数组。

高效解决方案：递推法（O(N)复杂度） #

直接对每个t_i计算梯形法则的话，时间复杂度是O(N²)，对于N=1e5来说完全不可行。我们可以利用积分的拆分性质和指数衰减的特点，推导递推公式来高效计算：

把A(t_i)拆分为两部分：

1. 从0到t_{i-1}的积分，通过指数衰减转换到t_i时刻
2. 从t_{i-1}到t_i的区间积分，用梯形法则近似

最终得到递推式：
$$A(t_i) = e^{-\gamma h} \cdot A(t_{i-1}) + h \cdot \left( \frac{1}{2}B(t_i) + \frac{1}{2}B(t_{i-1})e^{-\gamma h} \right)$$
其中h是步长（h=Tmax/N）。

对应的代码实现：

import numpy as np

gamma = 0.1  # parameter
Tmax = 100
N = 100000  # number of steps
t = np.linspace(0, Tmax, N)
h = Tmax / N  # 计算步长
B = np.random.rand(N)  # 生成对应每个s=t_j的随机数组

# 初始化A数组，长度和t一致
A = np.zeros_like(t)
A[0] = 0.0  # 初始条件：t=0时积分结果为0

# 预计算指数因子，避免重复计算提升效率
exp_gamma_h = np.exp(-gamma * h)

# 递推计算每个t_i对应的A(t_i)
for i in range(1, N):
    A[i] = exp_gamma_h * A[i-1] + h * (0.5 * B[i] + 0.5 * B[i-1] * exp_gamma_h)

这段代码运行后，A就是你需要的从t=0到t=Tmax的函数数组，每个元素对应一个时刻的积分结果。

为什么原代码不对？ #

原代码中：

• 你把积分上限固定成了Tmax，而不是当前的t_i
• 循环里累加的是针对Tmax的被积函数，最后只计算了一个标量结果
• 没有为每个t_i单独计算对应的积分区间

（可选）直观但低效的实现（仅用于理解） #

如果你想直观看到梯形法则的直接应用（不推荐用于N=1e5的情况，因为O(N²)太慢），可以参考这段代码：

A = np.zeros_like(t)
for i in range(N):
    if i == 0:
        A[i] = 0.0
        continue
    # 取0到t[i]对应的所有s点和B值
    s_points = t[:i+1]
    b_vals = B[:i+1]
    # 计算被积函数：e^(-gamma*(t[i]-s)) * B(s)
    integrand = np.exp(-gamma * (t[i] - s_points)) * b_vals
    # 梯形法则计算积分
    A[i] = h * (0.5 * integrand[0] + np.sum(integrand[1:-1]) + 0.5 * integrand[-1])

但注意，这段代码对于N=1e5来说会运行极久，所以优先使用递推法。

内容的提问来源于stack exchange，提问作者J.Agusti