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

Question

我想 1. 将辛普森法则表示为 python 中的一般积分函数，以及 2. 用它来计算和绘制函数 .[=20 的傅里叶级数系数=]

我盗用并改编了this code for Simpson's Rule, which seems to work fine for integrating simple functions such as ，或

给定周期，傅立叶级数系数计算如下：

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

我很难弄清楚如何表达 . I'm aware that 因为这个函数很奇怪，但我希望能够为其他函数计算它。

这是我目前的尝试：

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

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

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 if you want to use pre-builtin functionality. Similarly, there's also scipy.integrate.trapz

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

Plotting Fourier Series coefficients in Python using Simpson's Rule

python

numerical-integration

continuous-fourier

simpsons-rule