作者 : Fabian Pedregosa


目的

  • 从任意的精度评估表达式。
  • 在符号表达式上进行代数运算。
  • 用符号表达式进行基本的微积分任务 (极限、微分法和积分法)。
  • 求解多项式和超越方程。
  • 求解一些微分方程。

为什么是SymPy? SymPy是符号数学的Python库。它的目的是成为Mathematica或Maple等系统的替代品,同时让代码尽可能简单并且可扩展。SymPy完全是用Python写的,并不需要外部的库。

Sympy文档及库安装见http://www.sympy.org/


章节内容

  • SymPy第一步
    • 使用SymPy作为计算器
    • 练习
    • 符号
  • 代数运算
    • 展开
    • 化简
  • 微积分
    • 极限
    • 微分法
    • 序列扩展
    • 积分法
    • 练习
  • 方程求解
    • 练习
  • 线性代数
    • 矩阵
    • 微分方程

3.2.1 SymPy第一步

3.2.1.1 使用SymPy作为计算器

SymPy定义了三种数字类型:实数、有理数和整数。

有理数类将有理数表征为两个整数对: 分子和分母,因此Rational(1,2)代表1/2, Rational(5,2)代表5/2等等:

In [2]:
from sympy import *
a = Rational(1,2)
In [2]:
a
Out[2]:
1/2
In [3]:
a*2
Out[3]:
1

SymPy在底层使用mpmath, 这使它可以用任意精度的算术进行计算。这样,一些特殊的常数,比如e, pi, oo (无限), 可以被作为符号处理并且可以以任意精度来评估:

In [4]:
pi**2
Out[4]:
pi**2
In [5]:
pi.evalf()
Out[5]:
3.14159265358979
In [6]:
(pi + exp(1)).evalf()
Out[6]:
5.85987448204884

如你所见,将表达式评估为浮点数。

也有一个类代表数学的无限, 称为 oo:

In [7]:
oo > 99999
Out[7]:
True
In [8]:
oo + 1
Out[8]:
oo

3.2.1.2 练习

  • 计算 $\sqrt{2}$ 小数点后一百位。
  • 用有理数算术计算1/2 + 1/3 in rational arithmetic.

3.2.1.3 符号

与其他计算机代数系统不同,在SymPy你需要显性声明符号变量:

In [4]:
from sympy import *
x = Symbol('x')
y = Symbol('y')

然后你可以计算他们:

In [10]:
x + y + x - y
Out[10]:
2*x
In [11]:
(x + y)**2
Out[11]:
(x + y)**2

符号可以使用一些Python操作符操作: +, -, *, ** (算术), &, |, ~ , >>, << (布尔逻辑).


打印 这里我们使用下列设置打印


In [ ]:
sympy.init_printing(use_unicode=False, wrap_line=True)

3.2.2 代数运算

SymPy可以进行强大的代数运算。我们将看一下最常使用的:展开和化简。

3.2.2.1 展开

使用这个模块展开代数表达式。它将试着密集的乘方和相乘:

In [13]:
expand((x + y)**3)
Out[13]:
x**3 + 3*x**2*y + 3*x*y**2 + y**3
In [14]:
3*x*y**2 + 3*y*x**2 + x**3 + y**3
Out[14]:
x**3 + 3*x**2*y + 3*x*y**2 + y**3

可以通过关键词的形式使用更多的选项:

In [15]:
expand(x + y, complex=True)
Out[15]:
re(x) + re(y) + I*im(x) + I*im(y)
In [16]:
I*im(x) + I*im(y) + re(x) + re(y)
Out[16]:
re(x) + re(y) + I*im(x) + I*im(y)
In [17]:
expand(cos(x + y), trig=True)
Out[17]:
-sin(x)*sin(y) + cos(x)*cos(y)
In [18]:
cos(x)*cos(y) - sin(x)*sin(y)
Out[18]:
-sin(x)*sin(y) + cos(x)*cos(y)

3.2.2.2 化简

如果可以将表达式转化为更简单的形式,可以使用化简:

In [19]:
simplify((x + x*y) / x)
Out[19]:
y + 1

化简是一个模糊的术语,更准确的词应该是:powsimp (指数化简)、 trigsimp (三角表达式)、logcombine、radsimp一起。


练习

  • 计算$(x+y)^6$的展开。
  • 化简三角表达式$ \sin(x) / \cos(x)$

3.2.3 微积分

3.2.3.1 极限

在SymPy中使用极限很简单,允许语法limit(function, variable, point), 因此要计算f(x)类似$x \rightarrow 0$, 你应该使用limit(f, x, 0):

In [5]:
limit(sin(x)/x, x, 0)
Out[5]:
1

你也可以计算一下在无限时候的极限:

In [6]:
limit(x, x, oo)
Out[6]:
oo
In [7]:
limit(1/x, x, oo)
Out[7]:
0
In [8]:
limit(x**x, x, 0)
Out[8]:
1

3.2.3.2 微分法

你可以使用diff(func, var)微分任何SymPy表达式。例如:

In [9]:
diff(sin(x), x)
Out[9]:
cos(x)
In [10]:
diff(sin(2*x), x)
Out[10]:
2*cos(2*x)
In [11]:
diff(tan(x), x)
Out[11]:
tan(x)**2 + 1

你可以用下列方法检查是否正确:

In [12]:
limit((tan(x+y) - tan(x))/y, y, 0)
Out[12]:
tan(x)**2 + 1

可以用diff(func, var, n)方法来计算更高的导数:

In [13]:
diff(sin(2*x), x, 1)
Out[13]:
2*cos(2*x)
In [14]:
diff(sin(2*x), x, 2)
Out[14]:
-4*sin(2*x)
In [15]:
diff(sin(2*x), x, 3)
Out[15]:
-8*cos(2*x)

