python - 使用辛普森规则在 Python 中绘制傅里叶级数系数

Question

我想 1. 将辛普森规则表达为 python 中集成的通用函数，2. 用它来计算和绘制函数的傅里叶级数系数​​。

我已经为辛普森规则窃取并修改了这段代码，它似乎可以很好地集成简单的功能，例如, 或

给定 period ，傅立叶级数系数​​计算为：

其中 k = 1,2,3,...

我很难弄清楚如何表达。我知道，因为这个函数很奇怪，但我希望能够为其他函数计算它。

到目前为止，这是我的尝试：

import matplotlib.pyplot as plt
from numpy import *

def f(t):
    k = 1
    for k in range(1,10000): #to give some representation of k's span
        k += 1
    return sin(t)*sin(k*t)

def trapezoid(f, a, b, n):
    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for j in range(1, n):
        s += f(a + j*h)
    s += f(b)/2.0
    return s * h

print trapezoid(f, 0, 2*pi, 100)

这根本没有给出 0 的正确答案，因为它随着 k 的增加而增加，而且我确信我正在用 for 循环的隧道视觉接近它。特别是我的困难在于说明函数，以便将 k 读取为 k = 1,2,3,...

不幸的是，我遇到的问题没有指定要绘制的系数是什么，但我假设它是针对 k 的。

Answer 1

这是一种方法，如果您想运行自己的积分或傅立叶系数确定而不是使用numpyorscipy的内置方法：

import numpy as np

def integrate(f, a, b, n):
    t = np.linspace(a, b, n)
    return (b - a) * np.sum(f(t)) / n

def a_k(f, k):
    def ker(t): return f(t) * np.cos(k * t)
    return integrate(ker, 0, 2*np.pi, 2**10+1) / np.pi

def b_k(f, k):
    def ker(t): return f(t) * np.sin(k * t)
    return integrate(ker, 0, 2*np.pi, 2**10+1) / np.pi

print(b_k(np.sin, 0))

这给出了结果

0.0

---

另一方面，梯形积分对于均匀的时间间隔不是很有用。但如果你愿意：

def trap_integrate(f, a, b, n):
    t = np.linspace(a, b, n)
    f_t = f(t)
    dt = t[1:] - t[:-1]
    f_ab = f_t[:-1] + f_t[1:]
    return 0.5 * np.sum(dt * f_ab)

np.trapz如果您想使用预建功能，还有。同样，还有scipy.integrate.trapz