数值分析————插值积分

这里简单介绍几种数值积分的python实现，具体数学原理后面补上。

import math
import numpy as np

three_x = [-0.5773503, 0.5773503]
three_A = [1.0000000]

five_x = [-0.9061798, -0.5384693, 0, 0.5384693, 0.9061798]
five_A = [0.2369269, 0.4786287, 0.5688889]

def f_x(x):
    return 1.0 / x

def gauss_legendre5(a, b):
    addt = (a + b) / 2
    coef = (b - a) / 2
    return f_x(addt + coef * five_x[0]) * five_A[0] + f_x(addt + coef * five_x[1]) * five_A[1] + f_x(addt + coef * five_x[2]) * five_A[2] + f_x(addt + coef * five_x[3]) * five_A[1] + f_x(addt + coef * five_x[4]) * five_A[0]

def gauss_legendre2(a, b):
    addt = (a + b) / 2
    coef = (b - a) / 2
    return coef * (f_x(addt + coef * three_x[0]) * three_A[0] + f_x(addt + coef * three_x[1]) * three_A[0])

def composite_gauss_legendre2(a, b, n):
    h = b - a
    sn = 0.0
    left = a
    right = a + (h / n)
    for i in range(n):
        sn += gauss_legendre2(left, right)
        left = right
        right = left + (h / n)
    sn = sn * h / 2
    return sn
    


def simpson(left, h, n):
    '''
        simpson插值积分
    '''
    sn = 0.0
    for i in range(n):
        right = left + h
        middle = 0.5 * (left + right)
        sn += f_x(left) + 4.0 * f_x(middle) + f_x(right)
        left = right
        print(sn)
    return h * sn / 6.0;

def composite_trapezoidal(left, h, n):
    '''
        复化梯形插值积分
    '''
    sn = 0.0
    for i in range(n):
        right = left + h
        sn += f_x(left) + f_x(right)
        left = right
    return h * sn / 2.0;

def romberg(left, h, err):
    '''
        romberg插值积分，一种递推式的积分
    '''
    T = np.zeros((5,5))
    
    for i in range(5):
        n = pow(2, i)
        T[i][0] = composite_trapezoidal(left, h/n, n)
    
    flag = 0
    for i in range(5):
        for j in range(i):
            binary_times = pow(4, j+1)
            T[i][j+1] = (binary_times * T[i][j] - T[i-1][j]) / (binary_times - 1)
        if i > 0 and abs(T[i][i] - T[i-1][i-1]) <= err:
            flag = i
            break
    return T[flag][flag]