【Scipy高级计算】(1) 定积分、二重积分，附python完整代码

2024-06-13 15:28| 来源: 网络整理| 查看: 265

大家好，偶然接触到了python的Scipy高级科学计算库，可以看作是numpy的高级版，它可以很方便的求解线性代数、积分、插值、傅里叶变换等计算。今天和大家分享一下Scipy库的积分计算。

''' 安装： pip install scipy 主要模块 # ---------------------------------------------------------- # scipy.cluster 层次聚类模块，包含适量优化、k-means等 scipy.constants 数学科学常量 pi, e scipy.fft 快速傅里叶变换 scipy.integrate 积分模块，求多重积分、高斯积分、解常微分方程 scipy.interpolate 插值模块，提供样条插值、径向基函数插值 scipy.io 输入输出，文件读取 scipy.linalg 线性代数 scipy.misc 图像处理 scipy.optimize 优化算法，求解有约束或无约束的优化 scipy.signal 信号处理，卷积、差分等滤波方法和各种谱方法 scipy.stats 统计函数 # ---------------------------------------------------------- # ''' 1. scipy.integrate 库介绍

scipy.integrate 提供数值积分和常微分方程的求解算法，具体引用函数如下

''' 函数及其功能 # ---------------------------------------------------------- # quad 一重积分求解方法 dblquad 二重积分求解方法 tplquad 三重积分求解方法 nquad N重积分求解方法 fixed_quad 对func(x)做n维高斯积分 quadrature 求给定容限范围内的高斯积分 romberg 对函数做Romberg积分 # ---------------------------------------------------------- # ''' 2. 定积分计算

定积分的计算方法使用 scipy.integrate 模块下的 quad 函数。函数的参数和返回值如下：

''' 参数 # ---------------------------------------------------------- # func 被积函数，自定义函数 a 积分下限 b 积分上限 args 元组。func被积函数的未知变量X以外，需要的其他变量的参数， # ---------------------------------------------------------- # 返回值： (积分结果, 误差值) 误差用来辅助判断这个积分结果准不准 ''' 案例一：计算 $x^{2}$ 在 [0, 2] 积分区间的结果

计算公式为： $\int_{0}^{2}x^{2} dx$

被积函数中只有一个自变量x，积分函数指定参数args，给出非自变量的值，元组类型

函数返回结果是一个元组类型，包含两个元素，一个是积分计算结果，另一个是积分计算误差，一般只需要看计算结果就行。

from scipy.integrate import quad # 定义函数x的n次方，自变量x def myfunction(x, n): return x**n # 计算函数对x的积分，下限0上限2，被积函数中的n=2 data = quad(func=myfunction, a=0, b=2, args=(2,)) print(data) # 积分结果2.666, 计算误差2.9605e-14 # (2.666666666666667, 2.960594732333751e-14) 案例二：计算 $k\times x^{n}$ 在 [0, 2] 积分区间的结果

和案例一类似，被积函数中只有一个自变量x，而k和n在积分函数中指定，参数args指定(k, n)

下面计算 $3\times x^{2}$ 在[0,2]上下限内的积分结果

from scipy.integrate import quad # 定义函数k乘x的n次方，自变量x def myfunction(x, n, k): return k * x**n # 计算函数对x的积分，下限0上限2，被积函数中的n=2,k=3 data = quad(func=myfunction, a=0, b=2, args=(2,3)) print(data) # 积分结果8.0, 误差8.881e-14 # (8.0, 8.881784197001252e-14) 案例三：不同阶数对定积分结果的影响

现有被积函数 $k\times x^{n}$ ，积分上下限[0, 1]，若k=3，研究不同的阶数n对定积分结果的影响，n=0,1,2，...。如下，我们使用一个循环，比较 $3\times x^{0}$ 到 $3\times x^{14}$ 的积分结果，将每次积分结果保存在result列表中，绘制曲线。

from scipy.integrate import quad import matplotlib.pyplot as plt # 定义函数k乘x的n次方，自变量x def myfunction(x, n, k): return k * x**n # 存放积分结果 result = [] # 循环，查看0到14阶的x^(n)积分的曲线比较 for i in range(15): # 求每阶的积分结果, 下限0上限2, n=i, k=3 data = quad(func=myfunction, a=0, b=2, args=(i,3)) # 将每次循环的结果保存, 不要积分误差 result.append(data[0]) # 打印结果 print(result) # 绘制积分结果，查看阶数变化导致的积分结果的变化 plt.plot(result) plt.show()

案例四：实际应用

问题：

