数值积分初步

复化求积公式、Gauss-Legendre求积公式 Python实现

数值分析课程作业，放在博客里纪念一下，方便以后查看，若要转载请注明出处。

1. 题目分析

简单地使用Python，以面型对象的代码风格实现复化梯形公式、复化Simpson公式、Gauss-Legendre公式来解决如下积分问题:

\[I = \int_0^\theta\frac{\cos{x}}{\sqrt{x}}\mathrm{d}x \]

可以注意到这是一个瑕积分，不能直接去使用以上数值积分方法，先进行换元，令\(t=x^2\)可得\(I=2\int_0^\sqrt{\theta}\cos{t^2}\mathrm{dt}\)，现在其不是广义积分了，故转而对其进行数值积分。

2. 公式选用：

2.1 复化梯形公式

\[\int_a^bf(x)\mathrm{d}x \approx \sum_{k=1}^n \frac{h}{2}[f(x_{k-1})+f({x_k})=\frac{h}{2}[f(a)+\sum_{k=1}^{n}f(x_k)+f(b)] \]

其误差估计为：\(R[f] <= |\frac{(b-a)^3}{12n^2}M_2|\)，\(M_2\)为函数二阶导数绝对值最大值

2.2 复化Simpson公式

\[\int_a^bf(x)\mathrm{d}x \approx\sum_{k=1}^{n}\frac{h}{6}[f(x_{k-1}+4f(x_{k-\frac{1}{2}})+f(x_k))]=\frac{h}{6}[f(a)+4\sum_{k=1}^{n}f(x_{k-\frac{1}{2}})+2\sum_{k=1}^{n-1}f(x_k)+f(b)] \]

其误差估计为：\(R[f]\leq \frac{(b-a)^5M_4}{2880n^4}\)，\(M_4\)为函数四阶导数绝对值最大值

2.3 Gauss-Legendre公式

利用查表计算

注意：不要忘记将[a, b]上函数换元成[-1, 1]上再使用

误差估计为：\(R[f] \leq \frac{(b-a)^{2n+1}(n!)^4}{[(2n!)^3](2n+1)}M_{(2n)}\)，\(M_{(2n)}\)为函数四阶导数绝对值最大值

3. 代码实现

1. 在误差估计时候（尤其是Gauss型求积公式）需要计算被积函数的高阶导数，我们借助Sympy库的diff函数进行求导并用matplotlib.pyplot进行绘图操作
2. Gauss-Legendre并没有继承自IntegrationCal，当然略加更改就可以做到。

from math import sin, cos, sqrt, pi
import numpy as np
import matplotlib.pyplot as plt
# 汉字显示
plt.rcParams['font.sans-serif'] = ['SimHei']


# 数值积分接口描述
class IntegrationCal:
    def __init__(self, a, b, n):
        self.a = a
        self.b = b
        self.n = n

    def composite_cal(self, func):
        pass

    def error(self):
        pass


# Simpson公式类，内含普通Simpson与复化Simpson
class SimpsonCal(IntegrationCal):

    def composite_cal(self, func):
        h = (self.b - self.a) / (2 * self.n)
        f0 = func(self.a) + func(self.b)
        f1 = 0
        f2 = 0
        for i in range(2 * self.n - 1):
            x = self.a + i * h
            if i % 2 == 0:
                f2 = f2 + func(x)
            else:
                f1 = f1 + func(x)
        sn = h * (f0 + 2 * f2 + 4 * f1) / 3
        return sn

    # cos(x**2)的四阶导数，为单调递增的
    @staticmethod
    def func_4(x):
        return 8 * (4 * x ** 4 * cos(x ** 2) + 12 * x ** 2 * (x ** 2) - 3 * (x ** 2))

    def error(self):
        m = SimpsonCal.func_4(self.b)
        return m * (self.b - self.a) ** 5 / (2880 * self.n ** 4)


# 梯形公式类,内含普通梯形与复化梯形公式
class TrapeziumCal(IntegrationCal):

    def composite_cal(self, func):
        h = (self.b - self.a) / self.n
        sn = 0
        for i in range(self.n - 1):
            x0 = self.a + i * h
            x1 = self.a + (i + 1) * h
            sn += h * (func(x0) + func(x1)) / 2
        return sn

    def error(self):
        if self.b - sqrt(pi / 6) < 1e-3:
            m = 5.627598728468436
        elif self.b - pi / 4 < 1e-3:
            m = 7.271310062904556
        else:
            m = 7.689760105047184
        return m * (self.b - self.a) ** 3 / (12 * self.n ** 2)


# 只支持n = 5的Gauss-Legendre公式计算
class GaussLegendre:
    def __init__(self, theta):
        self.theta = theta
        self.n = 5

    def cal(self, func):
        v1 = 1 / 3 * sqrt(5 - 2 * sqrt(10 / 7))
        v2 = 1 / 3 * sqrt(5 + 2 * sqrt(10 / 7))
        w1 = (322 + 13 * sqrt(70)) / 900
        w2 = (322 - 13 * sqrt(70)) / 900
        return func(-v2, self.theta) * w2 + func(-v1, self.theta) * w1 + func(0, self.theta) * 128 / 225 + \
            func(v1, self.theta) * w1 + func(v2, self.theta) * w2

    def error(self):
        def factorial_cal(n):
            f = 1
            while n > 0:
                f *= n
                n -= 1
            return f
        if self.theta - pi/3 < 1e-3:
            m = 177660.6193610414
        elif self.theta - pi/4 < 1e-3:
            m = 177660.61853863852
        else:
            m = 532721.352671143
        return m * (sqrt(self.theta))**self.n * factorial_cal(self.n)**4 / 
            (factorial_cal(2*self.n)**3 * (2*self.n+1))


# 换元后函数
def a_func(x):
    return 2 * cos(x ** 2)


# 调整后Gauss-Legendre适用的函数
def func_for_gauss(x, theta):
    return sqrt(theta) * cos(((sqrt(theta) * (x + 1)) / 2) ** 2)


def main():
    # n = 5, theta变化
    # 参数选择与初始化
    theta_list = [sqrt(pi / 6), sqrt(pi / 4), sqrt(pi / 3)]
    trapezium_list = [TrapeziumCal(0, theta, 5) for theta in theta_list]
    simpson_list = [SimpsonCal(0, theta, 5) for theta in theta_list]
    gauss_list = [GaussLegendre(theta) for theta in [pi/6, pi/4, pi/3]]
    # 积分计算，采用列表生成式形式
    trapezium_result = [t.composite_cal(a_func) for t in trapezium_list]
    simpson_result = [s.composite_cal(a_func) for s in simpson_list]
    gauss_result = [g.cal(func_for_gauss) for g in gauss_list]
    # 误差分析
    trapezium_error = [t.error() for t in trapezium_list]
    simpson_error = [s.error() for s in simpson_list]
    gauss_error = [g.error() for g in gauss_list]
    print("梯形、simpson、gauss结果:", trapezium_result, simpson_result, gauss_result, sep='\n')
    print("梯形、simpson、gauss误差:", trapezium_error, simpson_error, gauss_error, sep='\n')

    # theta = pi / 4, n = 1:50
    n_list = list(range(1, 50))
    x_n = np.array(n_list)
    simpson_list2 = [SimpsonCal(0, pi / 4, n) for n in n_list]
    simpson_result2 = [s.composite_cal(a_func) for s in simpson_list2]
    simpson_error2 = np.array([s.error() for s in simpson_list2])
    plt.title("数值解与误差随n变化")
    plt.subplot(121)
    plt.plot(x_n, simpson_result2)
    plt.xlabel('n')
    plt.ylabel('数值解')
    plt.grid()
    plt.subplot(122)
    plt.plot(x_n, simpson_error2)
    plt.xlabel('n')
    plt.ylabel('误差估计')
    plt.grid()
    plt.tight_layout()
    plt.show()


def test():
    simpson = SimpsonCal(0, sqrt(pi / 3), 5)
    print(simpson.composite_cal(a_func))


if __name__ == '__main__':
    main()