首页 > 编程语言 >【数值计算方法】2&3维高斯积分的python实现

【数值计算方法】2&3维高斯积分的python实现

时间:2024-08-05 17:29:18浏览次数:14  
标签:upLimit 高斯 python 积分 float lowLimit 计算方法 lambda

目录

本文只给出pythont实现和例题,数学推导见【数值计算方法】数值积分&微分-python实现 - FE-有限元鹰 - 博客园

二维高斯积分

python实现二维高斯积分:


def Integ2dGuassLegendre(f,lowLimit:List[float]=[-1,-1], 
                            upLimit:List[float]=[1,1], 
                            n:int=3)->float:
    """给定积分区域[lowLimit,upLimit]和高斯积分点个数n(>=1),计算二维高斯-勒让德积分公式"""
    a,b,c,d=lowLimit[0],upLimit[0],lowLimit[1],upLimit[1]
    if n<=0:
        raise ValueError("高斯-勒让德积分时,n必须大于0")
    if n==1:
        return 4*f(0,0)
    if a==-1 and b==1 and c==-1 and d==1:
        # 标准型积分
        #计算权重和积分点位置
        x_is,w_is=legendre_gauss_points_and_weights(n)
        y_js,w_js=legendre_gauss_points_and_weights(n)

        return np.sum([w_is[ind_x]*w_js[ind_y]*f(xi,yj) for ind_x,xi in enumerate(x_is) for ind_y,yj in enumerate(y_js) ])
    else:
        # 非标准型积分,积分区域:[a,b]x[c,d]
        xt1=lambda t1: 0.5*(b-a)*t1+0.5*(b+a)
        yt2=lambda t2: 0.5*(d-c)*t2+0.5*(d+c)
        
        t1_is,w_is=legendre_gauss_points_and_weights(n)
        t2_js,w_js=legendre_gauss_points_and_weights(n)
        
        return np.sum([w_is[ind_x]*w_js[ind_y]*f(xt1(t1i),yt2(t2j))*(b-a)*(d-c)*0.25 for ind_x,t1i in enumerate(t1_is) for ind_y,t2j in enumerate(t2_js) ])

三维高斯积分

python实现:

def Integ3dGuassLegendre(f,lowLimit:List[float]=[-1,-1,-1], 
                            upLimit:List[float]=[1,1,1], 
                            n:int=3)->float:
    """给定积分区域[lowLimit,upLimit]和高斯积分点个数n(>=1),计算二维高斯-勒让德积分公式"""
    a,b=lowLimit[0],upLimit[0]
    c,d=lowLimit[1],upLimit[1]
    g,h=lowLimit[2],upLimit[2]
    if n<=0:
        raise ValueError("高斯-勒让德积分时,n必须大于0")
    if n==1:
        return 8*f(0,0)
    if a==-1 and b==1 and c==-1 and d==1:
        # 标准型积分
        #计算权重和积分点位置
        x_is,w_is=legendre_gauss_points_and_weights(n)
        y_js,w_js=legendre_gauss_points_and_weights(n)
        z_js,w_ks=legendre_gauss_points_and_weights(n)

        return np.sum([w_is[ind_x]*w_js[ind_y]*w_ks[ind_z]*f(xi,yj,zk) for ind_x,xi in enumerate(x_is) for ind_y,yj in enumerate(y_js) for ind_z,zk in enumerate(z_js)])
    else:
        # 非标准型积分,积分区域:[a,b]x[c,d]x[g,h]
        xt1=lambda t1: 0.5*(b-a)*t1+0.5*(b+a)
        yt2=lambda t2: 0.5*(d-c)*t2+0.5*(d+c)
        zt3=lambda t3: 0.5*(h-g)*t3+0.5*(h+g)
        
        t1_is,w_is=legendre_gauss_points_and_weights(n)
        t2_js,w_js=legendre_gauss_points_and_weights(n)
        t3_ks,w_ks=legendre_gauss_points_and_weights(n)
        
        return np.sum([w_is[ind_x]*w_js[ind_y]*w_ks[ind_z]*f(xt1(t1i),yt2(t2j),zt3(t3k))*(b-a)*(d-c)*(h-g)*0.125 for ind_x,t1i in enumerate(t1_is) for ind_y,t2j in enumerate(t2_js) for ind_z,t3k in enumerate(t3_ks)])

验证

from formu_lib import *
import numpy as np
from scipy.integrate import dblquad
# 定义被积函数
def integrand(x, y):
    return np.exp(x*x+y*y)

# 计算二重积分
result, error = dblquad(integrand, -1, 1, lambda x: -1, lambda x: 1)
print("numpy 二重积分结果:", result)
ans1=Integ2dGuassLegendre(integrand,[-1, -1],[1, 1],n=5)
print(f"本地二重积分结果:{ans1}")

from scipy.integrate import tplquad
# 定义被积函数
def integrand3(x, y, z):
    return np.exp(x*x+y*y+z*z)

