Python求解线性方程组
文章目录
- solve
- 带状矩阵
- 三角矩阵
- 循环矩阵
scipy.linalg
中封装了一系列用于求解矩阵方程的算法,列表如下
函数 | 方程 | 备注 |
---|---|---|
solve(a, b) | a x = b ax=b ax=b | |
solve_banded((l, u), ab, b) | a x = b ax=b ax=b | a为带状矩阵 |
solveh_banded(ab, b) | a x = b ax=b ax=b | a为厄密正定带状矩阵 |
solve_circulant(c, b) | c x = b cx=b cx=b | c是循环矩阵 |
solve_triangular(a, b) | a x = b ax=b ax=b | a是三角阵 |
solve
solve
是一个功能十分强大的函数,可用于求解形如a@x=b
之类的问题,其参数如下
solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a='gen', transposed=False)
其中assume_a
是对矩阵a
的特征描述,可选四个字符串
gen
普通矩阵(gen
eric)sym
对称矩阵(sym
metric)her
厄密特矩阵(her
mitian)pos
正定矩阵(pos
tive define)
其他参数含义如下
lower
当assume_a
不为gen
时可用,若为True
,则只计算a
的下三角部分,否则只计算上三角部分overwrite_a
为True
时允许覆盖a
的数据overwrite_b
为True
时允许覆盖b
的数据check_finite
为True
时,检查数据是否全为有限值transposed
为True
时,计算a.T@x=b
如果允许覆盖a
或者b
的数据,则可以节省一点内存和运算时间。
import numpy as np
from scipy.linalg import solve
a = np.random.rand(10,10)
x = np.arange(10)
b = a @ x
x1 = solve(a, b)
print(x)
[0 1 2 3 4 5 6 7 8 9]
print(x1)
[-1.28900588e-14 1.00000000e+00 2.00000000e+00 3.00000000e+00
4.00000000e+00 5.00000000e+00 6.00000000e+00 7.00000000e+00
8.00000000e+00 9.00000000e+00]
带状矩阵
如果某个方阵的所有非零元素都集中在对角线附近的带状区域,则为带状矩阵。对于带状矩阵 A A A,若在 i > j + l i>j+l i>j+l时, a i j = 0 a_{ij}=0 aij=0,则 A A A的下带宽(lower bandwidth)为 l l l;若在 j > i + u j>i+u j>i+u时, a i j = 0 a_{ij}=0 aij=0,则 A A A的上带宽(upper bandwidth)为 u u u。例如下面的 5 × 5 5\times5 5×5矩阵中,下带宽 l = 1 l=1 l=1,上带宽 u = 2 u=2 u=2。
[ ⋅ ⋅ ⋅ 0 0 ⋅ ⋅ ⋅ ⋅ 0 0 ⋅ ⋅ ⋅ ⋅ 0 0 ⋅ ⋅ ⋅ 0 0 0 ⋅ ⋅ ] \begin{bmatrix} ·&·&·&0&0\\ ·&·&·&·&0\\ 0&·&·&·&·\\ 0&0&·&·&·\\ 0&0&0&·&·\\ \end{bmatrix} ⋅⋅000⋅⋅⋅00⋅⋅⋅⋅00⋅⋅⋅⋅00⋅⋅⋅
在scipy.linalg
中,solve_banded
和solveh_banded
均用于
a
a
a为带状矩阵时的求解情况,二者均支持check_finite
参数,并且可以选择是否覆盖ab
, b
等,其他参数列表如下
solve_banded((l,u), ab, b,)
solveh_banded(ab, b, lower=False,)
三角矩阵
scipy.linalg
另外提供了对三角矩阵的求解方案solve_triangular
,亦支持check_finite
参数,以及可以选择是否覆盖a
, b
,其他参数列表如下
solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False)
其中,unit_diagonal
为True
时,则假定对角元素为1。
循环矩阵
所谓循环矩阵,其形式为
c = [ a 1 a 2 a 3 ⋯ a n a n a 1 a 3 ⋯ a n − 1 a n − 1 a 1 a 2 ⋯ a n − 2 ⋮ ⋮ ⋮ ⋮ a 2 a 3 a 4 ⋯ a 1 ] c=\begin{bmatrix} a_1&a_2&a_3&\cdots&a_n\\ a_n&a_1&a_3&\cdots&a_{n-1}\\ a_{n-1}&a_1&a_2&\cdots&a_{n-2}\\ \vdots&\vdots&\vdots&&\vdots&\\ a_2&a_3&a_4&\cdots&a_1\\ \end{bmatrix} c= a1anan−1⋮a2a2a1a1⋮a3a3a3a2⋮a4⋯⋯⋯⋯anan−1an−2⋮a1
形如
c
x
=
b
cx=b
cx=b这类求解问题,若
c
c
c是循环矩阵,则可通过快速傅里叶变换进行求解:
x
=
F
−
1
[
F
[
b
]
F
[
c
]
]
x=F^{-1}\big[\frac{F[b]}{F[c]}\big]
x=F−1[F[c]F[b]]。在scipy.linalg
中,对此进行了封装
scipy.linalg.solve_circulant(c, b, singular='raise', tol=None, caxis=-1, baxis=0, outaxis=0)
除了c,b
之外,其他参数含义如下
singularstr
为"raise"
时,则弹出奇异值错误;为lstsq
时,返回最小二乘解。tol
表示精度caxis
当c
的维度大于1时,caxis
表示循环向量所在轴。baxis
当b
的维度大于1时,baxis
表示计算时代求向量所在轴outaxis
若c,b
维度大于1,则返回值也是多维的,此时outaxis
为x
所在轴。