复合高斯积分（节点数小于等于3的版本Python实现）

- 日理万妓 2022-05-23 08:47 92阅读 0赞

## 被积函数 ##

y=11\+x2  y = 1 1 + x 2

x∈\[−1,1\]  x ∈ \[ − 1 , 1 \]

## 算法描述 ##

**缺点：** 在下面代码中，没有实现点数大于3的情况（**等于3的话可以使用**）  
因为计算对应的系数值的时候，这里就不是直接用python的sympy似乎就搞不定了。  
**简述：** 就是通过勒让德（正交多项式）的零点。来做积分的考量点。之后再结合高斯积分的特点，得到对应的代数精度下方程，解出对应区间内的系数值。  
当然也可以直接做伸缩变换。之后就需要乘以一个系数。（在积分的微分变量求导，导出的一个系数）

这个思想类似于之前做的在契比雪夫多项式零点的拉格朗日插值。

这些思想告诉我们，其实在做逼近之类的事情的时候，我们如果可以发现在某些特别重要的特征，重点将这些特征逼近，就可以大大的提高效率。类似于信息论里的不定长编码，通过概率来导出最短信源编码长一样，来使得整体信源符号熵最小的过程类似。

都是强调对重要特征的重点描述。

## 代码 ##

import numpy as np
    from sympy import *
    
    
    def t(begin, end):
        xs = ps * (end - begin) / 2 + (begin + end) / 2
        temp = 0
        for i in range(point):
            temp += al[i] * Y.subs(x, xs[i])
        temp = temp.subs(ans)
        return temp.evalf() * (end - begin) / 2
    
    
    def P():
        if point == 0:
            return 1
        elif point == 1:
            return x
        p0 = 1
        p1 = x
        for i in range(point - 1):
            temp = ((2 * i + 3) * x * p1 - (i + 1) * p0) / (i + 2)
            p0 = p1
            p1 = temp
        return p1
    
    
    def loss(begin=-1, end=1):
        T = sum([t(xl[i], xl[i + 1]) for i in range(n)])
        I = integrate(Y, (x, begin, end))
        print('%.18f' % (I - T).evalf())
    
    
    if __name__ == '__main__':
        x = symbols('x')
        point = 2
        n = 10
        xl = np.linspace(-1, 1, n + 1)
        Y = 1 / (1 + x ** 2)
        ps = np.array(solve(P()))
    
        funcs = []
        al = symbols('a0:%d' % point)
    
        for i in range(point):
            if i != 0:
                temp = integrate((x ** i), (x, -1, 1))
                for j in range(point):
                    temp -= al[j] * (ps[j] ** i)
            else:
                temp = integrate(1, (x, -1, 1))
                for j in range(point):
                    temp -= al[j]
            funcs.append(temp)
        ans = solve(funcs)
    
        loss()