# 计算三重积分
result3, error = tplquad(integrand3, -1, 1, lambda x: -1, lambda x: 1, lambda x, y: -1, lambda x, y: 1)
ans2=Integ3dGuassLegendre(integrand3,[-1,-1,-1],[1,1,1],n=5)

print("numpy三重积分结果:", result3)
print(f"本地三重积分结果:{ans2}")

输出:

  • numpy 二重积分结果: 8.557400519221307
  • 本地二重积分结果:8.557173227239266
  • numpy三重积分结果: 25.03299361973213
  • 本地三重积分结果:25.03199627931168

标签:upLimit,高斯,python,积分,float,lowLimit,计算方法,lambda
From: https://www.cnblogs.com/aksoam/p/18343674

相关文章

  • Python选择与循环
    条件语句Python中,选择(条件)语句可细分为3种形式,分别是if语句、ifelse语句和ifelifelse语句。标准格式:if表达式1:代码块1elif表达式2:代码块2elif表达式3:代码块3...//可以有零条或多条elif语句else:代码块n表达式可以是一个单纯的布尔值......
  • 【Python&GIS】Arcpy中常用出图函数详解
        出图是每个GISer都要经历的事,但有时候会有许多重复且多且无聊的出图任务,这个时候我们肯定想能不能自动化出图。ArcGIS中的模型创建就可以实现,但是我的数据大部分是在Python中处理的,所以就想能不能使用Python进行批量出图,正好今天跟大家分享一下。这里使用的mxd作为......
  • 使用Python 和 Selenium 抓取 酷狗 音乐专辑 附源码
    在这篇博客中,我将分享如何使用Python和Selenium抓取酷狗音乐网站上的歌曲信息。我们将使用BeautifulSoup解析HTML内容,并提取歌曲和专辑信息。依赖库requestsbeautifulsoup4selenium准备工作首先,我们需要安装一些必要的库:pipinstallrequestsbeautifulsoup4selenium......
  • 使用 Python和PyQt5 打造 你的专属文件查询工具! 附源码
    本文将介绍如何使用Python和PyQt5创建一个简单的文件查询工具。该工具允许用户选择一个目录,并在该目录中搜索特定的文件。依赖库首先,确保你已经安装了PyQt5库:pipinstallPyQt5步骤第一步:导入库我们需要导入必要的库,包括sys、os和PyQt5。importsysimportosfromP......
  • 用Python和PyQt5打造你的专属音乐播放器!轻松创建带封面的音乐列表
    在本文中,我们将介绍如何使用Python的PyQt5库创建一个简单的音乐播放器。这个音乐播放器可以显示歌曲的封面,并且点击封面就可以播放对应的歌曲。依赖库首先,我们需要安装PyQt5库。可以使用以下命令进行安装:pipinstallPyQt5功能如下:显示歌曲列表:程序会在指定的目录(在......
  • 了解 Databricks 文件系统 (DBFS) 中的文件访问与使用 Python 和 Spark 的卷的比较
    我当前正在尝试从Databricks文件系统(DBFS)读取和显示文件,但遇到了问题。这是我使用的代码:file_path="/dbfs/cluster-logs/use_case/default_job_cluster/cluster_id/init_scripts/cluster_id/20240801_proxy-init.sh.stderr.log"withopen(file_path,'r')asfile:......
  • 机器学习用python还是R,哪个更好?
    选择使用Python还是R进行机器学习取决于多个因素,包括您的具体需求、项目要求、个人偏好以及团队的技能水平。以下是一些关键点,可以帮助您做出决定:Python的优势广泛使用:Python是目前最流行的编程语言之一,特别是在数据科学和机器学习领域。它有一个庞大的社区和丰富的资源。......
  • python 元类:在调用“__set_name__”方法后编辑命名空间?
    假设我们用元类定义一个类。在类主体中,分配了对象,这些对象实现__set_name__以在类的数据结构中注册自身。是否可以在运行方法之后编辑命名空间?比如,分离填充的数据结构,将其分成两部分,然后在新属性下添加部分?__set_name__问题在于,在元类中调用之......
  • 无法加入进程,只能终止[Python 3.11,多处理]
    我有一个问题要问对Python的多处理库有更​​多经验的人,此时我几乎迷失了方向。我目前正在构建一个应该在Windows11和Windows11上运行的图像处理应用程序装有DebianLinux的OrangePi5。我的设置是,除了主程序之外,还有另外两个进程,一个用于处理不间断的......
  • Python 将Word转换为JPG、PNG、SVG图片
    将Word文档以图片形式导出,既能方便信息的分享,也能保护数据安全,避免被二次编辑。文本将介绍如何使用 Spire.DocforPython 库在Python程序中实现Word到图片的批量转换。Python将Word转换为JPG、JPEG、PNG、BMP等图片格式Python将Word文档转换为SVG格式 Python库安装: ......