导言

在前两篇文章中，我们学习了SymPy的输入输出、基本符号、表达式、函数的定义和使用，以及表达式的化简、展开与合并。

传送链接：
「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开

本篇介绍SymPy方程求解，包括：线性/非线性方程求解、线性方程组求解和非线性方程组求解，求解结果分为符号解和数值解¹。求解方程由两个主要函数：solve()和solveset()。

此外顺带学习一下求和式、连乘式、函数极限与数列极限的求解。

解方程（组）

有两个主要的函数都可以对方程进行求根：solve()和solveset()，对一般的方程均能求解，如

from sympy.abc import x, y
from sympy import solve, solvesetsolve(x**2 - y, x, dict=True)
# [{x: -sqrt(y)}, {x: sqrt(y)}]
solveset(x**2 - y, x)
# {-sqrt(y), sqrt(y)}

至于二者的区别，官方文档推荐我们根据以下需求进行选择

solve()
1. solve()函数的返回值可以相对容易地进行替换subs()而参与到其他的运算中
2. 希望显式地获得满足方程的符号表达式
solveset()
1. 可以获得较为精确的数学解（通过设置mathematical sets）
2. 可以展现所有的解，甚至无穷多个解
3. 可以自由地限制解的区间

solve函数

solve(f, *symbols, **flags) ，能够求解多项式、超越函数、分段函数、线性方程组等，举例：

x     = sympy.symbols('x')
expr1 = sympy.solve(x - 1, x)                   # 求满足方程 x - 1 = 0 的 x
print(expr1)
# [1]x, y = sympy.symbols('x y')
expr2 = sympy.solve([2*x-y-3,3*x+y-7],[x,y])    # 线性方程组
print(expr2)
# {x: 2, y: 1}

也可以获得符号解，如 y ( x ) = x + 1 x − 1 y(x) = \dfrac{x + 1}{x - 1} y(x)=x−1x+1，求 x ( y ) x(y) x(y)

x, y = symbols('x y')
expr = Eq((x + 1) / (x - 1), y)
print(solve(expr, x))                           # 直接将Eq投喂给solve函数
# [(y + 1)/(y - 1)]

当然*symbols参数也可以使用表达式或函数

solve(x + 2 + sqrt(3), x + 2)
# [-sqrt(3)]
solve(f(x) - x, f(x))
# [x]
solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
#  [x + f(x)]# 使用可索引的符号A[i]
from sympy import Indexed, IndexedBase, Tuple, sqrt
A = IndexedBase('A')
eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
solve(eqs, eqs.atoms(Indexed))
# {A[1]: 1, A[2]: 2}

*symbols参数也可以省略，此时SymPy自动判断哪些是待求解的符号

solve(x - 3)
# [3]
solve(x**2 - y**2)
# [{x: -y}, {x: y}]
solve(z**2*x**2 - z**2*y**2)
# [{x: -y}, {x: y}, {z: 0}]
solve(z**2*x - z**2*y**2)
# [{x: y**2}, {z: 0}]

如果不想得到完全的、显式的解，可以使用参数implicit=True，如求解 x + e x = 0 x+e^x=0 x+ex=0，可以得到解 x = − e x x = -e^x x=−ex

solve(x + exp(x), x)
# [-LambertW(1)]
solve(x + exp(x), x, implicit=True)
# [-exp(x)]

求系数 a , b a, b a,b使得方程 ( a + b ) x − b + 2 = 0 (a+b)x - b + 2 = 0 (a+b)x−b+2=0恒成立

solve((a + b)*x - b + 2, a, b)
# {a: -2, b: 2}

默认情况下，SymPy会在整个复数域进行求解，如果将解限制在实数域内，可以参考如下方法

from sympy import Symbol, solve
x = Symbol('x', real=True)
solve(x**4 - 256, x, dict=True)
[{x: -4}, {x: 4}]                       # 如果不限制real=True，则输出{-4, 4, -4*I, 4*I}

带不等式条件的方程

solve([x < 3, x - 2])
# Eq(x, 2)
solve([x > 3, x - 2])
# False

部分可变参数**flags及缺省值：

dict=False：返回格式为方程的解的映射列表
set=False：返回格式为符号列表与解组成的元组
check=True：对解进行检查，会排除使分式的分母为零的解
simplify=True：是否对解进行化简，设为False可以缩短运行时间
force=False：是否不考虑符号的正负而进行强行化简

solveset函数

函数solveset(f, symbol=None, domain=Complexes) 可以在指定的域内求解等式或不等式，如果在实数域求解，可以用solveset_real，复数域可以用solveset_comples

