实验内容

试建立8个节点的Gauss-Legendre求积公式，要求数值求积公式中的Gauss点以及求积系数都具有8位有效数字；用所建立的Gauss型求积公式计算积分

的近似值，并用四等分区间的复化两点Gauss型求积公式计算该积分的近似值，对二者的精度进行比较，并解释其原因。（提示：用非线性方程迭代解法求得Gauss 点，求积系数根据插值型求积公式的系数计算，具体的积分值可以用数值方法求得)。

算法原理

高斯——勒让德求积公式

插值型求积公式任意给定n + 1个节点，至少有 n 次代数精度。 牛顿-柯特斯公式是等距节点的插值求积公式，当 n 为 偶数时，有n + 1次代数精度。插值求积公式中有n + 1个节点，n + 1个求积系数，共有2n + 2个未知参数，适当选取这些参数可使求积公式具有2n + 1次代数精度，这时就构成了高斯型求积公式。高斯型求积公式就是最高次代数精度的求积公式。高斯公式不局限与等步长分区域，意思是不用二分的方法直接划分区域，在简单的计算下就可以得到高精度的积分式子。插值型求积公式的节点a < x0 < x1 < … < xn < b是高斯点的充分必要条件是以这些节点为零点的多项式。与任何次数不超过n的多项式P（x）带权正交，同时由于其系数都大于0，所以我们也可以得到该求积公式是收敛的。对于勒让德求积公式，就是说我们取权函数P（x）=1，同时区间为[-1, 1]，所以勒让德多项式的零点就是我们的求积的高斯点。使插值求积公式有2n+1次代数精度的节点xi,i=0,1,…,n称为高斯点，该插值求积公式称为高斯型求积公式，简称高斯公式。高斯点的确定原则上可以根据它们满足的代数精度要求、通过求解非线性方程组得出。对勒让德多项式取其零点作为求积节点，即可构造出高斯-勒让德(Gauss-Legendre)公式，也简称高斯公式。

复化高斯公式

将求积区间[a，b]分为 n 等分：n=(b-a)/h上使用低阶高斯公式(如一阶中矩形公式):
复化两点高斯型求积公式
高斯公式的值随节点数目的增加而收敛到准确积分值。

优点:

  • 收敛、稳定；

  • 计算量小，代数精度高。

缺点：

  • 高斯点难求（即多项式的根难求）；

  • 高斯点是无理数，求积系数也是无理数。

使用情况：

  • f(x) 赋值量大；

  • 计算的积分多。

程序框图

高斯——勒让德求积公式
复化高斯型求积公式

编程实现

本次编程采用Python的面向对象方法进行编写，其中对算法的初步理解以及求解方程的根的数值方法是本程序的难点所在。
利用所编程序，确定并计算题目中所给积分的数值解:
以下将用Python进行编程实验，分别利用高斯——勒让德求积公式与复化两点高斯型求积公式对积分要求积分进行求解，最后进行结果分析与讨论。
Python代码详解及运行结果如下：

代码详解：

gauss_legendre_int.py: 高斯——勒让德求积公式文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# -*- coding: utf-8 -*-
# @Time    : 2022/5/13 4:05 下午
# @Author  : Leo
# @FileName: gauss_legendre_int.py
# @Software: PyCharm
# @Blog    ：https://guojxblog.cn
# @GitHub  ：https://github.com/guojx0820
# @Email   ：[email protected]

import numpy as np
import sympy
import math


