hans

hans

【Python】四、Sympy——《用Python做科学计算》


SymPy 是 Python 的数学符号计算库,用它可以进行数学公式的符号推导。为了调用方便,下面所有的实例程序都假设事先从 sympy 库导入了所有内容:

>>> from sympy import *

* 表示所有的意思


1. 经典公式#

1668632086147.jpg
叫做欧拉恒等式,其中 e 是自然指数的底,i 是虚数单位,pi 是圆周率。此公式被誉为数学最奇妙的公式,它将 5 个基本数学常数用加法、乘法和幂运算联系起来。下面用 SymPy 验证一下这个公式。

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

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

欧拉恒等式可以下面的公式进行计算,
1668632127887.jpg
为了用 SymPy 求证上面的公式,我们需要引入变量 x。在 SymPy 中,数学符号是 Symbol 类的对象,因此必须先创建之后才能使用:

>>> x = Symbol('x')

expand 函数可以将公式展开,我们用它来展开 E**(I*x) 试试看:

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

没有成功,只是换了一种写法而已。这里的 exp 不是 math.exp 或者 numpy.exp,而是 sympy.exp,它是一个类,用来表述自然指数函数。


math.exp (x): 返回 x 的指数,x 必须是一个数。
numpy.exp (x): 返回 x 的指数,x 可以是数组,可以是一个数。

>>> import math
>>> import numpy as np

>>> x = 2
>>> y = np.array([1,2,3,4])
>>> z = np.array([2])

>>> math.exp(x)
7.38905609893065
>>> math.exp(z)
7.38905609893065
>>> math.exp(y)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: only length-1 arrays can be converted to Python scalars

>>> np.exp(x)
7.3890560989306504
>>> np.exp(z)
7.3890560989306504
>>> np.exp(y)
array([  2.71828183,   7.3890561 ,  20.08553692,  54.59815003])

numpy 就是 python 的数组运算包!


expand 函数有关键字参数 complex (带有虚部的参数类型,本质就是复数,NumPy 里讲过),当它为 True 时,expand 将把公式分为实数和虚数两个部分:

>>> 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 表示取虚数部分。显然这里的运算将符号 x 当作复数了。为了指定符号 x 必须是实数,我们需要如下重新定义符号 x:

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

python 中单引号 '' 和双引号 "" 用法:
它们其实并没有明确的区别和特别用法,绝大部分情况是可以互相替换的。但是有一种情况需要注意:

>>> str = "Hello,world!"
>>> str1 = 'Hello,world!'
>>> str2 = "Let's go!"
>>> str3 = 'Let\'s go!' # \是转意符号,告诉系统其后'的作用不是引号。
>>> str4 = "Say:\"Hello!\""
>>> str5 = 'Say:"Hello"'
>>> print str
Hello,world!
>>> print str1
Hello,world!
>>> print str2
Let's go!
>>> print str3
Let's go!
>>> print str4
Say:"Hello!"
>>> print str5
Say:"Hello"

下面不加转意符号会出现报错提示:

>>> str6 ='Let's go!'
  File "<stdin>", line 1
    str6 ='Let's go!'
               ^
SyntaxError: invalid syntax

2. 积分运算#

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 坐标函数为:

>>> 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))

很遗憾,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) 计算出。
1668632183810.jpg
因此我们需要对 circle_area 中的变量 r 进行替代:

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

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

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

用 subs 进行算式替换
subs 函数可以将算式中的符号进行替换,它有 3 种调用方式:
expression.subs (x, y) : 将算式中的 x 替换成 y
expression.subs({x,u}) : 使用字典进行多次替换
expression.subs ([(x,y),(u,v)]) : 使用列表进行多次替换
请注意多次替换是顺序执行的,因此:
expression.sub([(x,y),(y,x)])
并不能对两个符号 x,y 进行交换。


我这里举一些例子就明白了。

>>> from sympy import *
>>> a,b,c = symbols('a,b,c') #首字母小写
>>> y = a + 10*b + 100*c
>>> y = y.subs(a,b)
>>> y
11*b + 100*c
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,b),(b,c)])
>>> y
111*c
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,b),(b,a)])
>>> y
11*a + 100*c
>>> x=Symbol('x') #注意这个首字母大写
>>> y = a + 10*b + 100*c
>>> y = y.subs([(a,x),(b,a),(x,b)])
>>> y
10*a + b + 100*c

最后这个按顺序用笔写写。


想看完整内容,请去原文地址: http://old.sebug.net/paper/books/scipydoc/sympy_intro.html

読み込み中...
文章は、創作者によって署名され、ブロックチェーンに安全に保存されています。