Python科学计算利器——SymPy库

本文讲述的核心库:sympy
官方在线文档:http://docs.sympy.org/0.7.1/guide.html#guide

sympy是一个Python的科学计算库,用一套强大的符号计算体系完成诸如多项式求值、求极限、解方程、求积分、微分方程、级数展开、矩阵运算等等计算问题。虽然Matlab的类似科学计算能力也很强大,但是Python以其语法简单、易上手、异常丰富的三方库生态,个人认为可以更优雅地解决日常遇到的各种计算问题。

如果你遇到了一个难题,不要犹豫,来找Python,它几乎不会让你失望的。写本文的初衷也是妹子在做金融作业的时候遇到大量的计算,用普通计算器算真是工程量浩大,所以找到sympy这么个库,写代码帮忙辅助计算一下,然后在这里写博客记录下来。

安装sympy库

pip install sympy

常用的sympy内置符号

虚数单位i

In [13]: import sympy

In [14]: sympy.I
Out[14]: I

In [15]: sympy.I ** 2
Out[15]: -1

# 求-1的平方根
In [16]: sympy.sqrt(-1)
Out[16]: I

注:本文后面的示例都省略导包语句:import sympy

自然对数的底e

In [18]: sympy.E
Out[18]: E

# 求对数
In [20]: sympy.log(sympy.E)
Out[20]: 1

无穷大oo

In [26]: 1/sympy.oo
Out[26]: 0

In [27]: 1 + sympy.oo
Out[27]: oo

圆周率pi

In [60]: sympy.pi
Out[60]: pi

In [61]: sympy.sin(sympy.pi/2)
Out[61]: 1

用sympy进行初等运算

Python 2.x中用除号/做两个整数的除法,实际上是整除运算,为了防止这种情况的发生,避免不必要的麻烦,下文的所有示例一开始都加上一句:from __future__ import division,这个时候除号/本身就变成了真实除法,而//才是整除,比如:

# 导入division包之前
In [1]: 1/2
Out[1]: 0

In [2]: from __future__ import division

# 导入division包之后
In [3]: 1/2
Out[3]: 0.5

In [4]: 1//2
Out[4]: 0

求对数

# 自然对数
In [10]: sympy.log(sympy.E)
Out[10]: 1

In [11]: sympy.log(sympy.E ** 3)
Out[11]: 3

# 以10为底1000的对数
In [12]: sympy.log(1000,10)
Out[12]: 3

求平方根

In [13]: sympy.sqrt(4)
Out[13]: 2

In [14]: sympy.sqrt(-1)
Out[14]: I

求n次方根

# 求8的3次方根
In [15]: sympy.root(8,3)
Out[15]: 2

求k次方

In [21]: 2 ** 3
Out[21]: 8

In [22]: 16 ** (1/2)
Out[22]: 4.0

求阶乘

In [35]:  sympy.factorial(4)
Out[35]: 24

求三角函数

sin函数为例:

In [86]: sympy.sin(sympy.pi)
Out[86]: 0

In [87]: sympy.sin(sympy.pi/2)
Out[87]: 1

表达式与表达式求值

sympy可以用一套符号系统来表示一个表达式,如函数、多项式等,并且可以进行求值,比如:

# 首先定义x为一个符号,表示一个变量
In [96]: x = sympy.Symbol('x')

In [97]: fx = 2*x + 1

# 可以看到fx是一个sympy.core.add.Add类型的对象,也就是一个表达式
In [98]: type(fx)
Out[98]: sympy.core.add.Add

# 用evalf函数,传入变量的值,对表达式进行求值
In [101]: fx.evalf(subs={x:2})
Out[101]: 5.00000000000000

还支持多元表达式:

In [102]: x,y = sympy.symbols('x y')

In [103]: f = 2 * x + y

# 以字典的形式传入多个变量的值
In [104]: f.evalf(subs = {x:1,y:2})
Out[104]: 4.00000000000000

# 如果只传入一个变量的值,则原本输出原来的表达式
In [105]: f.evalf(subs = {x:1})
Out[105]: 2.0*x + y

用sympy解方程(组)

使用sympy.solve函数解方程,该函数通常传入两个参数,第1个参数是方程的表达式(把方程所有的项移到等号的同一边形成的式子),第2个参数是方程中的未知数。函数的返回值是一个列表,代表方程的所有根(可能为复数根)。

解最简单的方程

比如下面我们来求两个方程:

# 首先定义 `x`为一个符号,代表一个未知数
In [24]: x = sympy.Symbol('x')

# 解方程:x - 1 = 0
In [25]: sympy.solve(x - 1,x)
Out[25]: [1]

# 解方程:x ^ 2 - 1 = 0
In [26]: sympy.solve(x ** 2 - 1,x)
Out[26]: [-1, 1]

# 解方程:x ^ 2 + 1 = 0
In [27]: sympy.solve(x ** 2 + 1,x)
Out[27]: [-I, I]

把函数式赋给一个变量

有时候为了书写起来简洁,可以把一个函数式起个名字,比如:

In [30]: x = sympy.Symbol('x')

In [31]: f = x + 1

In [32]: sympy.solve(f,x)
Out[32]: [-1]

解方程组

比如要解这么个二元一次方程组:


代码如下:

# 一次性定义多个符号
In [28]: x,y = sympy.symbols('x y')

