本文主要是介绍Python数值微积分,摆脱被高数支配的恐惧,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
文章目录
- 差分和累加
- 积分
- 多重积分
Python科学计算:数组💯数据生成
差分和累加
微积分是现代科学最基础的数学工具,但其应用对象往往是连续函数,而其在非连续函数的类比,便是差分与累加。在【numpy】中,可通过【diff】和【cumsum】来完成这两项任务。
以 y = sin 2 x y=\sin 2x y=sin2x为例,其导数为 d y d x = 2 cos x \frac{\text dy}{\text dx}=2\cos x dxdy=2cosx,积分则为 ∫ y d x = − 1 2 cos 2 x + C \int y\text dx=-\frac{1}{2}\cos 2x+C ∫ydx=−21cos2x+C, C C C是某个常数。这三个函数的曲线分别为
绘图函数如下
import matplotlib.pyplot as plt
import numpy as np
dx = 0.1
x = np.arange(100)*dx
y = np.sin(2*x)
plt.plot(x, y, label="y=sin(2x)")
plt.plot(x[1:], np.diff(y)/dx, label="diff(y)/dx")
plt.plot(x, np.cumsum(y)*dx, label="cumsum(y)*dx")plt.grid()
plt.legend()
plt.show()
其中,diff用于求差分,其输入参数除了待差分数组之外,还有n和axis,比较常用,n为差分的阶数,默认为1;axis用于高维数组中,表示计算的方向,默认-1表示最后一个轴。
cumsum用于累加,对于输入数组 y y y,其返回数组为 S S S,则 S n = ∑ i = 0 n y i S_n=\sum_{i=0}^ny_i Sn=∑i=0nyi。
无论diff还是cumsum,均只针对输入数组进行操作,而不会考虑微积分中至关重要的 d x \text dx dx,所以绘图时对这一部分进行了补全。
此外,由于差分的实质是后一个减去前一个,所以元素个数必然会减少,所以在绘图时,令 x x x从1开始。这是一个在编程时很容易出错的地方,故而numpy还提供了另一个函数【ediff1d】,这是一个只做一阶差分计算的函数,但提供了to_end
和to_begin
参数,分别用于在diff计算结果的后面或前面补充数值。
积分
积分一开始被引入教材,是以梯形求和为示例的:将函数 y = f ( x ) y=f(x) y=f(x)无限分割,然后对相邻两点取平均,再乘以 d x \text dx dx之后进行求和,即 lim δ x → 0 ∑ y i + y i + 1 2 δ x \lim_{\delta_x\to0}\sum \frac{y_{i}+y_{i+1}}{2}\delta_x limδx→0∑2yi+yi+1δx。
【trapz】可实现上述过程,但要求 y y y是一个给定的数组,且 δ x \delta_x δx为1。很显然,这个过程只能称之为梯形求和,毕竟积分的要求是 δ x → 0 \delta_x\to0 δx→0, 1 1 1和 0 0 0有着本质的区别。
为此,【scipy.intergrate】作为顾名思义的积分模块,提供了真真正正的积分。为了行文简洁,后文将此模块简称为【si】模块。
【quad】是【si】中最常用的积分函数,以函数 x 2 x^2 x2和 sin x \sin x sinx为例,其使用流程如下
import numpy as np
from scipy.integrate import quadfunc = lambda x: x**2
quad(func, 0, 4) # (21.33, 2.37-13)
quad(np.sin, 0, np.pi) # (2.0, 2.22e-14)
其中,quad共输入了三个参数,分别是待积分函数、积分下界与积分上界,其返回值有二,分别为积分结果和计算误差。
这两个测试函数的解析形式如下,可见计算结果基本温和。
∫ 0 4 x 2 d x = 1 3 x 3 ∣ 0 4 = 64 3 ≈ 21.3 ∫ 0 π sin x d x = − cos x ∣ 0 π = 2 \int_0^4 x^2\text dx=\frac{1}{3}x^3\big|^4_0=\frac{64}{3}\approx 21.3\\ \int^\pi_0\sin x\text dx=-\cos x\big|^\pi_0=2 ∫04x2dx=31x3 04=364≈21.3∫0πsinxdx=−cosx 0π=2
除了三个必须输入的参数之外,下列参数也较为常用
args
为func
函数中,除待求积分参数之外的其他参数,默认为空epsabs, epsrel
分别为绝对和相对误差,默认为 1.49 × 1 0 − 8 1.49\times10^{-8} 1.49×10−8limit
自适应算法中子区间的个数,默认50points
断点位置,默认为Noneweight, wvar
定义域区间内的权重类型和权重,默认为Nonewopts, maxp1
切比雪夫矩及其上限,默认为None和50- full_output=0, limlst=50, complex_func=False
其中,weight
和wvar
参数的具体取值如下。
weight | wvar | 函数 |
---|---|---|
“cos” | w w w | cos w x \cos wx coswx |
“sin” | w w w | sin w x \sin wx sinwx |
“alg” | α , β \alpha, \beta α,β | g ( x ) g(x) g(x) |
“alg-loga” | α , β \alpha, \beta α,β | g ( x ) log ( x − a ) g(x)\log(x-a) g(x)log(x−a) |
“alg-logb” | α , β \alpha, \beta α,β | g ( x ) log ( b − x ) g(x)\log(b-x) g(x)log(b−x) |
“alg-log” | α , β \alpha, \beta α,β | g ( x ) log ( x − a ) log ( b − x ) g(x)\log(x-a)\log(b-x) g(x)log(x−a)log(b−x) |
“cauchy” | c c c | 1 x − c \frac{1}{x-c} x−c1 |
其中, g ( x ) = ( x − a ) α ∗ ( b − x ) β g(x)=(x-a)^\alpha*(b-x)^\beta g(x)=(x−a)α∗(b−x)β
设func
为 f ( x ) = x f(x)=x f(x)=x,若weight
参数为cos
,而wvar
取值为 w w w,则实际计算的积分表达式为
∫ a b cos w f ( x ) d x \int_a^b\cos wf(x)\text dx ∫abcoswf(x)dx
示例如下
func = lambda x : x
quad(func, 0, np.pi) # (4.935, 5.478e-14)
quad(func, 0, np.pi, weight='cos', wvar=1) # (-2.00, 1.926e-13)
多重积分
在【si】中,除了quad之外,还提供了二重、三重以及N重积分的API,分别是【dblquad, tplquad, nquad】,三者所需参数如下
MIN = 1.49e-08
dblquad(func, a, b, gfun, hfun, args=(), epsabs=MIN, epsrel=MIN)
tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=MIN, epsrel=MIN)
nquad(func, ranges, args=None, opts=None, full_output=False)
dblquad
以二重积分为例,其对应的问题可表述为下式
∫ a b ∫ y g ( x ) y h ( x ) f ( y , x ) d x d y \int^b_a\int^{y_h(x)}_{y_g(x)} f(y,x)\text dx\text dy ∫ab∫yg(x)yh(x)f(y,x)dxdy
在函数dblquad
中,func对应 f ( y , x ) f(y,x) f(y,x),a,b对那个上式的 a , b a,b a,b,gfun, hfun对应上式的 y g ( x ) , y h ( x ) y_g(x), y_h(x) yg(x),yh(x)。
接下来求解下面的积分
∫ 1 2 ∫ x 2 x 3 x y d y d x = ∫ 1 2 1 2 ( x y 2 ) ∣ x 2 x 3 d x = ∫ 1 2 1 2 ( x 7 − x 5 ) d x = 1 2 ( 1 8 x 8 − 1 6 x 6 ) ∣ 1 2 = 1 2 ( 2 8 8 − 2 6 6 ) + 1 48 = 513 48 \begin{aligned} &\int^2_1\int^{x^3}_{x^2} xy\text dy\text dx\\ =&\int^2_1 \frac{1}{2}(xy^2)\vert^{x^3}_{x^2}\text dx=&\int^2_1 \frac{1}{2}(x^7-x^5)\text dx\\ =&\frac1 2(\frac1 8x^8-\frac1 6x^6)\vert^2_1=&\frac1 2(\frac{2^8}{8}-\frac{2^6}{6})+\frac{1}{48}\\ =&\frac{513}{48} \end{aligned} ===∫12∫x2x3xydydx∫1221(xy2)∣x2x3dx=21(81x8−61x6)∣12=48513∫1221(x7−x5)dx21(828−626)+481
Python代码如下
from scipy.integrate import dblquad
func = lambda x,y : x*y
gf = lambda x: x**2
hf = lambda x: x**3
dblquad(func, 1, 2, gf, hf)
# (10.6875, 5.284867210146833e-13)
计算结果与 513 48 \frac{513}{48} 48513一致。
与二重积分相比,三重积分tplquad只是多了一组qfun和rfun,相当于z处于qfun(x,y)和rfun(x,y)之间。
【nquad】貌似不支持回调函数,其参数ranges是元组的列表,每个元组代表对应未知量的取值范围。若将其映射为三重积分函数,则ranges可表示为 ( ( a 1 , b 1 ) , ( a 2 , b 2 ) , ⋯ , ( a n , b n ) ) ((a_1,b_1), (a_2, b_2),\cdots,(a_n, b_n)) ((a1,b1),(a2,b2),⋯,(an,bn))
下面仍以函数func为例,用nquad得出结果
from scipy.integrate import nquad
nquad(func, [[1,2], [3, 4]])
#(0.39276170758930756, 4.91851540406507e-15)
这篇关于Python数值微积分,摆脱被高数支配的恐惧的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!