需要用到的工具有:sympy
,numpy
和scipy
库。
文章目录
一、用SymPy库求符号解
1. 求极限
求极限使用limit
函数,其函数调用格式为limit(e, z, z0, dir='+')
,用于计算函数e(z)
在z0
处的极限。
参数表:
e
:表达式,需要求极限的式子。z
:用于求极限的符号变量,表达式中其他符号变量视为常量。z0
:z
趋向的数值,无穷大表示为oo, -oo
。dir
:可选字符串类型参数,默认为+
,表示极限的方向。如果dir=+-
则计算双侧极限,dir=+
计算右极限,dir=-
计算左极限。
使用举例:
- 求
lim
x
→
0
s
i
n
x
x
\lim_{x\to 0}\frac{\mathrm{sin}x}{x}
limx→0xsinx的值,代码如下:
输出为1。from sympy import * from sympy.abc import x print(limit(sin(x)/x, x, 0))
- 求
lim
x
→
0
+
1
x
\lim_{x\to0^+}\frac{1}{x}
limx→0+x1与
lim
x
→
0
−
1
x
\lim_{x\to0^-}\frac{1}{x}
limx→0−x1,代码如下:
输出为from sympy import * from sympy.abc import x print(limit(1/x, x, 0, '+')) print(limit(1/x, x, 0, '-'))
oo
与-oo
。倘若将dir
参数改为+-
,则会弹出错误:ValueError: The limit does not exist since left hand limit = -oo and right hand limit = oo
。
2.求导数/偏导
求导数使用diff
函数,其调用格式为diff(f, *symbols, **kwargs)
,求f
关于symbols
的导数。也可以求高阶导数,一般的使用方式是diff(func, x, n)
,这里n
是可选的阶数,默认为一阶。
使用举例:
- 求
f
′
(
x
)
,
f
(
x
)
=
c
o
s
2
x
s
i
n
x
f'(x), f(x)=\mathrm{cos}^2x\mathrm{sin}x
f′(x),f(x)=cos2xsinx,代码如下:
输出为from sympy import * from sympy.abc import x f = cos(x)**2*sin(x) print(diff(f, x))
-2*sin(x)**2*cos(x) + cos(x)**3
。 - 已知
z
=
s
i
n
x
+
x
2
e
y
z = \mathrm{sin}x+x^2e^y
z=sinx+x2ey,求
∂
2
z
∂
x
∂
y
\frac{\partial^2{z}}{\partial{x}\partial{y}}
∂x∂y∂2z与
∂
2
z
∂
x
2
\frac{\partial^2{z}}{\partial{x}^2}
∂x2∂2z,代码如下:
输出为from sympy import * from sympy.abc import x, y z = sin(x) + x**2*exp(y) print(diff(z, x, y)) print(diff(z, x, 2))
2*x*exp(y)
与2*exp(y) - sin(x)
。
3.求积分
求积分包括定积分与不定积分,使用integrate
函数,其调用格式为integrate(f, var, ...)
,这里f
是被积函数,var
可以是单个函数变量(不定积分),也可以是一个给出上下限的元组(symbol, a, b)
(定积分)。
使用举例:
- 求
∫
x
y
d
x
\int xy \mathrm{d}x
∫xydx,代码如下:
输出为from sympy import * from sympy.abc import x, y f = x*y print(integrate(f, x))
x**2*y/2
。 - 求
∫
1
a
l
n
x
d
x
\int_1^a \mathrm{ln} x\mathrm{d}x
∫1alnxdx,代码如下:
输出为from sympy import * from sympy.abc import x, a f = log(x) print(integrate(f, (x, 1, a)))
a*log(a) - a + 1
。 - 求
∫
e
x
2
d
x
\int e^{x^2}\mathrm{d}x
∫ex2dx,代码如下:
输出为from sympy import * from sympy.abc import x f = exp (x**2) print(integrate(f, x))
sqrt(pi)*erfi(x)/2
。
4.求泰勒展开
求泰勒展开可以用series
函数,其调用格式为series(expr, x=None, x0=0, n=6, dir='+')
,参数释义如下:
expr
:被展开的表达式。x
:符号变量,表达式中对其展开。x0
:在何处展开,不能是oo
。n
:展开的阶数。dir
:可选字符串类型,默认为+
,表示泰勒展开的方向。
使用举例:求 s i n x \mathrm{sin}x sinx在 x = 0 x=0 x=0处的1,3,5,7阶展开。代码如下:
from sympy import *
from sympy.abc import x
y = sin(x)
for i in range(1, 8, 2):
print(series(y, x, 0, i))
输出结果为
O(x)
x + O(x**3)
x - x**3/6 + O(x**5)
x - x**3/6 + x**5/120 + O(x**7)
5.级数求和
级数求和使用summation
函数,调用格式为summation(f, (i, a, b))
,等价于计算
∑
i
=
a
b
f
(
i
)
。
\sum_{i=a}^bf(i)。
i=a∑bf(i)。
使用举例:
- 计算
∑
k
=
1
∞
1
k
2
\sum_{k=1}^{\infty}\frac{1}{k^2}
∑k=1∞k21,代码如下:
输出为from sympy import * from sympy.abc import k f = 1/k**2 print(summation(f, (k, 1, oo)))
pi**2/6
。 - 计算
∑
k
=
1
n
k
3
\sum_{k=1}^n k^3
∑k=1nk3,代码如下:
输出为from sympy import * from sympy.abc import k, n f = k**3 print(summation(f, (k, 1, n)))
n**4/4 + n**3/2 + n**2/4
。
6.多项式的处理
对于多项式,有一些处理的函数如下。
函数 | 作用 |
---|---|
expand | 展开多项式 |
factor | 因式分解 |
collect | 合并同类项 |
cancel | 有理分式化简 |
apart | 分式展开为真分式 |
使用举例:
from sympy import *
from sympy.abc import x, y
#将多项式展开
expr1 = (x**2+3*x)*(x**2+4)
print("expr1 = ", expand(expr1))
#因式分解
expr2 = x**4 + 2*x**3 + x**2
print("expr2 = ", factor(expr2))
#合并同类项
expr3 = x*y + 2*x**2 + x**2*y**2+3*x
print("expr3 = ", collect(expr3, x))
#有理分式化简
expr4 = (x**3 + 3*x**2 + 3*x + 1)/(x+1)
print("expr4 = ", cancel(expr4))
#化为真分式
expr5 = (x**4 + 4*x**3 + 4*x**2 + 2*x + 1)/(x**2 + x - 1)
print("expr5 = ", apart(expr5))
输出结果为:
expr1 = x**4 + 3*x**3 + 4*x**2 + 12*x
expr2 = x**2*(x + 1)**2
expr3 = x**2*(y**2 + 2) + x*(y + 3)
expr4 = x**2 + 2*x + 1
expr5 = x**2 + 3*x + 3*(x + 1)/(x**2 + x - 1) + 2
7.解方程(符号解)
解方程使用solve
函数,一般调用格式为solve(f, symbols)
,解出方程f(x)=0
的根。此函数也用于解多元函数组,这时候f
和symbols
均以列表的形式传入。
使用举例:
- 解方程
x
3
−
3
x
2
+
x
=
1
x^3-3x^2+x=1
x3−3x2+x=1,代码如下:
输出为from sympy import * from sympy.abc import x f = x**3 - 3**2 + x - 1 print(solve(f, x))
[2, -1 - 2*I, -1 + 2*I]
。 - 解方程组
{ x 2 + y 2 = 1 , x − y = 0. \left\{ \begin{array}{l} x^2+y^2=1,\\ x-y=0. \end{array} \right. {x2+y2=1,x−y=0.
代码如下:
输出为from sympy import * from sympy.abc import x, y f = [x**2+y**2-1, x-y] print(solve(f, [x, y]))
[(-sqrt(2)/2, -sqrt(2)/2), (sqrt(2)/2, sqrt(2)/2)]
。
8.求微分方程的符号解
求微分方程的符号解使用dsolve
函数,在声明符号变量时,要将函数指定为符号变量,这可以使用Function()
函数。如果是初值问题,还要用ics
参数指定初值。
使用举例:
- 求非齐次微分方程
y
′
′
−
5
y
′
+
6
y
=
x
e
2
x
y''-5y'+6y=xe^{2x}
y′′−5y′+6y=xe2x的解,代码如下:
输出为from sympy import * x = symbols('x') y = Function('y') #或y = symbols('y', cls=Function) eq = diff(y(x), x, 2) - 5*diff(y(x), x) + 6*y(x) - x*exp(2*x) print(dsolve(eq, y(x)))
Eq(y(x), (C1 + C2*exp(x) - x**2/2 - x)*exp(2*x))
,即 y ( x ) = ( C 1 + C 2 e x − x 2 2 − x ) e 2 x y(x)=(C_1+C_2e^x-\frac{x^2}2-x)e^{2x} y(x)=(C1+C2ex−2x2−x)e2x。 - 求初值问题
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 * x = symbols('x') y = Function('y') eq = diff(y(x), x, 2) - 5*diff(y(x), x) + 6*y(x) result = dsolve(eq, y(x), ics={y(0):1, diff(y(x),x).subs(x,0):0}) print(result)
Eq(y(x), (3 - 2*exp(x))*exp(2*x))
,即 y ( x ) = ( 3 − 2 e x ) e 2 x y(x)=(3-2e^x)e^{2x} y(x)=(3−2ex)e2x。注意,f.subs(x,0)
指将f
中的x
赋值为0。
二、用scipy库求数值解
使用scipy
库相比sympy
库的符号形式要显得麻烦一些,因为它依赖于numpy
的数组。
1.一重积分
求一重积分,用到scipy.integrate
库中的quad
函数,其调用方式为quad(f, a, b)
,这里f
是被积函数,a,b
是其积分下上限,返回一个元组,其中第一个数是积分值,第二个数是误差范围。注意f
要传入一个函数(lambda
函数或def
函数)而不是一个表达式。
使用示例:求
s
i
n
x
x
\frac{\mathrm{sin}x}{x}
xsinx在
[
0
,
1
]
[0,1]
[0,1]上的积分。代码如下:
import numpy as np
from scipy.integrate import quad
f = lambda x: np.sin(x)/x #这里使用的是numpy的sin函数
result = quad(f, 0, 1)
print(result)
输出结果为(0.9460830703671831, 1.0503632079297089e-14)
。
2.二重积分
求二重积分,使用的是scipy.integrate
中的dblquad()
函数,一般的调用格式为dblquad(func, a, b, gfun, hfun)
,这里被积函数func
的声明格式为func(y,x)
(内部参数的顺序按照累次积分的顺序),a,b,gfun,hfun
均为积分限,且按照积分的书写顺序来。
使用示例:求二重积分
I
=
∬
x
2
+
y
2
≤
1
e
−
x
2
2
s
i
n
(
x
2
+
y
)
d
x
d
y
I=\iint \limits_{x^2+y^2\leq1}e^{-\frac{x^2}{2}}\mathrm{sin}(x^2+y)dxdy
I=x2+y2≤1∬e−2x2sin(x2+y)dxdy
我们要先将其化成累次积分,为
I
=
∫
−
1
1
d
x
∫
−
1
−
x
2
1
−
x
2
f
(
x
)
d
y
I = \int_{-1}^{1}dx\int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}f(x)dy
I=∫−11dx∫−1−x21−x2f(x)dy
代码如下:
from scipy.integrate import dblquad
import numpy as np
f = lambda y,x: np.exp(-x**2/2)*np.sin(x**2+y) #注意参数顺序
bd = lambda x: np.sqrt(1-x**2) #积分限
result = dblquad(f, -1, 1, lambda x:-bd(x), bd) #这里不能直接写-bd
print(result)
运行结果为:(0.5368603826989582, 3.696155159715886e-09)
。
3.三重积分
三重积分使用scipy.integrate
中的tplquad
函数,其调用方式与dblquad
类似。
使用举例:求
I
=
∫
0
2
d
x
∫
−
2
x
−
x
2
2
x
−
x
2
d
y
∫
0
6
z
x
2
+
y
2
+
1
d
z
I = \int_0^2dx\int_{-\sqrt{2x-x^2}}^{\sqrt{2x-x^2}}dy\int_0^6z\sqrt{x^2+y^2+1}dz
I=∫02dx∫−2x−x22x−x2dy∫06zx2+y2+1dz
代码如下:
from scipy.integrate import tplquad
import numpy as np
f = lambda z,y,x: z*np.sqrt(x**2+y**2+1)
bd = lambda x: np.sqrt(2*x-x**2)
result = tplquad(f, 0, 2, lambda x:-bd(x), bd, 0, 6)
print(result)
输出为(87.45019779526699, 8.742462398458883e-08)
。
4.求非线性方程组的数值解
非线性方程组的求解,使用scipy.optimize
中的fsolve
函数,其调用格式为fsolve(func, x0, ...)
,求解的是func=0
。这里x0
是一个n维数组,代表方程的初始猜测值。
使用示例:求解非线性方程组
{
5
x
2
+
3
=
0
4
x
1
−
2
s
i
n
(
x
2
x
3
)
=
0
x
2
x
3
−
1.5
=
0
\left\{ \begin{array}{l} 5x_2+3=0\\ 4x_1-2\mathrm{sin}(x_2x_3)=0\\ x_2x_3-1.5=0 \end{array} \right.
⎩⎨⎧5x2+3=04x1−2sin(x2x3)=0x2x3−1.5=0
代码如下:
from scipy.optimize import fsolve
import numpy as np
f = lambda x: [5*x[1]+3, 4*x[0]-2*np.sin(x[1]*x[2]), x[1]*x[2]-1.5] #这里x是一个列表,f也是一个列表
result = fsolve(f, [1.0, 1.0, 1.0]) #猜测值列表是一个与x同维度的列表
print(result)
输出结果为[ 0.49874749 -0.6 -2.5 ]
。
5.求一元函数的极值点
求一元函数的极值点,用的是scipy.optimize
库中的fminbound
函数,其调用方式为fminbound(f, a, b)
,求函数f
在
[
a
,
b
]
[a,b]
[a,b]上的极小值点。
使用示例:求
x
2
−
2
x
+
1
x^2-2x+1
x2−2x+1在
[
−
5
,
5
]
[-5,5]
[−5,5]上的极小值点与极小值,代码如下:
from scipy.optimize import fminbound
import numpy as np
f = lambda x: x**2-2*x+1
result = fminbound(f, -5, 5)
print(result)
print(f(result))
输出结果为1.0, 0.0
。
6.求多元函数的极值点
求多元函数的极值点,用的是scipy.optimize
库中的minimize
函数,其调用方式为minimize(fun, x0, args=(), method=None)
。这里fun
是要求极小值的函数,x0
是初始猜测值,method
表示使用的算法。最后使用结果的x
和fun
属性得到结果。
使用示例:求
f
(
x
)
=
100
(
x
2
−
x
1
2
)
2
+
(
1
−
x
1
)
2
f(x)=100(x_2-x_1^2)^2+(1-x_1)^2
f(x)=100(x2−x12)2+(1−x1)2的极小值点与极小值,代码如下:
from scipy.optimize import minimize
import numpy as np
f = lambda x: 100*(x[1]-x[0]**2)**2+(1-x[0])**2
result = minimize(f, [2.0, 2.0])
print(result.x, result.fun)
输出结果为[0.9999956 0.9999912] 1.9394621611053384e-11
,已经很接近极小值点[1, 1]
。