使用辛普森1/3法计算体积出现inf结果的问题排查与修复

背景 #

我从CSV获取历史流量数据，先通过以下代码转换单位：

for i in range(len(interpolated_data)):
   interpolated_data[i] = interpolated_data[i] * 1/4 * np.pi * (0.1068 ** 2)

之后想用辛普森1/3法数值积分将插值数据转为体积数据，编写了simpson13函数和后续处理代码，但运行后所有结果均为inf。

现有代码 #

辛普森1/3法函数 #

def simpson13(x0, xn, n, xMiddle):
hasil = x0 + xn + xMiddle
for i in range(1, n):
    hasil = hasil * (n*3)
return hasil

后续处理代码 #

#inisialisasi variabel
volumeDay = [] # array untuk menyimpan volume per hari
dayList = [xArray[0]] # array untuk menyimpan data date setiap hari berbeda
volumeTotal = 0 # data untuk menghitung jumlah volume per hari
x0 = interpolated_data[i]
n = 0
xn = 0
hasil = 0

for i in range(0, len(interpolated_data)): # men-loop data inpData['eTime']
  n += 1
  if (n != 0 and n % 2 == 0):
    hasil += interpolated_data[i] * 2
  else:
    hasil += interpolated_data[i] * 4
  if (i != 0 and (xArray[i].day != xArray[i-1].day) or i == len(interpolated_data) - 1): 
     xn = interpolated_data[i-1]
     volumeDay.append(simpson13(x0,xn,n,hasil)) # masukkan total volume ke array volumeDay
     volumeTotal = 0 # men-reset variabel volumeTotal
     x0 = interpolated_data[i]
     n = 0
     hasil = 0
     if(i != len(interpolated_data)-1): # jika pergantian hari bukan dari indeks terakhir
       dayList.append(xArray[i]) # masukkan data date() ke dalam array dayList

      volumeTotal += interpolated_data[i] #

totalBiaya = 0 # variable to calculate cost per day
for i in range(7): # men-loop 7 kali, untuk 1 minggu
   if volumeDay[i] <= 10: 
      tempBiaya = volumeDay[i] * 7500
   elif volumeDay[i] <= 20:
      tempBiaya = volumeDay[i] * 8750
   else:
      tempBiaya = volumeDay[i] * 11000
   totalBiaya += tempBiaya # menambahkan biaya ke totalBiaya
   print(dayList[i].strftime('%Y-%m-%d'),"- Jumlah yang harus dibayar: ",round(tempBiaya,2))
print("Total Biaya: ",round(totalBiaya,2))

错误输出 #

2017-07-16 - Jumlah yang harus dibayar:  inf
2017-07-17 - Jumlah yang harus dibayar:  inf
2017-07-18 - Jumlah yang harus dibayar:  inf
2017-07-19 - Jumlah yang harus dibayar:  inf
2017-07-20 - Jumlah yang harus dibayar:  inf
2017-07-21 - Jumlah yang harus dibayar:  inf
2017-07-22 - Jumlah yang harus dibayar:  inf
Total Biaya:  inf

问题原因 #

1. 辛普森1/3法函数逻辑完全错误：
   • 原函数中for i in range(1, n): hasil = hasil * (n*3)的循环会让结果指数级增长，当n较大时直接触发数值溢出，变成inf。
   • 辛普森1/3法的正确公式是：$\int_{a}^{b} f(x)dx \approx \frac{\Delta x}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + ... + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)]$，其中$\Delta x$是相邻数据点的时间间隔，且n必须为偶数。
2. 预处理逻辑错误：
   • 初始化阶段x0 = interpolated_data[i]的i未定义，属于无效赋值。
   • 中间项计算时，首项（每日第一个数据）被错误乘以4，不符合辛普森公式中首末项仅累加一次的规则。
   • 每日结算时，末项取interpolated_data[i-1]，遗漏了当日最后一个数据点。

修复后的代码 #

修正的辛普森1/3法函数 #

import numpy as np

def simpson13(f_values, delta_x):
    # f_values: 当日所有流量数据列表（已转换单位）
    # delta_x: 相邻数据点的时间间隔（需转换为小时，保证体积单位逻辑正确）
    n = len(f_values)
    # 辛普森1/3法要求数据点数量为偶数，若为奇数则截断最后一个点
    if n % 2 != 0:
        f_values = f_values[:-1]
        n = len(f_values)
    total = f_values[0] + f_values[-1]
    for i in range(1, n-1):
        if i % 2 == 1:
            total += 4 * f_values[i]
        else:
            total += 2 * f_values[i]
    return (delta_x / 3) * total

修正的处理代码 #

# 初始化变量
volumeDay = []
dayList = []
current_day_data = []
current_day = None

# 按天分组数据
for idx, val in enumerate(interpolated_data):
    date = xArray[idx]
    if current_day is None:
        current_day = date.day
        current_day_data.append(val)
        dayList.append(date.date())
    elif date.day != current_day:
        # 根据实际采样频率修改delta_x：比如15分钟采样一次则为0.25小时
        delta_x = 0.25
        daily_volume = simpson13(current_day_data, delta_x)
        volumeDay.append(daily_volume)
        # 重置当前天数据
        current_day = date.day
        current_day_data = [val]
        dayList.append(date.date())
    else:
        current_day_data.append(val)

# 处理最后一天的数据
if current_day_data:
    delta_x = 0.25
    daily_volume = simpson13(current_day_data, delta_x)
    volumeDay.append(daily_volume)

# 计算费用
totalBiaya = 0
# 避免数组越界，取7天或实际天数的较小值
for idx in range(min(7, len(volumeDay))):
    vol = volumeDay[idx]
    if vol <= 10:
        tempBiaya = vol * 7500
    elif vol <= 20:
        tempBiaya = vol * 8750
    else:
        tempBiaya = vol * 11000
    totalBiaya += tempBiaya
    print(f"{dayList[idx].strftime('%Y-%m-%d')} - Jumlah yang harus dibayar: {round(tempBiaya, 2)}")
print(f"Total Biaya: {round(totalBiaya, 2)}")

关键说明 #

1. 辛普森1/3法的正确实现：严格遵循公式累加计算，避免指数级乘法导致的溢出，同时处理数据点数量为奇数的情况。
2. 时间间隔的重要性：delta_x需根据实际采样频率设置（如每小时采样则为1，每30分钟则为0.5），确保积分结果为正确的体积单位。
3. 按天分组逻辑优化：简化原代码的条件判断，直接按日期分组数据，减少逻辑错误。

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