3 使用 SymPy 做数学研究基础

本章导航:

  1. 介绍 SymPy 的基础语法。
  2. 举了两个实例介绍如何使用 SymPy 做数学研究。

9.1 Python 中的数学:符号计算

Python 有许多优秀的处理数学的工具,本章仅仅介绍符号计算库:SymPy。

9.1.1 使用 SymPy 表示数

Sympy 库完全由 Python 写成,支持符号计算、高精度计算、模式匹配、绘图、解方程、微积分、组合数学、离散 数学、几何学、概率与统计、物理学等方面的功能。

复数一般使用 a + bi 的形式进行表示。其中 a, b \in \mathbb{R},而 i 是虚数的单位。在 SymPy 中针对这些特定的“数”,有专门的符号表示。

使用 I 表示 i,使用 Rational 表示分数。

from sympy import I, Rational, sqrt, pi, E, oo
3 + 4I, Rational(1, 3), sqrt(8), sqrt(-1), pi, E**2, oo

输出的结果是:

3+4i, \frac{1}{3}, 2\sqrt{2}, i, \pi, e^2, \infty

可以看出,SymPy 可以完美的表示我们熟悉的各种数学中定义的“数”。

9.1.2 使用 SymPy 提升您的数学计算能力

为什么要使用 SymPy 呢?直接使用 Python 的 math 库不就可以了吗?您也许会有诸如此类的疑问。在解释原因之前,我们来看一个例子:

import math
math.sqrt(8)

输出的结果是:

2.8284271247461903

2.8284271247461903?您没有看出,由于计算机是以二进制方式进行数学运算的,故而计算虚数的值总是不精确的,而 SymPy 的计算结果便是精确的。

import math, sympy
math.sqrt(8) ** 2, sympy.sqrt(8) ** 2

输出结果是:

(8.000000000000002, 8)

9.1.2.1 对数学表达式进行简化与展开

在 SymPy 中使用 symbols 定义数学中一个名词:变量。使用 symbols 创建多个变量名时,请使用空格或者逗号进行隔开,即 symbols('x y') 或者 symbols('x,y') 的形式。

from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
expr

输出结果:

x + 2y

这样,我们便可以将 expr 当作数学中的表达式了:

expr + 1, expr - x

输出结果:

x+2y+1, 2y

既然谈到数学表达式,总会涉及的其的展开与化简。比如 x(x+2y) = x^2 +2xy,可以这样:

from sympy import expand, factor
expanded_expr = expand(x*expr)
expanded_expr

输出结果是:

x^2 + 2xy

我们还可以将其展开式转换为因子的乘积的形式:

factor(expanded_expr)

输出的结果是:

x(x + 2y)

SymPy 不仅仅如此,它还可以支持 LATEX 公式:

alpha, beta, nu = symbols('alpha beta nu')
alpha, beta, nu

输出结果是:

\alpha, \beta, \nu

需要注意的是 symbols 定义的变量 x, y 等不再是字符串,您将其看作是数学中的变量即可。

9.1.2.2 求导,微分,定积分,不定积分,极限

先载入一些会用到的函数,方法,符号:

from sympy import symbols, sin, cos, exp, oo
from sympy import expand, factor, diff, integrate, limit
x, y, t = symbols('x y t')

在 SymPy 中使用 diff 对数学表达式求导,比如:(\sin(x)e^x)' = e^x \sin(x) + e^x \cos(x),使用 SymPy,则有:

diff(sin(x)*exp(x), x)

输出结果是:

e^x \sin(x) + e^x \cos(x)

在 SymPy 中使用 integrate 求解不定积分与定积分。比如,\int e^x\sin(x) + e^x\cos(x) \text2cccccccc8 x = \sin(x)e^x,即:

integrate(exp(x)*sin(x) + exp(x)*cos(x), x)

输出结果:

\sin(x)e^x

定积分 \int_{-\infty}^{\infty} \sin(x^2)\text2cccccccc8 x = \frac{\sqrt{2 \pi}}{2}

integrate(sin(x**2), (x, -oo, oo))

输出结果:

\frac{\sqrt{2} \sqrt{\pi}}{2}

在 SymPy 中使用 limit 求极限,比如 \lim_{x \to 0} \frac{\sin(x)}{x} = 1,即:

limit(sin(x)/x, x, 0)

输出结果:

1

9.1.2.3 解方程

from sympy import solve, dsolve, Eq, Function

可以使用 SymPy 的 solve 求解 f(x) = 0 这种类型的方程。比如,求解 x^2 - 2 = 0

solve(x**2 - 2, x)

输出结果是:

[-\sqrt{2}, \sqrt{2}]

您还可以使用 SymPy 的 dsolve 函数求解微分方程:

y = Function('y')
dsolve(Eq(y(t).diff(t, t) - y(t), exp(t)), y(t))

输出 y'' - y = e^t 的结果:

C_2 e^{-t} +(C_1 + \frac{t}{2})e^t

小贴士:您可以使用 sympy.latex 将您的运算结果使用 LATEX 的形式打印出来:

from sympy import latex
latex(Integral(cos(x)**2, (x, 0, pi)))