A公司的每个月销售额为100万元，平均利润为销售额的10%。根据公司以往的运营经验，广告宣传期间的月销售额的变化率近似服从增长曲线 $1000000*e^{0.02t}$ ，其中 t 为广告持续长度，以月份为单位。A公司正在决定是否要投入一次持续12个月，总成本为130000元的广告活动。依照公司规定，如果新增销售额的利润超过广告投资的10%，则可以投入该广告，请问该公司是否应当进行此次广告活动。

答：被积函数为 $\int_{0}^{12} 10000000\times e^{0.002t} dt$

已知广告宣传期间每个月的销售额变化率，可根据变化率求解广告期间的销售总量，平均利润占销售总量的10%，求得广告期间的总利润，减去不打广告期间的总利润，即可求得打广告期间的净收入利润。

from scipy.integrate import quad import math # 定义月销售额函数 def sales(t): # 自变量t代表做了几个月广告 # 定义销售额随广告时间的变化函数 data = 1000000 * (math.e**(0.02*t)) return data # 查看12个月的销售额增长曲线 result = [] for t in range(12): # 计算每个月的销售额 sell = 1000000 * (math.e**(0.02*t)) # 保存积分结果 result.append(sell) # 绘图查看增长率 plt.plot(result) plt.show() # 通过增长曲线计算总销售额 # 积分计算打广告期间从开始到结果的总销售额, 下限0上限12，不需要查看误差 revenue = quad(func=sales, a=0, b=12)[0] print('打广告带来的收益：', revenue) # 13562457.516070236 # 计算打广告和不打广告相比，新增的利润 revenue_add = revenue - 1000000 * 12 # 总销售额增加了多少 profit_add = revenue_add * 0.1 # 打广告带来的新增利润 print('打广告带来的新增的利润:', profit_add) # 156245.75160702356 # 计算打广告带来的新增利润是否超过广告投入成本 net_profit = profit_add - 130000 print('广告净收入:', net_profit) # 26245.751607023558 3. 二重积分

二重积分是二元函数在空间上的积分。本质是求曲顶柱体的体积，重积分可以用来计算曲面面积，平面薄片重心等。

定积分的计算方法使用 scipy.integrate 模块下的 dblquad 函数。函数的参数和返回值如下，gfun积分区间下限在图中表示 $\Psi 1(x)$ ， hfun积分区间下限在图中表示 $\Psi 2(x)$

''' 参数 # ---------------------------------------------------------- # func 被积函数 a 积分下限 b 积分上限 gfun 积分区间下限 hfun 积分区间上限 args 元组。func被积函数的未知变量xy以外，需要的其他变量的参数， # ---------------------------------------------------------- # 返回值：(积分结果, 误差值) 误差用来辅助判断这个积分结果准不准 ''' 案例一：计算空间斜平面的体积

现在有一个空间斜平面 $z=4-x-y$ ，计算该平面在 $x\in [0, 2], y\in [0,1]$ 范围内的体积，如下图。

计算公式为： $\int_{0}^{2}dx\int_{0}^{1}(4-x-y)dy$

先对 y 在 [gfun, hfun] 区间内求积分，再对 x 在 [a, b] 区间内求积分

from scipy.integrate import dblquad # 定义空间平面 def myfunction(x, y): # z = 4-x-y的平面 return 4-x-y # 计算二重积分, x的下限为0上限为2，y的下限为0上限为1 # 先对y积分，再对x积分 value, error = dblquad(func=myfunction, a=0, b=2, gfun=0, hfun=1) print('value:', value, 'error:', error) # value: 5.0 error: 5.551115123125783e-14 案例二：计算扇形区域面积

计算 $\int \int 3x^{2}y^{2}d\sigma$ ，其中D是由 x 轴，y 轴和抛物线 $y=1-x^{2}$ 在第一象限内所围成的区域。如下图所示。

由于此时的积分上下限 [a, b] 和积分区间上下限 [gfun, hfun] 相互依赖。因此，先对 y 在 [x,1-x2] 区间内求积分，后对 x 在 [0, 1] 区间内求积分，计算公式如下：

$\int_{0}^{1}3x^{2}dx\int_{x}^{1-x^{2}}y^{2}dy$

对y积分的积分区域上限，与对x积分的积分区间有关，因此y的积分区域 hfun 应当用一个函数来表示

# 案例二：计算围成扇形面积的积分 from scipy.integrate import dblquad # 定义被积函数 def myfunction(x, y): data = 3 * (x**2) * (y**2) return data # 如果外围积分区间和内部积分区间是相互依赖的就要写成函数形式 # y的定义域函数 def y_area(x): return 1-x**2 # 计算二重积分, 先对y积分，后对x积分 # # x定义域[0, 1], y定义域[0, 1-x^2] value, error = dblquad(func=myfunction, a=0, b=1, gfun=0, hfun=y_area) print(value) # 0.050793650793650794

【本文地址】

公司简介

联系我们