3 使用 SymPy 做数学研究基础

2019-10-27 本文已影响0人水之心

本章导航：

介绍 SymPy 的基础语法。
举了两个实例介绍如何使用 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 不仅仅如此，它还可以支持 L^AT^E_X 公式：

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) \text{d} 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)\text{d} 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 将您的运算结果使用 L^AT^E_X 的形式打印出来：

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 等。但是 x1 与 x2 代表的点的维度必须一样。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))

我们可以计算 ab 与 ac 的方向向量的夹角余弦：

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.p1 与 ab.p2。

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

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

题设：假设有一张图片 $I$ 的尺寸为 $(W_I, H_I)$ ， $I$ 存在一个目标物 $O$ ， $O$ 的尺寸为 $(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 之后的图片

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

3 使用 SymPy 做数学研究基础

9.1 Python 中的数学：符号计算

9.1.1 使用 SymPy 表示数

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

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

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

9.1.2.3 解方程

9.1.3 使用 SymPy 求值

9.1.4 使用 SymPy 做向量运算

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

猜你喜欢

热点阅读