【Scipy高级计算】

大家好，偶然接触到了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以外，需要的其他变量的参数，
# ---------------------------------------------------------- #

返回值： (积分结果, 误差值)
误差用来辅助判断这个积分结果准不准
'''

---

案例一：计算  在 [0, 2] 积分区间的结果

计算公式为： 

被积函数中只有一个自变量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)

---

案例二：计算  在 [0, 2] 积分区间的结果

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

下面计算  在[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)

---

案例三：不同阶数对定积分结果的影响

现有被积函数  ，积分上下限[0, 1]，若k=3，研究不同的阶数n对定积分结果的影响，n=0,1,2，...。如下，我们使用一个循环，比较  到  的积分结果，将每次积分结果保存在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%。根据公司以往的运营经验，广告宣传期间的月销售额的变化率近似服从增长曲线  ，其中 t 为广告持续长度，以月份为单位。A公司正在决定是否要投入一次持续12个月，总成本为130000元的广告活动。依照公司规定，如果新增销售额的利润超过广告投资的10%，则可以投入该广告，请问该公司是否应当进行此次广告活动。

答：被积函数为  

已知广告宣传期间每个月的销售额变化率，可根据变化率求解广告期间的销售总量，平均利润占销售总量的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积分区间下限在图中表示， hfun积分区间下限在图中表示 

'''
参数
# ---------------------------------------------------------- #
func     被积函数
a        积分下限
b        积分上限
gfun     积分区间下限
hfun     积分区间上限
args     元组。func被积函数的未知变量xy以外，需要的其他变量的参数，
# ---------------------------------------------------------- #

返回值：(积分结果, 误差值)
误差用来辅助判断这个积分结果准不准
'''

---

案例一：计算空间斜平面的体积

现在有一个空间斜平面 ，计算该平面在  范围内的体积，如下图。

计算公式为：

先对 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

---

案例二： 计算扇形区域面积

计算   ，其中D是由 x 轴，y 轴和抛物线  在第一象限内所围成的区域。如下图所示。

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

对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