使用Python求解方程

已剪辑自: https://zhuanlan.zhihu.com/p/24893371

新年第一篇，搞起.

这回写一个好久之前想做，一直搁着没做的东西—— Python 解方程（其实是放假回家，趁着家里电脑重装 LOL 的时间过来写一篇）. 咱这回用三种不同的方法，来应对平常碰到的简单方程.

Numpy 求解线性方程组

例如我们要解一个这样的二元一次方程组：

x + 2y = 3
4x ＋ 5y = 6
1
2

当然我们可以手动写出解析解，然后写一个函数来求解，这实际上只是用 Python 来单纯做“数值计算”. 但实际上，numpy.linalg.solve 可以直接求解线性方程组.

一般地，我们设解线性方程组形如 Ax=b，其中 A 是系数矩阵，b 是一维(n 维也可以，这个下面会提到)，x 是未知变量. 再拿上面地最简单的二元一次方程组为例，我们用 numpy.linalg.solve 可以这样写：

In [1]: import numpy as np
   ...: A = np.mat('1,2; 4,5')    # 构造系数矩阵 A
   ...: b = np.mat('3,6').T       # 构造转置矩阵 b （这里必须为列向量）
   ...: r = np.linalg.solve(A,b)  # 调用 solve 函数求解
   ...: print r
   ...:
Out[1]: [[-1.]
         [ 2.]]
1
2
3
4
5
6
7
8

那么前面提到的“ n 维”情形是什么呢？实际上就是同时求解多组形式相同的二元一次方程组，例如我们想同时求解这样两组：

x + 2y = 3
4x ＋ 5y = 6
1
2

和

x + 2y = 7
4x ＋ 5y = 8
1
2

就可以这样写：

In [2]: import numpy as np
   ...: A = np.mat('1,2; 4,5')          # 构造系数矩阵 A
   ...: b = np.array([[3,6], [7,8]]).T  # 构造转置矩阵 b （这里必须为列向量），
   ...: 注意这里用的是 array
   ...: r = np.linalg.solve(A,b)        # 调用 solve 函数求解
   ...: print r
   ...:
Out[2]: [[-1.         -6.33333333]
         [ 2.          6.66666667]]
1
2
3
4
5
6
7
8
9

SciPy 求解非线性方程组

先看官方文档的介绍：

scipy.optimize.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)[source]
1

一般来说，我们只需要用到 func 和 x0 就够了. func 是自己构造的函数，也就是需要求解的方程组的左端（右端为 0），而 x0 则是给定的初值.

我们来看一个具体的例子，求解：

x + 2y + 3z - 6 = 0
5 * (x ** 2) + 6 * (y ** 2) + 7 * (z ** 2) - 18 = 0
9 * (x ** 3) + 10 * (y ** 3) + 11 * (z ** 3) - 30 = 0
1
2
3

就可以这么写：

In [3]: from scipy.optimize import fsolve
   ...:
   ...: def func(i):
   ...:     x, y, z = i[0], i[1], i[2]
   ...:     return [
   ...:             x + 2 * y + 3 * z - 6,
   ...:             5 * (x ** 2) + 6 * (y ** 2) + 7 * (z ** 2) - 18,
   ...:             9 * (x ** 3) + 10 * (y ** 3) + 11 * (z ** 3) - 30
   ...:            ]
   ...:
   ...: r = fsolve(func,[0, 0, 0])
   ...: print r
   ...:
Out[3]: [ 1.00000001  0.99999998  1.00000001]
1
2
3
4
5
6
7
8
9
10
11
12
13
14

当然，SciPy 也可以用来求解线性方程组，这是因为 scipy.optimize.fsolve 本质上是最小二乘法来逼近真实结果.

SymPy 通吃一切

例如求解一个：

x + 2 * (x ** 2) + 3 * (x ** 3) - 6 = 0
1

直接就是：

In [4]: from sympy import *
   ...: x = symbols('x')
   ...: solve(x + 2 * (x ** 2) + 3 * (x ** 3) - 6, x)
