从例子开始

在详细介绍SymPy的语法结构和各种运算功能之前,让我们先通过两个实例说明用SymPy解决符号运算问题的一般步骤。

封面上的经典公式

请读者看一下本书封面左上角的那个数学公式:

e^{\mathrm{i} \pi} + 1 = 0

此公式被称为欧拉恒等式,其中e是自然常数,\mathrm{i}是虚数单位, \pi是圆周率。此公式被誉为数学中最奇妙的公式,它将5个基本数学常数用加法、乘法和幂运算联系起来。下面我们用SymPy验证此公式。

从sympy库载入的符号中,E表示自然常数,I表示虚数单位,pi表示圆周率,因此上面的公式可以直接如下计算:

>>> E**(I*pi)+1
0

SymPy除了可以直接计算公式的值,还可以帮助我们做数学公式的推导和证明。欧拉恒等式可以将 \pi代入下面的欧拉公式得到:

e^{{\mathrm{i}}x}=\cos x+ {\mathrm{i}} \sin x

在SymPy中可以使用expand()将表达式展开,我们用它展开e^{{\mathrm{i}}x}试试看:

>>> expand( E**(I*x) )
exp(I*x)

很遗憾没有成功,只是换了一种写法而已。这里的exp是SymPy中表示自然指数函数的类。当expand()的complex参数为True时,表达式将被分为实数和虚数两个部分:

>>> expand(exp(I*x), complex=True)
I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x))

这次将表达式展开了,但是得到的结果相当复杂。其中sin、cos、re、im都是SymPy中定义的表示数学函数的类:re是取实数值的函数,im是取虚数值的函数。显然expand()将x当作复数了。为了指定x为实数,我们需要如下重新定义x

>>> x = Symbol("x", real=True)
>>> expand(exp(I*x), complex=True)
I*sin(x) + cos(x)

终于得到了我们需要的公式,那么如何证明它呢?我们可以用泰勒多项式对其进行展开:

>>> tmp = series(exp(I*x), x, 0, 10)
>>> tmp

(1)1 + i x - \frac{1}{2} x^{2} - \frac{1}{6} i x^{3} + \frac{1}{24} x^{4} + \frac{1}{120} i x^{5} - \frac{1}{720} x^{6} - \frac{1}{5040} i x^{7} + \frac{1}{40320} x^{8} + \frac{1}{362880} i x^{9} + \mathcal{O}\left(x^{10}\right)

将数学表达式转换为LaTex格式

为了让读者更清晰、直观地查看SymPy的计算结果,本书使用数学公式显示较为复杂的结果。幸好NumPy提供了latex()将表达式自动转换为LaTex的数学公式格式,因此这些复杂的数学公式不需要笔者手工输入。例如:

>>> latex(exp(I*x))
$e^{\mathbf{\imath} x}$

series()对表达式进行泰勒级数展开。我们看到展开之后虚数项和实数项交替出现。根据欧拉公式,虚数项的和应该等于\sin x的泰勒级数,而实数项的和应该等于\cos x的泰勒展开。

下面获得tmp的实部:

>>> re(tmp)

(2)1 - \frac{1}{2} x^{2} + \frac{1}{24} x^{4} - \frac{1}{720} x^{6} + \frac{1}{40320} x^{8} + \Re\left(\mathcal{O}\left(x^{10}\right)\right)

下面对\cos x进行泰勒展开,我们看到其各项和上面的结果是一致的。

>>> series(cos(x), x, 0, 10)

(3)1 - \frac{1}{2} x^{2} + \frac{1}{24} x^{4} - \frac{1}{720} x^{6} + \frac{1}{40320} x^{8} + \mathcal{O}\left(x^{10}\right)

下面获得tmp的虚部:

>>> im(tmp)

(4)x  - \frac{1}{6} x^{3} + \frac{1}{120} x^{5} - \frac{1}{5040} x^{7} + \frac{1}{362880} x^{9} + \Im\left(\mathcal{O}\left(x^{10}\right)\right)

下面对\sin x进行泰勒展开,其各项也和上面的结果一致。

>>> series(sin(x), x, 0, 10)

(5)x - \frac{1}{6} x^{3} + \frac{1}{120} x^{5} - \frac{1}{5040} x^{7} + \frac{1}{362880} x^{9} + \mathcal{O}\left(x^{10}\right)

由于e^{{\mathrm{i}}x}的展开式的实部和虚部分别等于\cos x\sin x,因此验证了欧拉公式的正确性。

球体体积

球的体积介绍了如何使用数值定积分计算球体的体积,而SymPy中的integrate()则可以帮助我们进行符号积分。例如下面的语句用integrate()进行不定积分运算:

>>> integrate(x*sin(x), x)
-x*cos(x) + sin(x)

如果指定变量x的取值范围,integrate()则进行定积分运算:

>>> integrate(x*sin(x), (x, 0, 2*pi))
-2*pi

为了计算球体体积,首先看看如何计算圆形面积,假设圆形的半径为r,则圆上任意一点的Y坐标函数为:

y(x) = \sqrt{r^2 - x^2}

因此可以直接对函数y(x)在-r到r区间上进行定积分得到半圆面积:

>>> x, y, r = symbols('x,y,r')
>>> 2 * integrate(sqrt(r*r-x**2), (x, -r, r))
2*Integral((r**2 - x**2)**(1/2), (x, -r, r))

首先需要定义运算中所需的符号,这里用symbols()一次创建多个符号。integrate()没有计算出积分结果,而是直接返回了我们输入的算式。这是因为SymPy不知道r是大于0的,如下重新定义r,就可以得到正确答案了:

>>> r = symbols('r', positive=True)
>>> circle_area = 2 * integrate(sqrt(r**2-x**2), (x, -r, r))
>>> circle_area
pi*r**2

接下来对此面积公式进行定积分,就可以得到球体的体积,但是随着X轴坐标的变化,对应的切面的的半径会发生变化。假设X轴的坐标为x,球体的半径为r,则x处球的切面的半径可以使用前面的公式y(x)计算出。因此需要对圆的面积公式circle_area中的变量r进行替代:

>>> circle_area = circle_area.subs(r, sqrt(r**2-x**2))
>>> circle_area
pi*(r**2 - x**2)

用subs进行算式替换

subs()可以将算式中的符号进行替换,它有3种调用方式:

  • expression.subs(x, y):将算式中的x替换成y
  • expression.subs({x:y,u:v}):使用字典进行多次替换
  • expression.subs([(x,y),(u,v)]):使用列表进行多次替换

请注意多次替换是顺序执行的,因此:

expression.sub([(x,y),(y,x)])

并不能对符号x和y进行交换。

然后对circle_area中的变量x在区间-r到r上进行定积分,就可以得到球体的体积公式:

>>> integrate(circle_area, (x, -r, r))
4*pi*r**3/3

內容目录

上一个主题

SymPy-符号运算好帮手

下一个主题

数学表达式

本页

loading...