x = Symbol('x')     # x = Symbol('x', real=True)时，求解结果一样（必须在solveset函数中指定求解域）
pprint(solveset(exp(x) - 1, x), use_unicode=False)
# {2*n*I*pi | n in Integers}R = S.Reals
x = Symbol('x')
solveset(exp(x) - 1, x, R)
# {0}solveset_real(exp(x) - 1, x)
# {0}# 求解不等式
solveset(exp(x) > 1, x, R)
# Interval.open(0, oo)

求解单个方程时，solveset函数返回一个解集（solveset）、有限集（FiniteSet）、区间（Interval）或虚数集（ImageSet），当解不存在时，返回EmpySet空集，如果没在指定的符号限制之外找到了解，返回ConditionSet。

i = Symbol('i', imaginary=True)
solveset(Eq(i**2 + i*sin(i), 1), i, domain=S.Reals)
# ConditionSet(_R, Eq(_R**2 + _R*sin(_R) - 1, 0), Reals)

solveset函数不会特别地给出重根，需要使用roots函数获得多项式的重根

solveset(x**3 - 6*x**2 + 9*x, x)
# {0, 3}
roots(x**3 - 6*x**2 + 9*x, x)
# {0: 1, 3: 2}

不等式只在实数域内求解，如果在复数域内给出不等式，将返回错误NotImplementedError

solveset(exp(x) > 1, x, R)
# Interval.open(0, oo)

求解反函数invert_real(f_x, y, x) （复数域可以用函数invert_complex(f_x, y, x, domain=Complexes)）

invert_real(exp(x), y, x)
# (x, Intersection({log(y)}, Reals))# When does exp(x) == 1?
invert_real(exp(x), 1, x)
# (x, {0})

在指定域内求解方程

from sympy import Interval, pi, sin, solveset
from sympy.abc import x
solveset(sin(x), x, Interval(-pi, pi))
# {0, -pi, pi}

线性方程组可以专门用函数linsolve(system, *symbols)求解

linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))
# {(-y - 1, y, 2)}# 增广矩阵形式
linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))
# {(-y - 1, y, 2)}# Ax=b形式
M = Matrix(((1, 1, 1, 1), (1, 1, 2, 3)))
system = A, b = M[:, :-1], M[:, -1]
linsolve(system, x, y, z)
# {(-y - 1, y, 2)}

非线性方程组可以用专门的函数nonlinsove求解（对于三角函数求根，最好用solve函数）

from sympy import sqrt
system = [x**2 - 2*y**2 -2, x*y - 2]
vars = [x, y]
nonlinsolve(system, vars)
# {(-2, -1), (2, 1), (-√2⋅ⅈ, √2⋅ⅈ), (√2⋅ⅈ, -√2⋅ⅈ)}

求和 ∑ \sum ∑