In [29]: sympy.solve([x + y - 1,x - y -3],[x,y])
Out[29]: {x: 2, y: -1}

计算求和式

计算求和式可以使用sympy.summation函数,其函数原型为:sympy.summation(f, *symbols, **kwargs)

话不多少,举个栗子,比如求下面这个求和式子的值:


我们用初中的知识可以知道,这个式子的结果为:5050 * 2 = 10100

下面用代码来求:

In [37]: n = sympy.Symbol('n')

In [38]: sympy.summation(2 * n,(n,1,100))
Out[38]: 10100

可见结果是正确的。

如果sympy.summation函数无法计算出具体的结果,那么会返回求和表达式。

解带有求和式的方程

比如求这么一个方程:


代码如下:

In [43]: x = sympy.Symbol('x')

In [44]: i = sympy.Symbol('i',integer = True)

In [46]: f =  sympy.summation(x,(i,1,5)) + 10 * x - 15

In [47]: sympy.solve(f,x)
Out[47]: [1]

求极限

求极限用sympy.limit函数,其函数文档如下:

Signature: sympy.limit(e, z, z0, dir='+')
Docstring:
Compute the limit of e(z) at the point z0.

z0 can be any expression, including oo and -oo.

For dir="+" (default) it calculates the limit from the right
(z->z0+) and for dir="-" the limit from the left (z->z0-).  For infinite
z0 (oo or -oo), the dir argument is determined from the direction
of the infinity (i.e., dir="-" for oo).

函数文档中已经说得很清楚了,下面用代码示例来求几个极限。

如果学过微积分,就会知道微积分中有3个重要的极限:




下面就用sympy.limit函数来分别求这3个极限:

In [53]: x = sympy.Symbol('x')

In [54]: f1 = sympy.sin(x)/x

In [55]: sympy.limit(f1,x,0)
Out[55]: 1

In [56]: f2 = (1+x)**(1/x)

In [57]: sympy.limit(f2,x,0)
Out[57]: E

In [58]: f3 = (1+1/x)**x

In [59]: sympy.limit(f3,x,sympy.oo)
Out[59]: E

可见三个极限的计算结果都完全正确。

求导

求导使用sympy.diff函数,传入2个参数:函数表达式和变量名,举例如下:

In [63]: x = sympy.Symbol('x')

In [64]: f = x ** 2 + 2 * x + 1

In [65]: sympy.diff(f,x)
Out[65]: 2*x + 2

In [66]: f2 = sympy.sin(x)

In [67]: sympy.diff(f2,x)
Out[67]: cos(x)

# 多元函数求偏导
In [68]: y = sympy.Symbol('y')

In [70]: f3 = x**2 + 2*x + y**3

In [71]: sympy.diff(f3,x)
Out[71]: 2*x + 2

In [72]: sympy.diff(f3,y)
Out[72]: 3*y**2

求定积分

使用sympy.integrate函数求定积分,其功能比较复杂,非常强大,下面仅仅举几个比较简单的例子。

先来求一个最简单的积分:


牛顿-莱布尼兹公式可以立马口算出上面这个式子的结果是1,用代码计算如下:

n [74]: x = sympy.Symbol('x')

n [75]: f = 2 * x

# 传入函数表达式和积分变量、积分下限、上限
n [76]: sympy.integrate(f,(x,0,1))
ut[76]: 1

下面来算一个复杂一点的多重积分:



其中:


我们通过口算可以求出f(x)

所以:


下面用代码来计算上述过程:

In [82]: t,x = sympy.symbols('t x')

In [83]: f = 2 * t

In [84]: g = sympy.integrate(f,(t,0,x))

In [85]: sympy.integrate(g,(x,0,3))
Out[85]: 9

求不定积分

同样也是使用sympy.integrate函数求不定积分,下面仅仅举几个比较简单的例子。

比如求下面这个不定积分:


通过观察我们知道它的结果是:


下面用代码来计算这个不定积分的结果:

In [79]: x = sympy.Symbol('x')

In [80]: f = sympy.E ** x + 2 * x

In [81]: sympy.integrate(f,x)
Out[81]: x**2 + exp(x)

总结

从上面的一系列计算可以看出,sympy是个非常强大的科学计算库,本文所讲到的用法仅仅是它强大功能的冰山一角,还需以后在实际使用中进一步发掘。

本文所有较为复杂的数学公式都是先在MS Word的公式编辑器中编辑完之后截图到这里的。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,240评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,328评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,182评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,121评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,135评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,093评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,013评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,854评论 0 273
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,295评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,513评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,678评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,398评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,989评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,636评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,801评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,657评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,558评论 2 352

推荐阅读更多精彩内容

  • Android 自定义View的各种姿势1 Activity的显示之ViewRootImpl详解 Activity...
    passiontim阅读 172,062评论 25 707
  • 作品:罗马记忆房间思维导图 作者:李梓琦 日期:2017年5月13日 谢谢你...
    娜琦琦阅读 299评论 2 4
  • 时间到了, 一切的努力, 变戏法似的, 成为一个谜……
    小剧在成长阅读 206评论 0 2
  • 有一種諾言 是首獨奏的舞曲 像灑了鹽的蜜糖 永遠沾不上 心靈的味蕾 有一種謊言 像把無骨的陽傘 永遠撐不開 被遺棄...
    秋鳶子阅读 259评论 0 1