高斯勒让德(Gauss-legendre)求解多重积分(python,数值积分)

第四十四篇 高斯勒让德求解多重积分多重积分在工程分析中,经常需要在一个面积或体积上对函数进行积分。多重积分的解析方法在有限的情况下是可能的,但在这一篇中使用数值积分去求解。考虑两个变量的函数在二维区域R上的积分,如下图所示。函数f(x, y)可以被认为是在R区域上以直角从平面上伸出的第三维空间。需要的积分为通过前面篇章中描述的一维空间函数积分得到,涉及两个变量的二重积分将产生下面形式的积分法则样本点位于R的平面上,坐标为(xi, yi), i = 1,…N是函数的值。每个函数的计算都用wi
摘要由CSDN通过智能技术生成

第四十四篇 高斯勒让德求解多重积分

多重积分

在工程分析中,经常需要在一个面积或体积上对函数进行积分。多重积分的解析方法在有限的情况下是可能的,但在这一篇中使用数值积分去求解。一维的函数积分详见重复牛顿-科特斯积分重复高斯勒让德积分
考虑两个变量的函数在二维区域R上的积分,如下图所示。函数f(x, y)可以被认为是在R区域上以直角从平面上伸出的第三维空间。
在这里插入图片描述
需要的积分为
在这里插入图片描述
通过前面篇章中描述的一维空间函数积分得到,涉及两个变量的二重积分将产生下面形式的积分法则
在这里插入图片描述
样本点位于R的平面上,坐标为(xi, yi), i = 1,…N是函数的值。每个函数的计算都用wi加权,加权之和必须等于R的面积。
显然,在明确地定义不规则区域的积分计算起来有些困难。在实践中,将形状不规则的区域细分为若干更小的简单子区域(如矩形)就足够了,在这些子区域上可以很容易地进行数值积分。然后将每个子区域的解相加,得到整个区域R上的最终近似结果。这种方法类似于前面篇章中提到的“重复”规则。
首先,考虑在平行于笛卡尔坐标方向的矩形区域上的积分,因为这极大地简化了极限的定义。随后,这些概念被扩展到xy平面上一般的四边形和三角形区域的积分。本节中描述的所有方法都可以很容易地推广到三重积分。

矩形区域积分

在下图中,考虑函数f(x, y)在所示的矩形区域的积分。由于矩形的边界与笛卡尔坐标方向平行,变量之间没有关系,可以直接应用前面描述的积分方法。
在这里插入图片描述

牛顿柯特

例如,在每个方向上应用梯形规则将导致4个样本点(n = 4),如上图所示,矩形的每个角上都有一个样本点,从而产生该规则
在这里插入图片描述
在每个方向上应用Simpson规则,得到如下图6所示的9个样本点(n = 9),从而得到该法则为
在这里插入图片描述
在这里插入图片描述
其实可以猜到,上面得到的的加权系数来自于在每一个方向分别考虑的加权系数的简单乘积

计算实例

在每个方向使用梯形法则去计算下面积分
在这里插入图片描述
在一维中,梯形规则有2个样本点,所以在二重积分中,我们将使用n = 2*2 = 4个样本点。由上图知,h = 1, k = 2,有梯形法则的二维公式得到
在这里插入图片描述
精确解为15.3333

计算实例

在每个方向使用辛普森法则去计算下面的积分
在这里插入图片描述
在一维中,辛普森法则有3个样本点,所以在二重积分中,使用n = 3**2 = 9个采样点。由上图得,h = 0.5, k = 1,根据三点辛普森法则得到,二维积分得
在这里插入图片描述
在这个例子中,在两个方向上都使用辛普森法是低效的。在x方向上用辛普森法则和在y方向上用梯形法则也可以得到精确的解。可以使用n = 3 × 2=6个样本点,得到
#

高斯-勒让德

高斯-勒让德法则也可以应用于这种类型的多重积分,但必须注意找到样本点的正确位置。一种方法是进行坐标变换,使每个方向上的积分极限变为±1。就可以使用这篇中提到得权值和采样点值。详细的表格对应数值可见高斯勒让德求积分
矩形区域的另一种方法是将样本点在每个坐标方向上的范围中点左右进行加权和比例系数相乘。
考虑下图所示的矩形区域上的两点高斯-勒让德积分。
在这里插入图片描述
在x方向上,样本点坐标为
在这里插入图片描述
对应的加权系数w1 = w2 = h。同样在y方向上,样本点为
在这里插入图片描述
对应加权系数w1 = w2 = k,则法则为
在这里插入图片描述

计算实例