sympy.summation(f, *symbols, **kwargs) ，用法：

                            b____\   `
summation(f, (i, a, b)) =  )    f/___,i = a

直接上栗子：求 ∑ n = 1 100 3 n \sum_{n=1}^{100} 3n ∑n=11003n

n     = sympy.Symbol('n', integer=True)
expr1 = sympy.summation(3 * n,(n,1,100))
print(expr1)
# 15150

求 ∑ i = 1 n 2 i − 1 \sum_{i=1}^n 2i-1 ∑i=1n2i−1

from sympy import summation, oo, symbols, log
i, n, m = symbols('i n m', integer=True)
summation(2*i - 1, (i, 1, n))
# n**2

求 ∑ i = 0 ∞ ( 1 2 ) i \sum_{i=0}^\infty \left(\frac{1}{2}\right)^i ∑i=0∞(21)i

summation(1/2**i, (i, 0, oo))
# 2

求 ∑ n = 0 m ∑ i = 0 n i , i = 0 , . . . , n , n = 0 , . . . , m \sum_{n=0}^m \sum_{i=0}^n i, \quad i=0, ..., n,\quad n=0, ..., m ∑n=0m∑i=0ni,i=0,...,n,n=0,...,m

summation(i, (i, 0, n), (n, 0, m))
# m**3/6 + m**2/2 + m/3

连乘 ∏ \prod ∏

连乘函数product(*args, **kwargs)，用法：

                             b_____
product(f(n), (i, a, b)) = |   | f(n)|   |i = a

直接上栗子， ∏ i = 1 k i \prod_{i=1}^k i ∏i=1ki：

from sympy import product, symbols
i, n, m, k = symbols('i n m k', integer=True)product(i, (i, 1, k))
# factorial(k)

求 ∏ i = 1 k m \prod_{i=1}^k m ∏i=1km

product(m, (i, 1, k))
# m**k

求 ∏ k = 1 n ∏ i = 1 k i \prod_{k=1}^n \prod_{i=1}^k i ∏k=1n∏i=1ki

product(i, (i, 1, k), (k, 1, n))
# Product(factorial(k), (k, 1, n))

求函数极限

函数为limit(e, z, z0, dir='+') ，返回e(z)函数在z0处的极限。逼近方向dir有+右逼近、-左逼近、+-左右逼近。栗子

from sympy import limit, sin, oo
from sympy.abc import x
limit(sin(x)/x, x, 0)
# 1
limit(1/x, x, 0) # default dir='+'
# oo
limit(1/x, x, 0, dir="-")
# -oo
limit(1/x, x, 0, dir='+-')
# zoo
limit(1/x, x, oo)
# 0

求数列极限

函数为limit_seq(expr, n=None, trials=5)，返回表达式expr在n趋于无穷大时的极限。trials为递归次数，如果返回了None可以尝试提高递归次数获得结果。栗子：

from sympy import limit_seq, Sum, binomial
from sympy.abc import n, k, m
limit_seq((5*n**3 + 3*n**2 + 4) / (3*n**3 + 4*n - 5), n)
# 5/3
limit_seq(binomial(2*n, n) / Sum(binomial(2*k, k), (k, 1, n)), n)        # binomial 二项式
# 3/4
limit_seq(Sum(k**2 * Sum(2**m/m, (m, 1, k)), (k, 1, n)) / (2**n*n), n)   # Sum (有限)求和
# 4

Meurer A, Smith CP, Paprocki M, Čertík O, Kirpichev SB, Rocklin M, Kumar A, Ivanov S, Moore JK, Singh S, Rathnayake T, Vig S, Granger BE, Muller RP, Bonazzi F, Gupta H, Vats S, Johansson F, Pedregosa F, Curry MJ, Terrel AR, Roučka Š, Saboo A, Fernando I, Kulal S, Cimrman R, Scopatz A. (2017) SymPy: symbolic computing in Python. PeerJ Computer Science 3:e103 https://doi.org/10.7717/peerj-cs.103 ↩︎

「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限相关推荐

「SymPy」符号运算(5) Vector向量及运算
目录导言坐标系和向量笛卡尔坐标系与向量坐标原点向量运算四则运算点乘.叉乘并矢/外积(Dyadic) ∇ \nabla ∇算子(Hamiltonian算子) 梯度旋度散度向量操作 ...
「SymPy」符号运算(4) 微积分与有限差分
目录导言积分不定积分定积分多重积分求导一阶导数高阶导数偏导数有限差分常微分差分差分系数高阶差分偏微分差分导言在前几篇中,我们学习了SymPy的基本语法.方程求解等基础知 ...
「SymPy」符号运算(6) 矩阵Matrix及基础运算
目录导言创建矩阵列表初始化行向量列向量维度和数集二元函数 `lambda`函数特殊矩阵基本操作索引增删基础运算向量运算导言在前几篇文章中,我们学习了SymPy基础/高级用 ...
「SymPy」符号运算(7) Matrix进阶运算与线性代数
目录 0 导言 1 矩阵一般属性秩逆迹转置共轭伴随行列式特征值特征向量特征多项式 2 矩阵运算与线性代数范数标准化条件数矩阵分解黑塞矩阵雅克比矩阵 Jordan标准型 ...
「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开
目录导言输出替换.演化化简.合并与展开化简展开合并 `cancel`函数 `apart`函数 `rewrite`函数 `expand_func`函数导言在前一篇文章中,我们简单学习了 ...
「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符
目录 1 简介 2 导入库 3 定义符号 4 定义函数 5 表达式 6 等式/不等式 7 `SymPy`假设与限制 8 运算符/函数常用运算符/函数数学常数三角函数复杂函数指数运算 1 简介 ...
「SymPy」实战之Maxwell分布律分子最概然、均方根与平均速率
目录 0 导言 1 Maxwell速度分布律 2 SymPy推导特征速度 2.1 导入库和方法 2.2 初始化输出 2.3 定义符号 2.4 设置符号替换 2.5 定义分布律 2.6 分布律作图 2. ...
页面元素之「￥」符号的使用原则和技巧
转载自公众号:一鸣说作者:一鸣十二赞产品经理在做页面设计的时候,难免需要展示「金额」.「价格」等跟钱有关的元素.有的时候需要在表示「钱」的数字前放上￥符号,有的时候不放,那么在使用的时候有什 ...
sympy —— Python 符号运算
0. 素因子分解 >> import sympy >> from sympy.ntheory import factorint >> factorint(100) ...

「SymPy」符号运算(3) (非)线性方程(组)求解、数列求和、连乘、求极限

目录

导言