FFT和Romberg求积

FFT：

代码实现：

import numpy as np

# fx为初始的函数结点值
fx = np.array(list(map(float, input("f(xk):").split(' '))))
# N为函数结点个数
N = len(fx)
# p为计算次数
p = np.log2(N)
X = []
for i in fx:
    X.append(i / N)

A1 = np.zeros(N, dtype=complex)
A2 = np.zeros(N, dtype=complex)
omiga = np.zeros(N, dtype=complex)

# A1和omiga赋初值
for j in range(N):
    A1[j] = X[j]
    omiga[j] = complex(np.cos(2 * np.pi * j / N), np.sin(2 * np.pi * j) / N)

# 分奇偶计算
for q in range(1, int(p + 1)):
    # 奇数情况下
    if q % 2 == 1:
        for m in range(int(pow(2, p - q))):
            for n in range(int(pow(2, q - 1))):
                A2[m * int(pow(2, q)) + n] = A1[m * int(pow(2, q - 1)) + n] + A1[m * int(pow(2, q - 1)) + n + int(pow(2, p - 1))]
                A2[m * int(pow(2, q)) + n + int(pow(2, q - 1))] = (A1[m * int(pow(2, q - 1)) + n] - A1[m * int(pow(2, q - 1)) + n + int(pow(2, p - 1))]) * omiga[m * int(pow(2, q - 1)) % N]
    else:
        for m in range(int(pow(2, p - q))):
            for n in range(int(pow(2, q - 1))):
                A1[m * int(pow(2, q)) + n] = A2[m * int(pow(2, q - 1)) + n] + A2[m * int(pow(2, q - 1)) + n + int(pow(2, p - 1))]
                A1[m * int(pow(2, q)) + n + int(pow(2, q - 1))] = (A2[m * int(pow(2, q - 1)) + n] + A2[m * int(pow(2, q - 1)) + n + int(pow(2, p - 1))]) * omiga[m * int(pow(2, q - 1)) % N]

# 输出FFT后结果
if p % 2 == 0:
    print(A1)
else:
    print(A2)

实验结果：

Romber求积：

代码实现：

import numpy as np


# 被积函数
def fun(x):
    f = pow(x, 1.5)
    return f


# 求梯形值
def T2n(a, b, n, Tn):
    h = (b - a) / n  # 步长
    sum = 0.
    for k in range(n):
        sum += fun(a + (k + 0.5) * h)
    T2n = Tn / 2. + sum * h / 2.
    return T2n


# 求加速值
def romberg(max, a, b, epsilon):  # max为计算最大次数，a、b为积分下、上限，epsilon为限定误差
    tm = np.zeros(max, dtype=float)  # 第m行元素
    tm1 = np.zeros(max, dtype=float)  # 第m+1行元素
    tm[0] = (b - a) * (fun(a) + fun(b)) / 2.  # 初始值
    print(tm)
    k = 0
    err = 1
    while (err > epsilon and k < max - 1):  # 当误差小于预定误差，或计算次数大于最大次数时停止
        n = 2 ** k
        m = 1
        tm1[0] = T2n(a, b, n, tm[0])
        while (err > epsilon and m <= (k + 1)):  # 当误差小于预定误差，或本行全部计算完毕后停止
            tm1[m] = tm1[m - 1] + (tm1[m - 1] - tm[m - 1]) / (4. ** m - 1)
            result = tm1[m]
            err1 = abs(tm1[m] - tm[m - 1])  # 对角元素差
            err2 = abs(tm1[m] - tm1[m - 1])  # 同行前后两元素差
            err = min(err1, err2)
            m += 1
        tm = np.copy(tm1)  # 下移一行
        k += 1
        print(tm)
    return result


if __name__ == '__main__':
    f1 = romberg(6, 0, 1, 1.e-10)
    print(f1)

实验结果：