在两个方向使用高斯-勒让德计算下面的积分
在这里插入图片描述
在这个二重积分中,我们使用n = 2*22 = 4个样本点。由上图可知,h = 1, k = 0.5,根据上面得式子,样本点位于
在这里插入图片描述
因此得到
在这里插入图片描述
精确解为0.8
在这个例子中,两点法则能够在x方向上精确积分,但不能在y方向上。
通过类似的推理,由图下图所示的三点高斯勒让德法则可以得到x中得样本点和加权系数
x方向为
在这里插入图片描述
在y方向为
在这里插入图片描述
因此法则将变成
在这里插入图片描述
在这里插入图片描述

计算实例

在每个方向使用3点高斯-勒让德去计算下面的积分
在这里插入图片描述
在这个二重积分中,我们将使用n = 3**2 = 9个采样点。由上图可知,h = 1, k = 0.5,然后根据样本点坐标值得方程求得样本点为
在这里插入图片描述
因此根据积分方程,得到I = 0。8,这是精确解。
在这种情况下,四次项在y上的积分的精确解可以通过三点法则得到,在x方向上使用两点法则,总共有6个样本点。尽管在不同的方向上实现不同的数量积分法则是很容易的,但是大多数数值积分软件包在所有方向上都使用相同的数量法则。

程序如下

其中有一个主程序,和两个子程序,分别为高斯勒让德的样本点值和权值的子程序gauss_legendre;和局部坐标的形状函数的子程序fun_der。

#高斯拉盖尔法则
import numpy as np
import B
import math
ndim=2;nsp=9
if ndim==1:
    nod=2
elif ndim==2:
    nod=4
elif ndim==3:
    nod=8
coord=np.zeros((nod,ndim));der=np.zeros((ndim,nod));fun=np.zeros(nod)
coord[:,0]=(0.0,1.0,4.0,6.0);coord[:,1]=(0.0,2.0,3.0,1.0)
samp=np.zeros((nsp,ndim));wt=np.zeros((nsp))
B.gauss_legendre(samp,wt)
res=0
def f65(point):
    x=point[0];y=point[1]#;z=point[2]
    f65=x**2*y**2
    return f65
for i in range(1,nsp+1):
    B.fun_der(fun,der,samp,i)
    res=res+np.linalg.det(np.dot(der,coord))*wt[i-1]*f65(np.dot(np.transpose(coord),fun))
print('高斯-勒让德法则的多重积分')
print('维数',ndim)
for i in range(1,nod+1):
    print('坐标 (x,y[,z])',coord[i-1,:])
print('样本点的数量',nsp)
print('计算结果   ','{:13.4e}'.format(res))
gauss-legendre
def gauss_legendre(samp,wt):
    nsp=samp.shape[0]
    ndim=samp.shape[1]
    if ndim==1:
        if nsp==1:
            samp[0,0]=  0.0 
            wt[0]  =    2.0 
        elif nsp==2:
            samp[0,0]= -0.57735026918962576449 
            samp[1,0]=  0.57735026918962576449 
            wt[0]=      1.0 
            wt[1]=      1.0 
        elif nsp==3:
            samp[0,0]= -0.77459666924148337704 
            samp[1,0]=  0.0 
            samp[2,0]=  0.77459666924148337704 
            wt[0]=      0.55555555555555555556 
            wt[1]=      0.88888888888888888889 
            wt[2]=      0.55555555555555555556 
        elif nsp==4:
            samp[0,0]= -0.86113631159405257524 
            samp[1,0]= -0.33998104358485626481 
            samp[2,0]=  0.33998104358485626481 
            samp[3,0]=  0.86113631159405257524 
            wt[0]=      0.34785484513745385737 
            wt[1]=      0.65214515486254614271 
            wt[2]=      0.65214515486254614271 
            wt[3]=      0.34785484513745385737 
        elif nsp==5:
            samp[0,0]= -0.90617984593866399282 
            samp[1,0]= -0.53846931010568309105 
            samp[2,0]=  0.0 
            samp[3,0]=  0.53846931010568309105 
            samp[4,0]=  0.90617984593866399282 
            wt[0]=      0.23692688505618908751 
            wt[1]=      0.47862867049936646804 
            wt[2]=      0.56888888888888888889 
            wt[3]=      0.47862867049936646804 
            wt[4]=      0.23692688505618908751    
        else:
            print('积分点数量错误')
    elif ndim==2:
        if nsp==1:
            samp[0,0]=  0.0  
            samp[0,1]=  0.0 
            wt[0]=      4.0 
        elif nsp==4:
            samp[0,0]= -0.57735026918962576449 
            samp[1,0]=  0.57735026918962576449 
            samp[2,0]= -0.57735026918962576449 
            samp[3,0]=  0.57735026918962576449 
            samp[0,1
  • 10
    点赞
  • 52
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 8
    评论
评论 8
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

深渊潜航

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值