输出结果:

\int\limits_{0}^{\pi} \cos^{2}{\left(x \right)}\, dx

9.1.3 使用 SymPy 求值

在 9.4.2 中我们已经介绍了 SymPy 的大多数符号运算机制,接下来将讨论如何通过“赋值”(数学中的变量赋值)来获取表达式的值。如果您想要为表达式求值,可以使用 subs 函数(Substitution,即置换):

x = symbols('x')
expr = x + 1
expr.subs(x, 2)

输出的结果为:

3

正是我们想要得到的预期结果。

也可以使用 subs 函数作为参数置换,比如:

在 9.1.2.3 中出现的符号 Eq 是用来表示两个数学表达式的相等关系。比如, x + 1 = 4,可以这样:

Eq(x+1, 4)

输出结果:

x+1=4

又比如,(x+1)^2 = x^2 + 2x + 1,可以这样:

Eq((x + 1)**2, x**2 + 2*x + 1)

输出结果:

(x+1)^2 = x^2 + 2x + 1

如果您想要判断两个数学表达式是否相等,您有两种方式。下面我们来判断 \cos^2(x) - \sin^2(x) = \cos(2x)

方式一:使用 equals。

a = cos(x)**2 - sin(x)**2
b = cos(2*x)

a.equals(b)

输出结果为 True 符合预期。

方式二:使用 simplify 化简等式:

simplify(a - b)

输出结果为 0 也符合预期。

当然,也可以作图:

9.1.4 使用 SymPy 做向量运算

在 SymPy 中分别使用 Point,Segment 来表示数学中的线段实体(entity)。

from sympy import Segment, Point
from sympy.abc import x
a = Point(1, 2, 3)
b = Point([2, 3])
c = Point(0, x)
a, b, c

输出结果是:

(Point3D(1, 2, 3), Point2D(2, 3), Point2D(0, x))

Point 用来表示数学中的点(向量):x=(x_1,\cdots,x_n) \in \mathbb{R}^{n},其中 x_j \in \mathbb{R},\, 1\leq j \leq n。对于 Segment(x1,x2)中的参数 x1,x2 可以是 Point 实例,也可以是元组或者列表,array 等。但是 x1x2 代表的点的维度必须一样。Segment(x1,x2) 表示一个有向的线段,即矩形的对角线方向为 x_2 - x_1。比如,定义下图9.4.1的实体:

x1 = Point(1,1)
x2 = Point(2,2)
s = Segment(x1,x2)
图9.1.1 Segment 示意图

图中的有向对角线(向量)可以通过 s.direction 进行计算:

s.direction

输出结果为:

??????????2??(1,1)

即向量 x_2 - x_1。在许多领域中,一般将水平向右作为 X 轴,水平向下作为 Y 轴。左上角作为原点 (0,0)。

有了对角线的方向,便可以计算矩形框的宽与高:

w, h = s.direction

由于 Point 指代数学中的向量,故而其存在加法、减法、数乘以及内积运算,它们均被重载为 Python 的运算符。比如,我们定义向量 a,b,c,且 a+b=c,可有:

a = Point(0,1)
b = Point(1,0)
c = Point(1,1)

下面我们对它们进行加法运算:

a + b == c

输出结果为 True 符合预期。同样有:

a - b, c * 2, a.dot(c)

输出结果为:

(Point2D(-1, 1), Point2D(2, 2), 1)

所有结果均符合预期。下面我们再来看看 Segment

ab = Segment(a, b)
ac = Segment(a, c)
ab.direction, ac.direction

输出的结果是:

(Point2D(1, -1), Point2D(1, 0))

我们可以计算 abac 的方向向量的夹角余弦:

ab.direction.dot(ac.direction) / (ab.length * ac.length)

输出结果为 \frac{\sqrt{2}}{2},即 \cos(\theta) = \frac{a \cdot b}{|a||b|}。很容易求得 \theta = \frac{\pi}{4}。更方便的是,Segment 提供了函数直接计算夹角:

ab.angle_between(ac)

输出结果也是 \frac{\pi}{4}。是不是很方便?如果您想要获取 Segment 实例的全部点可以调用 ab.points,分别调取各点则可以使用 ab.p1ab.p2。

9.2 一个案例:实现图像局部 resize 保持高宽比不变

下面我们利用 9.4 节学习的东西构造一个实现图像局部 resize 保持高宽比不变的模型。

  • 题设:假设有一张图片 I 的尺寸为 (W_I, H_I),I 存在一个目标物 OO 的尺寸为 (W, H)。
  • 目标:将 O resize 为尺寸是 (w,h),而约束条件是保证 resize 之后的边界为 w_m,h_m(注意:这里的 h_m 表示下边界长度)。需要尽可能保持 resize 之后的图片块宽高比不变,即满足:

\frac{W}{H} = \frac{w - 2w_m}{h - h_m - h_p} \;\; (\text{公式 9.5.1})

其中 h_p 指的是上边界的长度。这里只有 h_p 的未知量,我们可以直接求出:

h_p = h - h_m - \frac{W - 2w_m}{W} H \;\; (\text{公式 9.5.2})

