数值分析-数值积分

数值积分上机课写的 Python 代码，实现数值积分。

import numpy as np
from scipy import integrate


def data_segmentation(x_range, n):
    X = np.linspace(x_range[0], x_range[1], n + 1)
    return X
    
    
def compound_trapezoidal_quadrature(f, x_range, n):
    if x_range[0] >= x_range[1]:
        raise ValueError('x_range[0] must be less than x_range[1]', x_range)
    X = data_segmentation(x_range, n)
    Y = f(X)
    for i in range(len(Y)):
        if np.isnan(Y[i]):
            Y[i] = f(X[i] + 1e-5)
    h = (x_range[1] - x_range[0]) / n
    return h / 2 * (Y[0] + 2 * sum(Y[1: -1]) + Y[-1])
    

def compound_simpson_quadrature(f, x_range, n):
    if x_range[0] >= x_range[1]:
        raise ValueError('x_range[0] must be less than x_range[1]', x_range)
    X = data_segmentation(x_range, n)
    Y = f(X)
    for i in range(len(Y)):
        if np.isnan(Y[i]):
            Y[i] = f(X[i] + 1e-5)
    h = (x_range[1] - x_range[0]) / n
    return h / 6 * (Y[0] + 4 * sum(f(X[: -1] + h/2)) + 2 * sum(Y[1: -1]) + Y[-1])


def romberg_quadrature(f, x_range, err):
    if x_range[0] >= x_range[1]:
        raise ValueError('x_range[0] must be less than x_range[1]', x_range)
    h = x_range[1] - x_range[0]
    T = []
    T.append([compound_trapezoidal_quadrature(f, x_range, 1)])
    k = 1
    while True:
        T.append([compound_trapezoidal_quadrature(f, x_range, 2 ** k)])
        for i in range(1, k + 1):
            a = 4 ** i
            T[k].append(a / (a - 1) * T[k][i - 1] - 1 / (a - 1) * T[k - 1][i - 1])
        if abs(T[k][-1] - T[k - 1][-1]) < err:
            return [T[k][-1], k]
        else:
            k += 1