【数值分析实验】数值积分：复化梯形、复化辛普森、龙贝格求积公式、自适应求积公式（python）

数值积分

调包

import math
import numpy as np
import matplotlib.pyplot as plt

fStr为被积函数函数名（str）

梯形求积

#梯形求积
def Trapezium(a,b,fStr):
    return (b - a) / 2 * (globals()[fStr](a) +  globals()[fStr](b))

复化梯形

#复化梯形
def TrapeziumComplex(a,b,epsilon,fStr):
    def F(x):
        return globals()[fStr](x)
    n = 0
    h0 = b - a
    T0 = h0 / 2 * (F(a) + F(b))
    while True:
        n += 1
        h1 = h0 / 2
        m = 2**(n - 1)
        sum = 0
        for k in range(1,m + 1):
            sum += F(a + (2 * k - 1) * h1)
        T1 = T0 / 2 + h1 * sum
        if abs(T1 - T0) < epsilon:
            break
        T0 = T1
        h0 = h1
    return T1,n

辛普森求积

#辛普森求积
def Simpson(a,b,fStr):
    return (b - a) / 6 * (globals()[fStr](a) +  globals()[fStr](b) + 4 * globals()[fStr]((a + b) / 2))

复化辛普森

#复化辛普森
def SimpsonComplex(a,b,epsilon,fStr):
    def F(x):
        return globals()[fStr](x)
    n = 0
    h0 = (b - a) / 2
    S0 = h0 / 3 * (F(a) + 4 * F((a + b) / 2) + F(b))
    while True:
        n += 1
        h1 = h0 / 2
        m = 2**(n - 1)
        sum1 = 0
        sum2 = 0
        for k in range(1,m * 2 + 1):
            sum1 += F(a + (2 * k - 1) * h1)
            if k <= m:
                sum2 += F(a + (4 * k - 2) * h1)
        S1 = S0 / 2 + h1 / 3 * (4 * sum1 - 2 * sum2)
        if abs(S1 - S0) < epsilon:
            break
        S0 = S1
        h0 = h1
    return S1,n

龙贝格求积

#龙贝格算法
def Romberg(a,b,epsilon,fStr):
    def F(x):
        return globals()[fStr](x)
    # 存放求积分范围
    h = b - a
    # 用于存放 T S C R  的计算结果
    T = [[j for j in [0, 0, 0, 0]] for i in range(4)]
    # 用于计行数
    i = 1
    T_1 = h * (F(b) + F(a)) / 2
    T[0][0] = T_1
    T[1][0] = T_1 / 2 + F((b - a) / 2) * h / 2
    T[1][1] = (T[1][0] * 4 ** i - T[0][0]) / (4 ** i - 1)
    # 精度达不到设定值时继续执行
    # 当行数小于四时，根据行数选择计算 T S C R 中的哪几个
    while abs(T[i][i]-T[i][i-1]) > epsilon:
        i += 1
        sum = 0
        if i == 4 :
            break
        for j in range(2 ** (i - 1)):
            sum += F(a + (2 * j + 1) * h / 2 ** i)
        T[i][0] = T[i - 1][0] / 2 + sum * h / 2 ** i
        for j in range(1, i + 1):
            T[i][j] = (T[i][j - 1] * 4 ** i - T[i - 1][j - 1]) / (4 ** i - 1)
    # 当行数大于四时，计算全部的 T S C R
    while abs(T[i-1][-1] - T[i-1][-2]) > epsilon:
        sum = 0
        T.append([])
        for j in range(2 ** (i - 1)):
            sum += F(a + (2 * j + 1) * h / 2 ** i)
        T[i].append(T[i - 1][0] / 2 + sum * h / 2 ** i)
        for j in range(1, 4):
            T[i].append((T[i][j - 1] * 4 ** i - T[i - 1][j - 1]) / (4 ** i - 1))
        i += 1
    return T[-1][-1],i

自适应求积

#自适应求积
def Adaptive(a,b,epsilon,fStr):
    def F(x):
        return globals()[fStr](x)
    h = (a + b) / 2
    ST = Simpson(a,b,fStr)
    SL = Simpson(a,h,fStr)
    SR = Simpson(h,b,fStr)
    if(abs(SL + SR - ST) <= 15.0 * epsilon):
        return SL + SR + (SL + SR - ST) / 15.0
    return Adaptive(a,h,epsilon/2.0,fStr) + Adaptive(h,b, epsilon/2.0,fStr)