基于Python辛普森法则实现傅里叶级数系数计算与绘图

我来帮你一步步解决这两个问题——先搞定通用的辛普森法则积分函数，再用它计算奇函数的傅里叶级数系数并可视化结果。

辛普森法则通过将积分区间划分为偶数个小段，用二次曲线拟合每个小段来近似积分值。下面是一个可以适配任意可调用函数、自定义区间和分割数的通用实现：

def simpsons_rule(f, a, b, n):
    # 辛普森法则要求分割数n必须为偶数，先做校验
    if n % 2 != 0:
        raise ValueError("分割数n必须是偶数，请调整参数")
    h = (b - a) / n  # 每个小段的长度
    # 计算奇数位置和偶数位置的函数值之和
    sum_odd_terms = sum(f(a + i*h) for i in range(1, n, 2))
    sum_even_terms = sum(f(a + i*h) for i in range(2, n-1, 2))
    # 套用辛普森公式计算积分
    integral_result = (h/3) * (f(a) + f(b) + 4*sum_odd_terms + 2*sum_even_terms)
    return integral_result

测试验证 #

比如用它计算f(x) = x²在[0,1]上的积分（理论值为1/3≈0.3333）：

def test_func(x):
    return x**2

print(simpsons_rule(test_func, 0, 1, 100))  # 输出约0.33333333333333337，符合预期

因为你提到目标函数是奇函数，周期为T，根据傅里叶级数的性质：

• 直流分量a₀ = 0
• 所有余弦项系数aₖ = 0（k≥1）
• 仅需计算正弦项系数bₖ，公式为：

bₖ = (2/T) ∫_{-T/2}^{T/2} f(x)·sin(2πk x/T) dx

又因为奇函数乘奇函数是偶函数，积分区间可以简化为[0, T/2]，公式变为：

bₖ = (4/T) ∫_{0}^{T/2} f(x)·sin(2πk x/T) dx

下面是调用辛普森法则计算bₖ的代码，同时附带一个示例奇函数（周期为2的方波）：

import numpy as np
import matplotlib.pyplot as plt

# 定义示例奇函数：周期T=2的方波，你可以替换成自己的目标函数
def target_odd_func(x):
    x_mod = x % 2  # 取周期内的位置
    return 1 if x_mod < 1 else -1

# 计算单个b_k系数
def compute_fourier_bk(f, T, k, n=1000):
    integrand = lambda x: f(x) * np.sin(2 * np.pi * k * x / T)
    # 调用辛普森法则计算[0, T/2]区间的积分
    integral = simpsons_rule(integrand, 0, T/2, n)
    return (4 / T) * integral

# 计算前N个正弦项系数
N = 10
T = 2
b_coeffs = [compute_fourier_bk(target_odd_func, T, k) for k in range(1, N+1)]

# 打印结果
print(f"前{N}个傅里叶正弦项系数：")
for idx, bk in enumerate(b_coeffs, 1):
    print(f"k={idx}: {bk:.4f}")

用计算得到的bₖ重构原函数，对比原函数和近似曲线：

# 生成x轴数据，覆盖3个周期
x = np.linspace(-3, 3, 1000)
# 计算傅里叶级数近似值
f_approx = np.zeros_like(x)
for k in range(1, N+1):
    f_approx += b_coeffs[k-1] * np.sin(2 * np.pi * k * x / T)

# 绘图
plt.figure(figsize=(10, 6))
plt.plot(x, [target_odd_func(xi) for xi in x], label='原奇函数（方波）', linestyle='--', alpha=0.7)
plt.plot(x, f_approx, label=f'傅里叶级数近似（前{N}项）', color='crimson')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('奇函数的傅里叶级数近似效果')
plt.legend()
plt.grid(True)
plt.show()

• 辛普森法则的分割数n越大，积分精度越高；对于不连续的奇函数（比如方波），建议取较大的n（如1000）来保证系数计算准确。
• 如果你的目标奇函数不是方波，只需要替换target_odd_func函数即可，记得保持其奇函数性质和周期T。

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