Python应用梯形法数值计算及相邻结果累加实现方案咨询

现有代码问题说明 #

你原有代码存在下标逻辑问题：当i=0时i-1=-1会取数组末尾的元素，得到的第一个输出结果不符合区间计算逻辑，正确的区间遍历范围应为1 ≤ i ≤ len(P)-1。

实现方案 #

步骤1：正确计算所有梯形区间的面积值 #

你要计算的被积函数为 f(P) = 2P/(mi_g * z)，每个区间的梯形面积公式为 (f(P[i-1]) + f(P[i])) * dP / 2，先将所有区间的面积存入列表。

步骤2：遍历面积列表计算相邻累加值 #

遍历面积列表，将相邻两个元素相加得到你需要的mp序列。

可运行代码 #

import numpy as np

P=np.array([0.0,150,300,450,600,750,900,1050,1200,1350,1500,1650,1800,1950,2100,2250,2400,2550,2700,2850,3000,3150])
mi_g=np.array([0.0,0.01238,0.01254,0.01274,0.01303,0.01329,0.0136,0.01387,0.01428,0.01451,0.01485,0.0152,0.01554,0.01589,0.0163,0.01676,0.01721,0.01767,0.0183,0.01862,0.01911,0.01961])
z=np.array([0.0,0.9856,0.9717,0.9582,0.9453,0.9332,0.9218,0.9112,0.9016,0.8931,0.8857,0.8795,0.8745,0.8708,0.8684,0.8671,0.8671,0.8683,0.8705,0.8738,0.8780,0.8830])
dP=150

# 计算被积函数数组
f_arr = 2 * P / (mi_g * z)
# 存储所有梯形区间面积
trap_areas = []
for i in range(1, len(P)):
    area = (f_arr[i-1] + f_arr[i]) * dP / 2
    trap_areas.append(area)

# 计算相邻累加的mp序列
mp_list = []
for i in range(len(trap_areas)-1):
    mp = trap_areas[i] + trap_areas[i+1]
    mp_list.append(round(mp)) # 按你示例的精度取整

# 输出结果
for idx, val in enumerate(mp_list, 1):
    print(f"mp{idx} = {val}")

如果要简化代码可以用numpy的向量化操作，无需显式循环：

# 向量化计算所有梯形面积
trap_areas = np.convolve(f_arr, np.array([0.5*dP, 0.5*dP]), mode='valid')
# 向量化计算相邻元素和
mp_list = np.convolve(trap_areas, np.array([1, 1]), mode='valid').round().astype(int)

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