高斯积分和代码实现

RangoLei_Lzs

已于 2024-02-26 19:52:44 修改

阅读量687

点赞数 19

分类专栏：力学文章标签： python 算法

于 2024-02-26 19:50:21 首次发布

本文链接：https://blog.csdn.net/leizhengshenglzs/article/details/136306644

版权

力学专栏收录该内容

2 篇文章 0 订阅

订阅专栏

高斯积分理论

高斯积分是有限元中常用的数值积分方法，由于有限元插值函数一般为多项式，被积函数一般也是多项式，高斯积分使用很少的积分点，就能达到较高的积分精度 Legendre-Gauss Quadrature。

高斯积分使用N个积分点就能达到2N-1阶精度，为了定义统一的积分格式，我们一般将积分限[a,b]变换到[-1,1]，这样选定n，wi，xi后就能计算任意积分限的积分。

高斯积分在区间[-1,1]上的积分点可由由n阶勒让德多多项式的根计算 Legendre polynomial。

权值由如下公式计算

高级积分点和权值计算

计算一维积分高斯点和权重~

def gauss_1d(npts):
    """Return Gauss points and weights for Gauss quadrature in 1D

    Parameters
    ----------
    npts : int
      Number of quadrature points.

    Returns
    -------
    wts : ndarray
      Weights for the Gauss-Legendre quadrature.
    pts : ndarray
      Points for the Gauss-Legendre quadrature.
    """
    if npts == 2:
        pts = [-0.577350269189625764, 0.577350269189625764]
        wts = [1.00000000000000000, 1.00000000000000000]
    elif npts == 3:
        pts = [-0.774596669241483377, 0, 0.774596669241483377]
        wts = [0.555555555555555556, 0.888888888888888889,
               0.555555555555555556]
    elif npts == 4:
        pts = [-0.861136311594052575, -0.339981043584856265,
               0.339981043584856265, 0.861136311594052575]
        wts = [0.347854845137453857, 0.652145154862546143,
               0.652145154862546143, 0.347854845137453857]  
    else:
        msg = "The number of points should be in [2, 10]"
        raise ValueError(msg)

    return pts, wts

多维高斯积分点和权值

多维高斯积分点是一维的张量积 , 例如四边形，每个方向两个积分点P[2] = {-0.577350269189625764, 0.577350269189625764}，则 $P_{ij} = (P_i,\ P_j)$ , i,j下标为积分点索引，积分点的坐标值为相应一维值的组合。

def gauss_nd(npts, ndim=2):
    """
    Return Gauss points and weights for Gauss quadrature in
    an ND hypercube.

    Parameters
    ----------
    npts : int
      Number of quadrature points.

    Returns
    -------
    nd_wts : ndarray
      Weights for the Gauss-Legendre quadrature.
    nd_pts : ndarray
      Points for the Gauss-Legendre quadrature.
    """
    pts, wts = gauss_1d(npts)
    nd_pts = np.array(list(product(pts, repeat=ndim)))
    nd_wts = product(wts, repeat=ndim)
    nd_wts = [np.prod(nd_wt) for nd_wt in nd_wts]
    return nd_pts, nd_wts

上面代码中，函数product为张量积，得到的结果为积分点数组，积分点个数为 total_pts = pow(npts,ndim)，数组维数为[total_pts，ndim] ，可参考如下计算示例：

//2个积分点
pts, wts = gauss_1d(2)
print(pts,wts)
#[-0.5773502691896257, 0.5773502691896257] [1.0, 1.0]

//三维积分，三重张量积
nd_pts = np.array(list(product(pts, repeat=3)))
print(nd_pts)
# result
'''
[[-0.57735027 -0.57735027 -0.57735027]
 [-0.57735027 -0.57735027  0.57735027]
 [-0.57735027  0.57735027 -0.57735027]
 [-0.57735027  0.57735027  0.57735027]
 [ 0.57735027 -0.57735027 -0.57735027]
 [ 0.57735027 -0.57735027  0.57735027]
 [ 0.57735027  0.57735027 -0.57735027]
 [ 0.57735027  0.57735027  0.57735027]]
 '''
#
nd_wts = list(product(wts, repeat=3))
print(nd_wts)
#[(1.0, 1.0, 1.0), (1.0, 1.0, 1.0), (1.0, 1.0, 1.0), (1.0, 1.0, 1.0), 
#(1.0, 1.0, 1.0), (1.0, 1.0, 1.0), (1.0, 1.0, 1.0), (1.0, 1.0, 1.0)]

//最后一行将组合的积分点的weight相乘，得到积分点的weight
nd_wts = [np.prod(nd_wt) for nd_wt in nd_wts]
#[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]