3.2.3.3 序列展开

SymPy也知道如何计算一个表达式在一个点的Taylor序列。使用series(expr, var):

In [16]:
series(cos(x), x)
Out[16]:
1 - x**2/2 + x**4/24 + O(x**6)
In [17]:
series(1/cos(x), x)
Out[17]:
1 + x**2/2 + 5*x**4/24 + O(x**6)

练习

计算$\lim_{x\rightarrow 0} \sin(x)/x$

计算log(x)对于x的导数。


3.2.3.4 积分法

SymPy支持超验基础和特殊函数的无限和有限积分,通过integrate() 功能, 使用了强大的扩展的Risch-Norman算法和启发式和模式匹配。你可以积分基本函数:

In [18]:
integrate(6*x**5, x)
Out[18]:
x**6
In [19]:
integrate(sin(x), x)
Out[19]:
-cos(x)
In [20]:
integrate(log(x), x)
Out[20]:
x*log(x) - x
In [21]:
integrate(2*x + sinh(x), x)
Out[21]:
x**2 + cosh(x)

也可以很简单的处理特殊函数:

In [22]:
integrate(exp(-x**2)*erf(x), x)
Out[22]:
sqrt(pi)*erf(x)**2/4

也可以计算一下有限积分:

In [23]:
integrate(x**3, (x, -1, 1))
Out[23]:
0
In [24]:
integrate(sin(x), (x, 0, pi/2))
Out[24]:
1
In [25]:
integrate(cos(x), (x, -pi/2, pi/2))
Out[25]:
2

不标准积分也支持:

In [26]:
integrate(exp(-x), (x, 0, oo))
Out[26]:
1
In [27]:
integrate(exp(-x**2), (x, -oo, oo))
Out[27]:
sqrt(pi)

3.2.3.5 练习

3.2.4 方程求解

SymPy可以求解线性代数方程,一个或多个变量:

In [28]:
solve(x**4 - 1, x)
Out[28]:
[-1, 1, -I, I]

如你所见,第一个参数是假设等于0的表达式。它可以解一个很大的多项式方程,也可以有能力求解多个方程,可以将各自的多个变量作为元组以第二个参数给出:

In [29]:
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
Out[29]:
{x: -3, y: 1}

也直接求解超越方程(有限的):

In [30]:
solve(exp(x) + 1, x)
Out[30]:
[I*pi]

多项式方程的另一个应用是factorfactor将多项式因式分解为可化简的项,并且可以计算不同域的因式:

In [31]:
f = x**4 - 3*x**2 + 1
factor(f)
Out[31]:
(x**2 - x - 1)*(x**2 + x - 1)
In [32]:
factor(f, modulus=5)
Out[32]:
(x - 2)**2*(x + 2)**2

SymPy也可以解布尔方程,即,判断一个布尔表达式是否满足。对于这个情况,我们可以使用satisfiable函数:

In [33]:
satisfiable(x & y)
Out[33]:
{x: True, y: True}

这告诉我们(x & y)是真,当x和y都是True的时候。如果一个表达式不是True,即它的任何参数值都无法使表达式为真,那么它将返回False:

In [34]:
satisfiable(x & ~x)
Out[34]:
False

3.2.4.1 练习

  • 求解系统方程$x + y = 2$, $2\cdot x + y = 0$
  • 是否存在布尔值,使$(~x | y) & (~y | x)$为真?

3.2.5 线性代数

3.2.5.1 矩阵

矩阵通过Matrix类的一个实例来创建:

In [35]:
from sympy import Matrix
Matrix([[1,0], [0,1]])
Out[35]:
Matrix([
[1, 0],
[0, 1]])

与NumPy数组不同,你也可以在里面放入符号:

In [36]:
x = Symbol('x')
y = Symbol('y')
A = Matrix([[1,x], [y,1]])
A
Out[36]:
Matrix([
[1, x],
[y, 1]])
In [37]:
A**2
Out[37]:
Matrix([
[x*y + 1,     2*x],
[    2*y, x*y + 1]])

3.2.5.2 微分方程

SymPy可以解 (一些) 常规微分。要求解一个微分方程,使用dsolve。首先,通过传递cls=Function来创建一个未定义的符号函数:

In [38]:
f, g = symbols('f g', cls=Function)

f 和 g是未定义函数。我们可以调用f(x), 并且它可以代表未知的函数:

In [39]:
f(x)
Out[39]:
f(x)
In [40]:
f(x).diff(x, x) + f(x)
Out[40]:
f(x) + Derivative(f(x), x, x)
In [41]:
dsolve(f(x).diff(x, x) + f(x), f(x))
Out[41]:
f(x) == C1*sin(x) + C2*cos(x)

关键词参数可以向这个函数传递,以便帮助确认是否找到最适合的解决系统。例如,你知道它是独立的方程,你可以使用关键词hint=’separable’来强制dsolve来将它作为独立方程来求解:

In [42]:
dsolve(sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x), f(x), hint='separable')
Out[42]:
[f(x) == -asin(sqrt(C1/cos(x)**2 + 1)) + pi,
 f(x) == asin(sqrt(C1/cos(x)**2 + 1)) + pi,
 f(x) == -asin(sqrt(C1/cos(x)**2 + 1)),
 f(x) == asin(sqrt(C1/cos(x)**2 + 1))]

练习

  • 求解Bernoulli微分方程

    $x \frac{d f(x)}{x} + f(x) - f(x)^2=0$

  • 使用hint=’Bernoulli’求解相同的公式。可以观察到什么?