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

目录

• 导言
• 解方程（组）
• • solve函数
  • solveset函数
• 求和 ∑ \sum ∑
• 连乘 ∏ \prod ∏
• 求函数极限
• 求数列极限

导言

在前两篇文章中，我们学习了

SymPy

的输入输出、基本符号、表达式、函数的定义和使用，以及表达式的化简、展开与合并。

传送链接：

「SymPy」符号运算(1) 简介/符号/变量/函数/表达式/等式/不等式/运算符

「SymPy」符号运算(2) 各种形式输出、表达式的化简合并与展开

本篇介绍

SymPy

方程求解，包括：线性/非线性方程求解、线性方程组求解和非线性方程组求解，求解结果分为符号解和数值解1。求解方程由两个主要函数：

solve()

和

solveset()

。

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

解方程（组）

有两个主要的函数都可以对方程进行求根：

solve()

和

solveset()

，对一般的方程均能求解，如

from sympy.abc import x, y
from sympy import solve, solveset

solve(x**2 - y, x, dict=True)
# [{x: -sqrt(y)}, {x: sqrt(y)}]
solveset(x**2 - y, x)
# {-sqrt(y), sqrt(y)}
           

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

1. solve()

   1. solve()

      函数的返回值可以相对容易地进行替换

      subs()

      而参与到其他的运算中
   2. 希望显式地获得满足方程的符号表达式

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=1100​3n

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=1n​2i−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=0n​i,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=1k​i：

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=1k​m

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=1k​i

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
           

1. 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 ↩︎