class GaussLegendreIntegration:
    """
    高斯——勒让德求积公式，核心部分求解零点与系数
    """

    def __init__(self, int_fun, int_interval, zeros_num=10):
        self.int_fun = int_fun  # 被积函数，符号定义
        if len(int_interval) == 2:
            self.a, self.b = int_interval[0], int_interval[1]  # 被积区间
        else:
            raise ValueError("积分区间参数设置不规范，应为[a,b].")
        self.n = int(zeros_num)  # 正交多项式的零点数
        self.zeros_points = None  # 勒让德高斯零点
        self.int_value = None  # 积分值结果
        self.A_k = None  # 求积系数

    def cal_int(self):
        """
        高斯——勒让德求积公式
        :return:
        """
        self._cal_Ak_coef_()  # 求解求积系数Ak与零点
        fun_val = self.int_fun(self.zeros_points)  # 零点函数值
        self.int_value = np.dot(self.A_k, fun_val)  # 插值型求积公式
        return self.int_value

    def _cal_gauss_zeros_points_(self):
        """
        计算高斯零点
        :return:
        """
        t = sympy.Symbol('t')
        # 勒让德多项式构造
        p_n = (t ** 2 - 1) ** self.n / (math.factorial(self.n) * 2 ** self.n)
        diff_p_n = sympy.diff(p_n, t, self.n)  # 多项式的n阶导数
        # print(diff_p_n)
        # 求解高斯——勒让德多项式的全部零点
        self.zeros_points = np.asarray(sympy.solve(diff_p_n, t), dtype=np.float_)
        # print("高斯点:", self.zeros_points, sep="\n")
        return diff_p_n, t

    def _cal_Ak_coef_(self):
        """
        计算Ak系数
        :return:
        """
        diff_p_n, t = self._cal_gauss_zeros_points_()  # 求解高斯零点，传递函数和符号参数
        Ak_poly = sympy.lambdify(t, 2 / ((1 - t ** 2) * (diff_p_n.diff(t, 1) ** 2)))
        self.A_k = Ak_poly(self.zeros_points)  # 求解求积系数Ak
        # 区间转换，[a, b] --> [-1, 1]
        self.A_k = self.A_k * (self.b - self.a) / 2
        self.zeros_points = (self.b - self.a) / 2 * self.zeros_points + (self.a + self.b) / 2
        # print("高斯点:", self.zeros_points, sep="\n")
        # print("求积系数:", self.A_k, sep="\n")

test_gauss_legendre.py: 高斯——勒让德求积公式测试文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# -*- coding: utf-8 -*-
# @Time    : 2022/5/13 4:57 下午
# @Author  : Leo
# @FileName: test_gauss_legendre.py
# @Software: PyCharm
# @Blog    ：https://guojxblog.cn
# @GitHub  ：https://github.com/guojx0820
# @Email   ：[email protected]

import numpy as np
import sympy
import matplotlib.pyplot as plt
from NumericalIntegration.gauss_legendre_int import GaussLegendreIntegration


def fun(x):
    """
    定义积分的被积函数
    :param x:自变量
    :return:被积函数
    """
    return np.cos(x) * np.exp(x)


if __name__ == "__main__":
    gauss_zeros_num = np.arange(8, 20, 1)
    int_accurate = -0.5 * (np.exp(np.pi) + 1)
    print("精确值:", int_accurate)
    precision = []
    for num in gauss_zeros_num:
        legendre = GaussLegendreIntegration(fun, [0, np.pi], zeros_num=num)
        int_value = legendre.cal_int()
        precision.append(int_accurate - int_value)
        print("num:%d,积分值:%.15f,误差:%.15e" % (num, int_value, precision[-1]))  # append list中最后一个值用precision[-1]
    plt.figure(figsize=(8, 6))
    plt.title("Gauss-Legendre Quadrature Formula Error Curve", fontdict={"fontsize": 15})
    plt.xlabel("The Number of Gauss Zero-Points", fontdict={"fontsize": 12})
    plt.ylabel("Error", fontdict={"fontsize": 12})
    plt.plot(gauss_zeros_num, precision, 'ro-')
    plt.grid(axis='y', color='b', linestyle='--', linewidth=0.5)
    plt.savefig("/Users/leo/Desktop/Gauss-Legendre.png", dpi=300, bbox_inches="tight")
    plt.show()
    # legendre._cal_gauss_zeros_points_()

composite_2points_gauss_legendre.py:复化两点高斯型求积公式文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# -*- coding: utf-8 -*-
# @Time    : 2022/5/15 5:11 下午
# @Author  : Leo
# @FileName: composite_2points_gauss_legendre.py
# @Software: PyCharm
# @Blog    ：https://guojxblog.cn
# @GitHub  ：https://github.com/guojx0820
# @Email   ：[email protected]

class CompositeGaussLegendreIntegration:
    """
    复化两点高斯——勒让德求积公式，等距节点数自定义
    """

    def __init__(self, int_fun, int_interval, interval_num=4):
        self.int_fun = int_fun  # 被积函数，符号定义
        if len(int_interval) == 2:
            self.a, self.b = int_interval[0], int_interval[1]  # 被积区间
        else:
            raise ValueError("积分区间参数设置不规范，应为[a, b].")
        self.interval_num = interval_num
        self.int_value = None

    def cal_int(self):
        """
        复化两点高斯——勒让德求积公式
        :return:
        """
        fun_value = 0  # 初始化函数值
        interval_len = (self.b - self.a) / self.interval_num
        for k in range(self.interval_num - 1):
            fun_value += self.int_fun(self.a + (k + 1 / 2) * interval_len)
        self.int_value = interval_len * fun_value
        return self.int_value