Out[4]: [1, -5/6 - sqrt(47)*I/6, -5/6 + sqrt(47)*I/6]
1
2
3
4

另外，

@Wayne Shi

的这篇 使用 Python 解数学方程 ，就重点讲述了 SymPy 解线性方程组的方法，所以我也就不再赘述了。

其实 SymPy 能干的太多了，有兴趣的可以看一看 GitHub上的 Quick examples.

---

---

SymPy简介

SymPy的官方教程：

https://github.com/sympy/sympy/wiki/Quick-examples

https://docs.sympy.org/latest/tutorial/index.html

已剪辑自: https://blog.csdn.net/starter_____/article/details/81989835

SymPy是符号数学的Python库。它的目标是成为一个全功能的计算机代数系统，同时保持代码简洁、易于理解和扩展。

In [1]:from sympy import *
In [2]:x = Symbol('x')
In [3]:y = Symbol('y')123
1
2
3

展开与折叠

expand( )展开方程

In [8]: ((x+y)**2).expand()
Out[8]: x**2 + 2*x*y + y**212
1
2

facrot( )折叠方程

In [13]: factor(x**2 + 2*x*y + y**2)
Out[13]: (x + y)**212
1
2

分离与合并

apart( )分离整式

In [14]: together(1 + 2/(x - 1))
Out[14]: (x + 1)/(x - 1)12
1
2

together( )合并整式

In [10]: together(1/x+1/y+1/z)
Out[10]: (x*y + x*z + y*z)/(x*y*z)12
1
2

简化表达式

simplify( )普通的化简

In [15]: simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
Out[15]: x - 112
1
2

trigsimp( )三角化简

In [18]: trigsimp(sin(x)/cos(x))
Out[18]: tan(x)12
1
2

powsimp( )指数化简

In [21]: powsimp(x**a*x**b)
Out[21]: x**(a + b)12
1
2

solve( )解方程

第一个参数为要解的方程，要求右端等于0，第二个参数为要解的未知数

一元一次方程

In [7]:solve(x * 3 - 6, x)
[2]12
1
2

二元一次方程

In [8]: solve([2 * x - y - 3, 3 * x + y - 7],[x, y])
Out[8]: {x: 2, y: 1}12
1
2

limit( )求极限

dir=’+’表示求解右极限，dir=’-‘表示求解左极限

In [10]: limit(1/x,x,oo,dir='+')
Out[10]: 0
In [11]: limit(1/x,x,oo,dir='-')
Out[11]: 01234
1
2
3
4

integrate( )求积分

不定积分

In [12]: integrate(sin(x),x)
Out[12]: -cos(x)12
1
2

定积分

In [13]: integrate(sin(x),(x,0,pi/2))
Out[13]: 112
1
2

diff( )求导

In [14]: diff(x**3,x)
Out[14]: 3*x**2

In [15]: diff(x**3,x,2)
Out[15]: 6*x12345
1
2
3
4
5

dsolve( )解微分方程

以 y′=2xy 为例

In [17]: f =Function('f')
In [18]: dsolve(diff(f(x),x) - 2*f(x)*x,f(x))
Out[18]: Eq(f(x), C1*exp(x**2))123
1
2
3

矩阵化简

In [19]:x1,x2,x3 = symbols('x1 x2 x3')
In [20]:a11,a12,a13,a22,a23,a33 = symbols('a11 a12 a13 a22 a23 a33')
In [21]:m = Matrix([[x1, x2, x3]])
In [22]:n = Matrix([[a11, a12, a13], [a12, a22, a23], [a13, a23, a33]])
In [23]:v = Matrix([[x1], [x2], [x3]])
In [24]:f = m * n * v
In [25]:f[0].subs({x1:1, x2:1, x3:1})

Out[26]:Matrix([[x1*(a11*x1 + a12*x2 + a13*x3) + x2*(a12*x1 + a22*x2 + a23*x3) + x3*(a13*x1 + a23*x2 + a33*x3)]])
Out[27]:a11 + 2*a12 + 2*a13 + a22 + 2*a23 + a33
1234567891011
1
2
3
4
5
6
7
8
9
10
11

