微积分相关
2024年11月22日大约 4 分钟
积分
scipy.integrate.quad
用于计算定积分。这个函数可以用来计算一个函数在给定区间上的定积分值。并返回积分的结果和误差估计。
from scipy.integrate import quad
# 定义被积函数
def f(x):
return 3*x**2
print(quad(f, 0, 1))(1.0, 1.1102230246251565e-14)梯度
np.gradient
用于计算数组的梯度。梯度是一个矢量，其分量是函数在各个方向上的偏导数。np.gradient 函数可以计算数组的梯度，并返回一个包含各个方向梯度的数组。
import numpy as np
x = np.array([1, 2, 4, 7, 11])
print(np.gradient(x))[1. 1.5 2.5 3.5 4. ]可以用这个函数来计算函数的导数。
x = np.arange(0,2,0.01)
y = x**2
dydx = np.gradient(y,x)
print(dydx)
print(f"在x=1处的导数值是：{round(dydx[100],4)}")[0.01 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26
0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54
...
3.92 3.94 3.96 3.97]
在x=1处的导数值是：2.0微分方程
scipy.integrate.odeint
用来求微分方程数值解
odeint()函数需要至少三个变量，第一个是微分方程函数，第二个是微分方程初值，第三个是微分的自变量
import matplotlib.pyplot as plt
import scipy.integrate as spi
import numpy as np
dy = lambda y,x: 1/(1+x**2)-2*y**2 # 定义y'
x = np.arange(0,10.5,0.1)
sol = spi.odeint(dy,0,x)
print(sol[:10])[[0. ]
[0.09900996]
[0.19230775]
[0.2752294 ]
[0.34482762]
[0.40000004]
[0.44117651]
[0.46979872]
[0.48780495]
[0.49723764]]对于高阶微分方程
import numpy as np
def fvdp(y, x):
'''
要把y看出一个向量，y = [dy0,dy1,dy2,...]分别表示y的n阶导，那么
y[0]就是需要求解的函数，y[1]表示一阶导，y[2]表示二阶导，以此类推
'''
dy1 = y[1] # y[1]=dy/dt，一阶导 y[0]表示原函数
dy2 = 20*(1-y[0]**2) * y[1] - y[0] # y[1]表示一阶微分
# y[0]是最初始，也就是需要求解的函数
# 注意返回的顺序是[一阶导， 二阶导]，这就形成了一阶微分方程组
return [dy1, dy2]
x = np.arange(0, 0.25, 0.01) # 给x规定范围
y0 = [0.0, 2.0] # 初值条件
# 初值[3.0, -5.0]表示y(0)=3,y'(0)=-5
# 返回y，其中y[:,0]是y[0]的值，就是最终解，y[:,1]是y'(x)的值
y = spi.odeint(fvdp, y0, x)
print("y:")
print(y[:, 0])
print("y':")
print(y[:, 1])y:
[0. 0.02213973 0.04917542 0.08217383 0.12241129 0.17139008
0.23083061 0.30261249 0.38862014 0.49043034 0.6087762 0.7427756
0.8890802 1.041387 1.19095297 1.32843735 1.44644047 1.54133536
1.61342302 1.6657038 1.70228635 1.72720565 1.74383888 1.7547616
1.76182797]
y':
[ 2. 2.44261521 2.98225626 3.63866651 4.43386984 5.39064944
6.52902659 7.85925963 9.36943546 11.00614061 12.64985337 14.09520735
15.05991697 15.25250139 14.50038395 12.8698378 10.67018007 8.3161543
6.15422899 4.37240537 3.01227868 2.02805197 1.34134997 0.87382356
0.56063686]scipy.integrate.solve_ivp
Solve_ivp函数的用法与odeint非常类似
不过f定义为(t,y)参数
同时比odeint多了两个参数。一个是t_span参数，表示自变量的取值范围；
另一个是method参数，可以选择多种不同的数值求解算法。
常见的内置方法包括RK45, RK23, DOP853, Radau, BDF等多种方法，通常使用RK45多一些。
返回的值上也有区别
odeint 返回一个数组，其中包含时间点和对应的解。 solve_ivp 返回一个 Solution 对象，提供了更多关于解的信息，例如时间点、解的值、求解过程中的事件等。
def f(t,y):
dy1 = y[1]
dy2 = y[2]
dy3 = -y[0]+dy1-dy2-np.cos(t)
return [dy1, dy2,dy3]
t = np.linspace(0, 6, 100)
tspan = (0.0, 6.0)
y0 = [0.0, np.pi, 0.0]
y_ = spi.solve_ivp(f, t_span=tspan, y0=y0, t_eval=t)
print("y:\n",y_.y[0,:])
print("y':\n", y_.y[1, :])
print("y'':\n", y_.y[2, :])
y:
[ 0. 0.19047725 0.3813934 0.57312021 0.76597764
0.96023627 1.15612795 1.35380768 1.5533875 1.75495329
...
-8.93373883 -10.14640863 -11.40213181 -12.70065003 -14.04160799
-15.42455339 -16.84893696 -18.31411245 -19.81933663 -21.3637693 ]
y':
[ 3.14159265 3.14533218 3.155821 3.17205015 3.19308723
3.21807175 3.24619388 3.27677009 3.30917626 3.34281454
3.37711336 3.41152743 3.44553775 3.47866472 3.51050446
...
-10.9388097 -11.098204 -11.2427817 -11.37172317 -11.48447992
-11.57987618 -11.65671337 -11.713771 -11.74980666 -11.7635561
-11.75373312 -11.71902967 -11.65811577 -11.56963957 -11.45222732
-11.30448336 -11.12499015 -10.91230826 -10.66497635 -10.3815112 ]相关信息
odeint内置的原理也是4-5阶 Runge-Kutta 法，版本比较早所以求解也相对较为稳定。但solve_ivp则是后来新增的方法，有可能出现不太稳定的现象，内置方法较多所以也更加灵活。
使用的区别
def dmove(Point, t, sets):
p, r, b = sets
x, y, z = Point
return np.array([p*(y-x), x*(r-z), x*y-b*z])
def dmove2(t,yy,sets):
p,r,b = sets
x,y,z = yy
return [p*(y-x),x*(r-z)-y,x*y-b*z]
t = np.arange(0, 30, 0.001)
P1 = spi.odeint(dmove, (0., 1., 0.), t, args=([10., 28., 3.],))
P2 = spi.solve_ivp(dmove2, (0,30) ,(0., 1.01, 0.),t_eval=t, args=([10., 28., 3.],))
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(P1[:, 0], P1[:, 1], P1[:, 2],color="pink")
ax.plot(P2.y[0], P2.y[1], P2.y[2],color="purple")
plt.show()