test_composite_2points_gauss_legendre.py:复化两点高斯型求积公式测试文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
# -*- coding: utf-8 -*-
# @Time    : 2022/5/15 5:56 下午
# @Author  : Leo
# @FileName: test_composite_2points_gauss_legendre.py
# @Software: PyCharm
# @Blog    ：https://guojxblog.cn
# @GitHub  ：https://github.com/guojx0820
# @Email   ：[email protected]

import numpy as np
import matplotlib.pyplot as plt
import symbol
from NumericalIntegration.composite_2points_gauss_legendre import CompositeGaussLegendreIntegration


def fun(x):
    """
    定义积分被积函数
    :param x: 自由变量
    :return:
    """
    return np.cos(x) * np.exp(x)


if __name__ == "__main__":
    interval_num = np.arange(100000, 1000000, 100000)
    int_accurate = -0.5 * (np.exp(np.pi) + 1)
    print("精确解:", int_accurate)
    precision = []
    for num in interval_num:
        legendre = CompositeGaussLegendreIntegration(fun, [0, np.pi], interval_num=num)
        int_value = legendre.cal_int()
        precision.append(int_accurate - int_value)
        print("num=%d,积分值:%.15f,误差:%.15e" % (num, int_value, precision[-1]))  # append list中最后一个值用precision[-1]
    plt.figure(figsize=(8, 6))
    plt.title("Composite Gauss Quadrature Formula Error Curve", fontdict={"fontsize": 16})
    plt.xlabel("The Number of Integral Intervals", fontdict={"fontsize": 12})
    plt.ylabel("Error", fontdict={"fontsize": 12})
    plt.plot(interval_num, precision, 'ro-')
    plt.grid(axis='y', color='b', linestyle='--', linewidth=0.5)
    plt.savefig("/Users/leo/Desktop/CompositeGauss.png", dpi=300, bbox_inches="tight")
    plt.show()

运行结果

高斯——勒让德求积公式高斯点计算结果

高斯——勒让德求积公式求积系数计算结果

高斯——勒让德求积公式精确值与数值解及误差

高斯——勒让德求积公式误差曲线

复化两点高斯求积公式精确值与数值解及误差

复化高斯型求积公式误差曲线

结果分析与讨论

  • 有两种求积公式的结果可以看出，高斯——勒让德求积公式计算较为复杂，但收敛速度快且精度较高，在20个高斯零点以内可以收敛，且误差缩小到1e-14的个数量级；而复化两点高斯型积分公式收敛速度较慢，需要将区间分为1e6个小区间之后，精度只能达到1e-6个数量级，继续增加其区间，缩小步长时，程序运行速度明显变慢，说明其效率低于高斯——勒让德求积公式（由于n=4，四等分区间时，其结果差异很大，所以舍弃，只能用百万区间划分才能确保精度）。

  • 高斯——勒让德求积公式精度高，收敛速度快的主要原因是：高斯——勒让德求积公式的高斯节点的选取的特殊性，并不采用等距节点，而是利用勒让德多项式的插值节点作为高斯点进行积分，其结果更加精确，且收敛速度更快，也更加容易接近精确解。

  • 反观等距节点插值，对于高次具有强震荡曲线的函数，效果并不理想。主要是由于函数在小区间内的变化率太大，所引起的误差较大，因此需要不断增加其等距节点个数，缩小步长或区间范围，在小区间中用高斯型求积公式进行积分求解才更有效果。

参考

[1]《数值分析》欧阳洁等。

[2] https://www.bilibili.com/video/BV1Lq4y1U7Hj?p=29&spm_id_from=pageDriver

[3] https://wenku.baidu.com/view/b3fc62d74b7302768e9951e79b89680203d86bb7.html

[4] https://blog.csdn.net/poqw2010/article/details/79487919

[5] https://blog.csdn.net/catinlbb/article/details/103534439

[6] https://www.doc88.com/p-3901680450501.html

[7] https://zhuanlan.zhihu.com/p/439768015