---

---

已剪辑自: https://www.cnblogs.com/zyg123/p/10549354.html

目录

前言

（一）求解多元一次方程-solve()

1.说明：

解多元一次方程可以使用solve()，在sympy里，等式是用Eq()来表示，

例如：2x=4 表示为：Eq(x*2, 4)

2.源代码：

"""
 解下列二元一次方程
  2x-y=3
  3x+y=7
"""
# 导入模块
from sympy import *

# 将变量符号化
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')

# 解一元一次方程
expr1 = x*2-4
r1 = solve(expr1, x)
r1_eq = solve(Eq(x*2, 4), x)
print("r1:", r1)
print("r1_eq:", r1_eq)

# 解二元一次方程
expr2 = [2*x-y-3, 3*x+y-7]
r2 = solve(expr2, [x, y])
print("r1:", r2)

# 解三元一次方程
f1 = x+y+z-2
f2 = 2*x-y+z+1
f3 = x+2*y+2*z-3
r3 = solve([f1, f2, f3], [x, y, z])
print("r3:", r3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

3.输出：

（二）解线性方程组-linsolve（）

1.说明：

在sympy中，解线性方程组有三种形式：

1. 默认等式为0的形式：linsolve(eq, [x, y, z])
2. 矩阵形式：linsolve(eq, [x, y, z])
3. 增广矩阵形式：linsolve(A，b, x, y, z)

2.源代码：

"""
    x+y+z-2=0
    2x-y+z+1=0
    x+2y+2z-3=0
"""
from sympy import *


x, y, z = symbols("x y z")

# 默认等式为0的形式
print("======默认等式为0的形式 =======")
eq = [x+y+z-2, 2*x-y+z+1, x+2*y+2*z-3]
result = linsolve(eq, [x, y, z])
print(result)
print(latex(result))

# 矩阵形式
print("======矩阵形式 =======")
eq = Matrix(([1, 1, 1, 2], [2, -1, 1, -1], [1, 2, 2, 3]))
result = linsolve(eq, [x, y, z])
print(result)
print(latex(result))

# 增广矩阵形式
print("======增广矩阵形式 =======")
A = Matrix([[1, 1, 1], [2, -1, 1], [1, 2, 2]])
b = Matrix([[2], [-1], [3]])
system = A, b
result = linsolve(system, x, y, z)
print(result)
print(latex(result))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

3.输出：

（三）解非线性方程组-nonlinsolve（）

1.说明：

nonlinsolve()用于求解非线性方程组，例如二次方，三角函数，，，等方程

2.源代码：

"""
    x**2+y**2-2=0
    x**3+y**3=0
"""

import sympy as sy
x, y = sy.symbols("x y")

eq = [x**2+y**3-2, x**3+y**3]
result = sy.nonlinsolve(eq, [x, y])
print(result)
print(sy.latex(result))
1
2
3
4
5
6
7
8
9
10
11
12

3.输出：

{(−1,1),(−1,−12−3i2),(−1,−12+3i2),(1−i,−1+i),(1+i,−1−i),(1−i−63+122−i−23+42,12+32+2−2+32),(1−−12−632+−4−232,−32+12−−8+(−3+1)22),(1−−4−232+−12−632,−32+12+−8+(−3+1)22),(1+−4+232+−12+632,12+32−2−2+32)}

（四）求解微分方程-dsolve（）

1.说明：

求解微分方程使用dsolve（），注意：

f = symbols(‘f’, cls=Function)的作用是声明f()是一个函数。

2.源代码：

from sympy import *


# 初始化
x = symbols('x')
f = symbols('f', cls=Function)

# 表达式
expr1 = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))

# 求解微分方程
r1 = dsolve(expr1, f(x))

print(r1)
print("原式：", latex(expr1))
print("求解后：", latex(r1))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

3.输出：

原式：

f(x)−2ddxf(x)+d2dx2f(x)=sin⁡(x)

解微分后：

f(x)=(C1+C2x)ex+cos⁡(x)2