浅谈sympy库的实际运用
本博客参考:《python数学实验与建模 科学出版社》
sympy库常见模块简介
sympy包含诸多功能,包括多项式、微积分、求解方程、离散数学等
模块名 | 描述 |
---|---|
abc | 符号变量模块 |
calculus | 积分运算 |
core | 基本加乘指数运算 |
discrete | 离散数学 |
functions | 基本函数及特殊函数 |
galgebra | 几何代数 |
geometry | 几何实体 |
integrals | 符号积分 |
matrices | 线性代数和矩阵 |
ntheory | 数论函数 |
physics | 物理学 |
plotting | 二维三维绘图 |
stats | 统计学 |
高等数学
求极限
(1)
lim
x
→
0
t
a
n
3
x
2
x
\lim \limits_{x \to 0}\frac{tan 3x}{2x}
x→0lim2xtan3x
(2)
lim
x
→
∞
(
1
−
1
x
)
x
\lim \limits_{x \to \infty}(1-\frac{1}{x})^x
x→∞lim(1−x1)x
from sympy import symbols,limit,oo
import sympy as sp
x=symbols('x')
limit(sp.tan(3*x)/(2*x),x,0)
>>> 3/2
limit(pow(1+1/x,x),x,oo)
>>> e
注意:在使用sympy库,一般情况下不需要引入numpy库,因为sympy主要是符号运算,numpy主要是数据运算,就对于上式中的自变量x来说,它是一个类似于字符串的数据类型,如果上式变成:
limit(np.tan(3*x)/(2*x),x,0)
,就会引发下面错误:AttributeError: ‘Mul’ object has no attribute ‘tan’,TypeError: loop of ufunc does not support argument 0 of type Mul which has no callable tan method
求偏导
已知函数 z = x 3 y − y 3 x z=x^3y-y^3x z=x3y−y3x,求 ∂ z ∂ x 、 ∂ z ∂ y 、 ∂ 2 z ∂ x 2 \frac{\partial z}{\partial x}、\frac{\partial z}{\partial y}、\frac{\partial^2 z}{\partial x^2} ∂x∂z、∂y∂z、∂x2∂2z
from sympy import diff,symbols
x,y=symbols('x y')
z=x**3*y-y**3*x
diff(z,x),diff(z,y),diff(z,x,2)
>>> (3*x**2*y - y**3, x**3 - 3*x*y**2, 6*x*y)
级数求和
求下面式子的累加和: ∑ k = 1 n ( k − 1 ) 2 = n ( n − 1 ) ( 2 n − 1 ) 6 \sum \limits_{k=1}^n(k-1)^2=\frac{n(n-1)(2n-1)}{6} k=1∑n(k−1)2=6n(n−1)(2n−1)
from sympy import symbols,summation,factor
k,n=symbols('k n')
summation((k-1)**2,(k,1,n))
>>> 𝑛^3/3−𝑛^2/2+𝑛/6
factor(summation((k-1)**2,(k,1,n)))# 因式分解
>>> 𝑛(𝑛−1)(2𝑛−1)/6
泰勒展开式
泰勒公式的意义在于利用多项式函数逼近比较复杂的原函数,利用多项式函数可以任意次求导,易于计算,且便于求解极值或者判断函数的性质,更容易地获取原函数的信息
求xsinx在0点处的3,5,7阶泰勒展开式
解: 三阶: x 2 + O ( x 3 ) 五阶: x − x 3 6 + O ( x 5 ) 七阶: x − x 3 6 + x 5 120 + O ( x 7 ) 三阶:x^2+O(x^3)\\五阶:x-\frac{x^3}{6}+O(x^5)\\七阶:x-\frac{x^3}{6}+\frac{x^5}{120}+O(x^7) 三阶:x2+O(x3)五阶:x−6x3+O(x5)七阶:x−6x3+120x5+O(x7)
程序实现
from sympy import series,sin
x=symbols('x')
y=x*sin(x)
series(y,x,0,3)
>>> 𝑥^2+𝑂(𝑥^3)
series(y,x,0,5)
>>>𝑥−𝑥^3/6+𝑂(𝑥^5)
series(y,x,0,7)
>>>𝑥−𝑥^3/6+𝑥^5/120+𝑂(𝑥^7)
不定积分、定积分
(1)
∫
0
1
e
x
d
x
=
2
\int_0^1e^{\sqrt x}dx=2
∫01exdx=2
(2)
∫
0
+
∞
s
i
n
x
x
d
x
=
π
2
\int^{+\infty}_0\frac{sin x}{x}dx=\frac{\pi}{2}
∫0+∞xsinxdx=2π
from sympy import integrate,symbols,sin,oo,E,sqrt
x=symbols('x')
integrate(E**sqrt(x),(x,0,1))
>>> 2
integrate(sin(x)/x,(x,0,oo))
>>> 𝜋/2
求解代数方程
(1)求解方程组
{
x
2
+
y
2
=
8
x
+
y
=
0
\begin{cases} x^2+y^2=8\\ x+y=0 \end{cases}
{x2+y2=8x+y=0
from sympy.abc import x,y
from sympy import solve
s=solve([x**2+y**2-8,x+y],[x,y])
s
>>> [(-2, 2), (2, -2)]
(2)求函数 f ( x ) = 2 x 3 − 5 x 2 + x f(x)=2x^3-5x^2+x f(x)=2x3−5x2+x的驻点
函数的一阶导数为0的点(驻点也称为稳定点,临界点)。对于多元函数,驻点是所有一阶偏导数都为零的
求解微分方程
求下列微分方程的通解
(1)
3
x
2
+
5
x
−
5
y
′
=
0
3x^2+5x-5y'=0
3x2+5x−5y′=0
(2)
y
′
′
=
y
′
+
x
y''=y'+x
y′′=y′+x
(3)
y
′
′
−
5
y
′
+
6
y
=
0
,
其中
y
(
0
)
=
1
,
y
′
(
0
)
=
0
y''-5y'+6y=0,其中y(0)=1,y'(0)=0
y′′−5y′+6y=0,其中y(0)=1,y′(0)=0
from sympy import symbols,diff,Function,dsolve
x=symbols('x')
y=symbols('y',cls=Function)
ep1=3*x**2+5*x-5*diff(y(x),x)
dsolve(ep1,y(x))
>>> 𝑦(𝑥)=𝐶1+𝑥^3/5+𝑥^2/2
ep2=diff(y(x),x,2)-diff(y(x),x)-x
dsolve(ep2,y(x))
>>> 𝑦(𝑥)=𝐶1+𝐶2e^x−𝑥^2/2−𝑥
ep4=diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)
dsolve(ep4,y(x),cls={y(0):1,diff(y(x),x).subs(x,0):0})
>>> y(x)=(C1+C2e^x)e^(2x)
绘图(plotting子模块)
二维曲线
在同一坐标轴上画出 y 1 = 2 s i n x , x ∈ [ − 6 , 6 ] , y 2 = c o s ( x + π 2 ) , x ∈ [ − 5 , 5 ] y_1=2sinx,x\in [-6,6],y_2=cos(x+\frac{\pi}{2}),x\in[-5,5] y1=2sinx,x∈[−6,6],y2=cos(x+2π),x∈[−5,5]的图像
from sympy.plotting import plot
from sympy import sin,cos
plot((2*sin(x),(x,-6,6)),(cos(x+pi/2),(x,-5,5)))
运行结果
三维图像
绘制出三维曲面 f ( z ) = s i n ( x 2 + z 2 ) f(z)=sin(\sqrt {x^2+z^2}) f(z)=sin(x2+z2)
from sympy.plotting import plot3d
z=symbols('z')
from sympy import sqrt
plot3d(sin(sqrt(x**2+z**2)),(x,-10,10),(z,-10,10))
运行结果
隐函数
绘制出隐函数 ( x − 1 ) 2 + ( y − 2 ) 3 − 4 = 0 (x-1)^2+(y-2)^3-4=0 (x−1)2+(y−2)3−4=0的图像
from sympy import plot_implicit as pt,Eq
pt(Eq((x-1)**2+(z-2)**3,4),(x,-6,6),(z,-2,4))
运行结果
线性代数
矩阵运算
sympy库中矩阵由类sp.Matrix定义
矩阵类方法或属性 | 描述 | 数学符号 |
---|---|---|
A.rank() | 矩阵的秩 | R ( A ) R(A) R(A) |
A.T/A.transpose() | 矩阵A的转置矩阵 | A T A^T AT |
A.inv() | A的逆矩阵 | A − 1 A^{-1} A−1 |
sp.det(A) | 矩阵A的行列式 | ∣ A ∣ |A| ∣A∣ |
A.norm() | 矩阵A的模 | ∣ ∣ A ∣ ∣ 2 ||A||^2 ∣∣A∣∣2 |
解线性方程组
求下面齐次线性方程组的基础解系
{
2
x
1
+
x
2
−
2
x
3
+
3
x
4
=
0
3
x
1
+
2
x
2
−
x
3
+
2
x
4
=
0
x
1
+
x
2
+
x
3
−
x
4
=
0
\begin{cases} 2x_1+x_2-2x_3+3x_4=0\\ 3x_1+2x_2-x_3+2x_4=0\\ x_1+x_2+x_3-x_4=0 \end{cases}
⎩
⎨
⎧2x1+x2−2x3+3x4=03x1+2x2−x3+2x4=0x1+x2+x3−x4=0
C=sp.Matrix([[2,1,-2,3],[3,2,-1,2],[1,1,1,-1]])
C.nullspace()
>>>[Matrix([
[ 3],
[-4],
[ 1],
[ 0]]),
Matrix([
[-4],
[ 5],
[ 0],
[ 1]])]
求下面方程组的通解
{
x
1
+
2
x
2
−
x
3
+
3
x
4
=
2
2
x
1
+
4
x
2
−
2
x
3
+
5
x
4
=
1
−
x
1
−
2
x
2
+
x
3
−
x
4
=
4
\begin{cases} x_1+2x_2-x_3+3x_4=2\\ 2x_1+4x_2-2x_3+5x_4=1\\ -x_1-2x_2+x_3-x_4=4 \end{cases}
⎩
⎨
⎧x1+2x2−x3+3x4=22x1+4x2−2x3+5x4=1−x1−2x2+x3−x4=4
import sympy as sp
D=sp.Matrix([[1,2,-1,3],[2,4,-2,5],[-1,-2,1,-1]])
E=sp.Matrix([2,1,4])
F=D.row_join(E)
F.rref()# 此处笔者只将增广矩阵转化为行最简形,不对下面步骤进行赘述
>>> (Matrix([
[1, 2, -1, 0, -7],
[0, 0, 0, 1, 3],
[0, 0, 0, 0, 0]]),
(0, 3))
特征值及特征向量
求下列矩阵的特征值及特征向量
M
=
[
−
2
1
1
0
2
0
−
4
1
3
]
M=\begin{bmatrix} -2 & 1 & 1\\\\ 0 & 2 & 0\\\\ -4 & 1 & 3\\\\ \end{bmatrix}
M=⎣
⎡−20−4121103⎦
⎤
import sympy as sp
M=sp.Matrix([[-2,1,1],[0,2,0],[-4,1,3]])
M.eigenvals()# 特征值
>>> {2: 2, -1: 1}
M.eigenvects()# 特征向量
>>> [(-1, 1, [Matrix([
[1],
[0],
[1]])]),
(2,
2,
[Matrix([
[1/4],
[ 1],
[ 0]]),
Matrix([
[1/4],
[ 0],
[ 1]])])]
注意,eigenvals()函数返回值:{2:2,-1:1},它表示{特征值:特征值对应的特征向量数量},即 λ 1 = λ 2 = 2 \lambda_1=\lambda_2=2 λ1=λ2=2对应两个特征向量 [ 1 / 4 , 1 , 0 ] 、 [ 1 / 4 , 0 , 1 ] [1/4,1,0]、[1/4,0,1] [1/4,1,0]、[1/4,0,1],而 λ 3 = − 1 \lambda_3=-1 λ3=−1对应一个特征向量 [ 1 , 0 , 1 ] [1,0,1] [1,0,1]