Python学习-Scipy库求积分、复合积分、微分（半圆面积、球体体积、复合梯形积分、常微分方程求解）

Python学习-Scipy库求积分函数、微分函数

积分integrate 面积：通用积分函数quad()；体积：通用双积分函数dblquad()

目录

1、通用积分函数quad()，求半圆面积 x2+y2=r^2

2、通用双积分函数dblquad()，求球体体积

3、复合梯形积分 integrate.trapz()

4、求解常微分方程 solve_ivp()

导入库

from scipy import integrate
import matplotlib.pyplot as plt
import numpy as np

1、通用积分函数quad()，求半圆面积 x2+y2=r**2

参数介绍：
func: 指定被积函数
a,b: 积分上下限
args: 指定传递给func函数的参数，元组
full_output: 设置非0，返回积分操作相关信息，默认为0
…
返回值：
y: 积分的结果值
abserr: 估计结果中的绝对误差
message: 收敛信息
explain: 仅附加’cos’和’sin’加权和无限积分限制的相关解释信息

# 1、通用积分函数quad()
"""
参数介绍：
func: 指定被积函数
a,b: 积分上下限
args: 指定传递给func函数的参数，元组
full_output: 设置非0，返回积分操作相关信息，默认为0
...
返回值：
y: 积分的结果值
abserr: 估计结果中的绝对误差
message: 收敛信息
explain: 仅附加'cos'和'sin'加权和无限积分限制的相关解释信息
"""
# 求半圆面积 x^2+y^2=r^2
a = -2  # 积分下限
b = 2  # 积分上限


def m_circle(x, r):
    return np.sqrt(r**2 - x**2)


area1, _ = integrate.quad(m_circle, a, b, 2)
print(area1)
plt.figure(num=1, figsize=(5, 5))
plt.rc('font', family='simhei', size=15)  # 设置中文显示，字体大小
plt.rc('axes', unicode_minus=False)  # 该参数解决负号显示的问题
x = np.linspace(-2, 2, 180)
y = np.sqrt(4-x**2)
plt.plot(x, y, 'm', lw=2)
plt.fill(x, y, color='m', alpha=0.7)
plt.title('单位半圆面积为'+'{:.2f}'.format(area1))
plt.show()

输出

6.283185307179593

2、通用双积分函数dblquad()，求球体体积

参数介绍：
func: 含两个变量的被两重积分函数，第一参数y，第二参数x
a,b: 指定x中的积分限制a<b区域
gfun: 指定y的下边界曲线的函数，必须包含一个浮点型的x参数
hfun: 指定y的上边界曲线的函数，必须包含一个浮点型的x参数
args: 指定传递给func函数的额外参数（指第三个开始的参数）
epsabs: 直接传递到内部一维正交积分的绝对容差，默认1.49e-8
epsrel: 内部一维积分的相对容差，默认1.49e-8
…
返回值：
y: 积分的结果值
abserr: 估计结果中的绝对误差

# 2、通用双积分函数dblquad()
"""
参数介绍：
func: 含两个变量的被两重积分函数，第一参数y，第二参数x
a,b: 指定x中的积分限制a<b区域
gfun: 指定y的下边界曲线的函数，必须包含一个浮点型的x参数
hfun: 指定y的上边界曲线的函数，必须包含一个浮点型的x参数
args: 指定传递给func函数的额外参数（指第三个开始的参数）
epsabs: 直接传递到内部一维正交积分的绝对容差，默认1.49e-8
epsrel: 内部一维积分的相对容差，默认1.49e-8
...
返回值：
y: 积分的结果值
abserr: 估计结果中的绝对误差
"""
# 求球体体积


def g_func_c(x):
    return -np.sqrt(1-x**2)


def h_func_c(x):
    return np.sqrt(1-x**2)


def func_s(x,y):
    return np.sqrt(1-x**2-y**2)


volume, _ = integrate.dblquad(func_s, -1, 1, g_func_c, h_func_c)
volume, _ = integrate.dblquad(func_s, -1, 1, lambda x: -h_func_c(x), lambda x: h_func_c(x))

print(volume*2)
print(np.pi*4/3*1**3)

输出

4.188790204786397
4.1887902047863905

3、复合梯形积分 integrate.trapz()

参数说明：
y: 样本集合对象
x: 对应样本点的x值
dx: 指定采样点之间的间距
axis: 指定积分的轴

# 复合梯形积分  integrate.trapz(y,x=None,dx=1.0,axis=-1)
"""
参数说明：
y: 样本集合对象
x: 对应样本点的x值
dx: 指定采样点之间的间距
axis: 指定积分的轴
"""
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)
area = integrate.trapz(y=np.abs(y), x=x)
print(area)
plt.plot(x, y, 'g--', linewidth=2)
plt.fill(x, y, color='b', alpha=0.7)
plt.show()

输出

3.9996643277879853

4、求解常微分方程 solve_ivp()

参数说明：
fun: 调用fun(t,y)函数，t是一个数值常量；y是一个数组
t_span: 2值元组，提供积分区间范围（t0,tf）
y0: 初始状态
method: 求解方法：‘RK45’ ‘RK23’ ‘Radau’ ‘BDF’ ‘LSODA’
t_eval: 计算连续解决方案存储需要的时间，受t_span影响
dense_output: 是否计算连续解决方案
…
返回值：
t: 时间点
y: t处的解
…

# 求解常微分方程  solve_ivp()
"""
参数说明：
fun: 调用fun(t,y)函数，t是一个数值常量；y是一个数组
t_span: 2值元组，提供积分区间范围（t0,tf）
y0: 初始状态
method: 求解方法：'RK45' 'RK23' 'Radau' 'BDF' 'LSODA'
t_eval: 计算连续解决方案存储需要的时间，受t_span影响
dense_output: 是否计算连续解决方案
...
返回值：
t: 时间点
y: t处的解
...
"""
plt.rc('font', family='simhei', size=15)  # 设置中文显示，字体大小
plt.rc('axes', unicode_minus=False)  # 该参数解决负号显示的问题


def f1(t, y):
    return np.sin(t**2+2)


t1 = np.linspace(-10, 10, 1000)
y0 = np.array([9])  # 初始值
y1 = integrate.solve_ivp(f1, (-10, 10), y0, method='LSODA', t_eval=t1)  # 求解
plt.plot(y1.t, y1.y[0], 'b')
plt.title('求一阶常微分方程')
plt.show()

输出