为了保证 resize 的一致性,需要在原图补边 H_m, W_m, H_p 对应于 h_m, w_m, h_p。下面给出它们的计算公式:

\frac{w}{h} = \frac{W + 2W_m}{H + H_m + H_p} \;\; (\text{公式 9.5.3})

同时还需要保证 resize 前后的宽高与边界的比例不变:

\frac{W}{W_m} = \frac{w}{w_m} \;\; (\text{公式 9.5.4})\\ \frac{H}{H_m} = \frac{h}{h_m} \;\; (\text{公式 9.5.5})

联合上述公式,便可以求出所有未知量。至此,模型建立完毕。下面考虑使用 Python 实现。

from sympy import Point, Segment, symbols, Eq, solve
from dataclasses import dataclass

@dataclass
class Resize:
    W: int
    H: int
    w: int
    h: int
    w_m: int
    h_m: int

    def model(self):
        H_m, W_m, H_p, h_p = symbols("H_m, W_m, H_p, h_p")
        eq1 = Eq(self.W/self.H, (self.w - 2 * self.w_m)/(self.h - self.h_m - h_p))
        eq2 = Eq(self.w/self.h, (self.W + 2 * W_m)/(self.H + H_m + H_p))
        eq3 = Eq(self.W/W_m, self.w/self.w_m)
        eq4 = Eq(self.H/H_m, self.h/self.h_m)
        R =  solve([eq1, eq2, eq3, eq4], [H_m, W_m, H_p, h_p])
        return R

# 传入已知量
W, H = 100, 100
w, h = 32, 32
w_m, h_m = 5, 5
# 模型建立
r = Resize(W, H, w, h, w_m, h_m)
# 模型求解
R = r.model()
# 查看求解结果
R

输出结果为:

{H_m: 15.6250000000000,
 W_m: 15.6250000000000,
 H_p: 15.6250000000000,
 h_p: 5.00000000000000}

这样我们完成了设定的目标。为了更直观,从网上找一张猫的图片作为例子。为了可视化的方便,我们需要定义一个画边界框的函数:

def bbox_to_rect(bbox, color):
    # 将边界框(左上x, 左上y, 右下x, 右下y)格式转换成matplotlib格式:
    # ((左上x, 左上y), 宽, 高)
    return plt.Rectangle(
        xy=(bbox[0], bbox[1]), width=bbox[2]-bbox[0], height=bbox[3]-bbox[1],
        fill=False, edgecolor=color, linewidth=1)

该函数的参数 bbox 取值为 (x_{\text{左上}}, y_{\text{左上}}, x_{\text{右下}}, y_{\text{右下}}),color 指定框的颜色。

from matplotlib import pyplot as plt
%matplotlib inline
# 读取图片
img = plt.imread("cat.jpg")
bbox = [65, 25, 235, 175] # 框的坐标
fig = plt.imshow(img)
fig.axes.add_patch(bbox_to_rect(bbox, 'blue'))
plt.show()

输出结果为:


图9.2.1 一张猫的图片

数组 img 指代 I,而图中的“猫”(即 O)我们可以使用切片的方式获取,即 img[bbox[1]:bbox[3]+1,bbox[0]:bbox[2]+1](这里需要注意,切片的第一维度为高),即:

patch = img[bbox[1]:bbox[3]+1,bbox[0]:bbox[2]+1]
plt.imshow(patch)
plt.show()

输出结果为:

图9.2.2 猫的切片,简称**猫片**,即图像块

我们可以直接求得猫片的宽和高:h, w = patch.shape[:2],但是,patch 是由 bbox 获取的,故而,我们也可以直接求得:

def get_aspect(bbox):
    '''
    计算 bbox 的宽和高
    '''
    # 加 1 是因为像素点是的尺寸为 1x1
    w = bbox[2] - bbox[0] + 1
    h = bbox[3] - bbox[1] + 1
    return w, h
# 获取宽高
W, H = get_aspect(bbox)

我们想要将猫片 resize 为 (480, 480),则可以设定:

w, h = 480, 480
w_m, h_m = 5, 5

接着,计算边界:

r = Resize(W, H, w, h, w_m, h_m)
R = r.model()
H_m, W_m, H_p, h_p = R.values()

由于这里获取的值是浮点数,我们需要将其转换为整数:

def float2int(*args):
    return [int(arg)+1 for arg in args]
# 获取整数型数据
H_m, W_m, H_p, h_p = float2int(*R.values())

最后,我们原图与 resize 之后的图片展示如下:

import cv2
patch = img[bbox[1]-H_p:bbox[3]+1+H_m,bbox[0]-W_m:bbox[2]+1+W_m]
resize = cv2.resize(patch, (w, h))
plt.imshow(patch)
plt.title("origin")
plt.show()
plt.imshow(resize)
plt.title("resize")
plt.show()

输出结果为:


图9.2.3 原图
图9.2.4 resize 之后的图片

本章的代码逻辑并不是很严谨,主要是提供如何一个创建计算机视觉软件的思路,您仔细研究本章的软件设计思路将会收获很多的。

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

推荐阅读更多精彩内容