SymPy包介绍#
SymPy 是一个用于符号数学计算的 Python 库。与数值计算库（如 NumPy 和 SciPy）不同，SymPy 的主要目标是提供一个完全符号化的系统，用于执行代数、微积分、方程求解等。
为什么使用SymPy
能够精确地处理代数表达式，不会出现数值误差。
可以与其他 Python 库（如 Matplotlib、NumPy）无缝集成。
提供了一个直观的 API，用于创建和操作数学表达式。
符号，表达式以及方程的定义#
from sympy import *
#符号的定义
x, y = symbols('x y')
x
\[\displaystyle x\]
from sympy import *
#符号的定义
beta = symbols('beta_1:4')
beta[0]
\[\displaystyle \beta_{1}\]
from sympy import *
#符号的定义
x, y = symbols('x y')
beta = symbols('beta_1:4')
#表达式的定义
expr = x**2 + 2*y + beta[0] + beta[1]*x + beta[2]*y
expr
\[\displaystyle \beta_{1} + \beta_{2} x + \beta_{3} y + x^{2} + 2 y\]
from sympy import *
#符号的定义
x, y = symbols('x y')
beta = symbols('beta_1:4')
#表达式的定义
expr = x**2 + 2*y + beta[0] + beta[1]*x + beta[2]*y
#等式的定义，注意，是用Eq函数定义的，并不是=或者==
eq = Eq(expr, 0)
eq
\[\displaystyle \beta_{1} + \beta_{2} x + \beta_{3} y + x^{2} + 2 y = 0\]
替换与化简#
sympy中可以使用subs函数用于替换
from sympy import *
x, y, z = symbols('x y z')
expr = cos(x**2) + 1
expr.subs(x, y)
\[\displaystyle \cos{\left(y^{2} \right)} + 1\]
from sympy import *
x, y, z = symbols('x y z')
expr = cos(x**2) + 1
expr.subs(x, z**2)
\[\displaystyle \cos{\left(z^{4} \right)} + 1\]
from sympy import *
x, y, z = symbols('x y z')
expr = cos(x**2) + 1
expr.subs(x**2, y)
\[\displaystyle \cos{\left(y \right)} + 1\]
from sympy import *
x, y, z = symbols('x y z')
expr = cos(x**2) + sin(x**2) + tan(x**3) + x*x
expr.subs(x**2, y)
\[\displaystyle y + \sin{\left(y \right)} + \cos{\left(y \right)} + \tan{\left(x^{3} \right)}\]
from sympy import *
x, y, z = symbols('x y z')
expr = cos(x**2) + sin(x**2) + tan(x**3) + x*x
expr.subs(x, sqrt(y))
\[\displaystyle y + \sin{\left(y \right)} + \cos{\left(y \right)} + \tan{\left(y^{\frac{3}{2}} \right)}\]
SymPy中提供了多种化简以及变换的函数，以及一个智能的化简函数simplify
from sympy import *
x = symbols('x')
simplify(sin(x)**2 + cos(x)**2)
\[\displaystyle 1\]
from sympy import *
x = symbols('x')
simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
\[\displaystyle x - 1\]
simplify有一个缺陷，由于表达式“最简化”并没有一个良好的定义，SymPy只能使用库中已有的化简操作，使用启发式方法来决定其认为的“最简化”结果。各种化简函数繁多，这里不再单独介绍，请有需要的同学移步官方文档Simplification
求值#
如果我们想要将一个数值表达式转化为一个浮点数，可以使用evalf函数
from sympy import *
expr = sqrt(8)
expr
\[\displaystyle 2 \sqrt{2}\]
from sympy import *
expr = sqrt(8)
expr.evalf()
\[\displaystyle 2.82842712474619\]
from sympy import *
expr = sqrt(8)
expr.evalf(3)
\[\displaystyle 2.83\]
from sympy import *
x = symbols('x')
expr = cos(2*x)
expr.evalf(subs={x: 2.4})
\[\displaystyle 0.0874989834394464\]
from sympy import *
one = cos(1)**2 + sin(1)**2
print((one - 1).evalf())
print((one - 1).evalf(chop=True))
-0.e-124
0
如果我们想要再多个点对表达式进行求值，我们可以使用lambdify对其进行求值
from sympy import *
import numpy as np
from matplotlib import pyplot as plt
x = symbols('x')
x_n = np.linspace(0,2*np.pi,100)
expr = sin(x)
f = lambdify(x, expr, "numpy")
fig, ax = plt.subplots()
ax.plot(x_n, f(x_n))
plt.show()
求解代数方程#
在SymPy中，我们可以使用solveset函数求解代数方程。其。还有一个函数solve也可以用于求解方程，但是并不推荐。
from sympy import *
x = symbols('x')
##注意：以下3种方式都是可以的
print(solveset(Eq(x**2, 1), x))
print(solveset(Eq(x**2 - 1, 0), x))
print(solveset(x**2 - 1, x))
{-1, 1}
{-1, 1}
{-1, 1}
from sympy import *
x = symbols('x')
solveset(x - x, x, domain=S.Reals)
\[\displaystyle \mathbb{R}\]
from sympy import *
x = symbols('x')
solveset(sin(x) - 1, x, domain=S.Reals)
\[\displaystyle \left\{2 n \pi + \frac{\pi}{2}\; \middle|\; n \in \mathbb{Z}\right\}\]
from sympy import *
x = symbols('x')
solveset(exp(x), x)
\[\displaystyle \emptyset\]
from sympy import *
x = symbols('x')
solveset(cos(x) - x, x)
\[\displaystyle \left\{x\; \middle|\; x \in \mathbb{C} \wedge - x + \cos{\left(x \right)} = 0 \right\}\]
from sympy import *
x, y, z = symbols('x y z')
solveset(x + y * z - 1,x)
\[\displaystyle \left\{- y z + 1\right\}\]
对于多元线性方程，我们需要使用linsolve求解
from sympy import *
x, y, z = symbols('x y z')
linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))
\[\displaystyle \left\{\left( - y - 1, \ y, \ 2\right)\right\}\]
from sympy import *
x, y, z = symbols('x y z')
linsolve([x + y + z - 1, x + y + 2*z - 3, x + 2*y + z - 2], (x, y, z))
\[\displaystyle \left\{\left( -2, \ 1, \ 2\right)\right\}\]
对于非线性多元方程的求解，我们需要用到nonlinsolve
from sympy import *
x, y = symbols('x y')
nonlinsolve([x**2 + 1, y**2 + 1], [x, y])
\[\displaystyle \left\{\left( - i, \ - i\right), \left( - i, \ i\right), \left( i, \ - i\right), \left( i, \ i\right)\right\}\]
微分学#
我们可以使用diff函数，求一元函数的导数、多元函数的偏导数。
from sympy import *
x, y, z = symbols('x y z')
diff(cos(x), x)
\[\displaystyle - \sin{\left(x \right)}\]
from sympy import *
x, y, z = symbols('x y z')
diff(exp(x**2), x)
\[\displaystyle 2 x e^{x^{2}}\]
from sympy import *
x, y, z = symbols('x y z')
diff(x**4, x, x, x)
\[\displaystyle 24 x\]
from sympy import *
x, y, z = symbols('x y z')
diff(x**4, x, 3)
\[\displaystyle 24 x\]
from sympy import *
x, y, z = symbols('x y z')
expr = exp(x*y*z)
diff(expr, x, y, y, z, z, z, z)
\[\displaystyle x^{3} y^{2} \left(x^{3} y^{3} z^{3} + 14 x^{2} y^{2} z^{2} + 52 x y z + 48\right) e^{x y z}\]
# diff也可以被作为方法调用
from sympy import *
x, y, z = symbols('x y z')
expr = exp(x*y*z)
expr.diff(x, y, y, z, z, z, z)
\[\displaystyle x^{3} y^{2} \left(x^{3} y^{3} z^{3} + 14 x^{2} y^{2} z^{2} + 52 x y z + 48\right) e^{x y z}\]
from sympy import *
m, n, a, b = symbols('m n a b')
expr = (a*x + b)**m
expr.diff((x, n))
\[\displaystyle \frac{\partial^{n}}{\partial x^{n}} \left(a x + b\right)^{m}\]
积分学#
我们使用integrate来计算不定积分与定积分。
注意： 计算不定积分时，不带常数！在用SymPy做高数作业的时候别把常数漏了
from sympy import *
x = symbols('x')
integrate(cos(x), x)
\[\displaystyle \sin{\left(x \right)}\]
计算定积分时，增加参数(integration_variable, lower_limit, upper_limit)即可，例如，我们想要计算
$\(\int_0^\infty e^{-x}dx\)$
from sympy import *
x = symbols('x')
integrate(exp(-x), (x, 0, oo))
\[\displaystyle 1\]
如果我们想要计算：
\[\int_{-\infty}^\infty\int_{-\infty}^\infty e^{-x^2-y^2}dxdy\]
from sympy import *
x, y = symbols('x y')
integrate(exp(-x**2 - y**2), (x, -oo, oo), (y, -oo, oo))
\[\displaystyle \pi\]
# 如果SymPy遇到无法计算的积分，其会返回一个未计算的积分对象
from sympy import *
x = symbols('x')
integrate(x**x, x)
\[\displaystyle \int x^{x}\, dx\]
极限#
在SymPy中，我们可以使用limit函数计算符号极限。如果我们想要求解\(\lim_{x\to x_0}f(x)\)，我们可以使用limit(f(x), x, x0)
from sympy import *
x = symbols('x')
limit(sin(x)/x, x, 0)
\[\displaystyle 1\]
from sympy import *
x = symbols('x')
expr = x**2/exp(x)
print(expr.subs(x, oo))
print(limit(expr, x, oo))
nan
0
在后面增加参数+或-，可以指定求的是左极限还是右极限
from sympy import *
x = symbols('x')
limit(1/x, x, 0, '+') #左右极限不同
\[\displaystyle \infty\]
from sympy import *
x = symbols('x')
limit(1/x, x, 0, '-') #左右极限不同
\[\displaystyle -\infty\]
求解微分方程#
在SymPy中，我们可以使用dsolve求解微分方程。我们以之前见到过的微分方程为例
$\(m\frac{d^2x}{dt^2}+c\frac{dx}{dt}+kx=0\)$
from sympy import *
f = symbols('f', cls=Function)
x = symbols('x')
m, c, k = symbols('m c k')
diffeq = Eq(m*f(x).diff(x, x) + c*f(x).diff(x) + k*f(x), 0)
diffeq
\[\displaystyle c \frac{d}{d x} f{\left(x \right)} + k f{\left(x \right)} + m \frac{d^{2}}{d x^{2}} f{\left(x \right)} = 0\]
from sympy import *
f = symbols('f', cls=Function)
x = symbols('x')
m, c, k = symbols('m c k')
diffeq = Eq(m*f(x).diff(x, x) + c*f(x).diff(x) + k*f(x), 0)
sol = dsolve(diffeq, f(x), ics={f(0): 5, f(x).diff(x).subs(x, 0): 0})
sol
\[\displaystyle f{\left(x \right)} = \left(- \frac{5 c}{2 \sqrt{c^{2} - 4 k m}} + \frac{5}{2}\right) e^{- \frac{x \left(c + \sqrt{c^{2} - 4 k m}\right)}{2 m}} + \left(\frac{5 c}{2 \sqrt{c^{2} - 4 k m}} + \frac{5}{2}\right) e^{\frac{x \left(- c + \sqrt{c^{2} - 4 k m}\right)}{2 m}}\]
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
f = symbols('f', cls=Function)
x = symbols('x')
m, c, k = symbols('m c k')
diffeq = Eq(m*f(x).diff(x, x) + c*f(x).diff(x) + k*f(x), 0)
sol = dsolve(diffeq, f(x), ics={f(0): 5, f(x).diff(x).subs(x, 0): 0})
a = sol.args[1] #得到等式右边
a = a.subs([(m,1),(c,0.5),(k,4)])
t_eval = np.linspace(0,10,100) #返回解用的时间点
f = lambdify(x, a, "numpy")
t_n = np.linspace(0,10,100)
x_n = f(t_eval).real
#绘制图像
fig, ax = plt.subplots()
ax.plot(t_n, x_n)
plt.show()
