高斯消元法
引言
之前看到一篇博客:数列找规律的问题,这篇博客说的是如何用解五元一次方程组的方式来获取数列(长度为5)的拟合曲线。所以想到如何去解一个 n n n元一次方程组,遂有此文。
基本思想
通过一系列的加减消元运算,直到得到类似
k
x
=
b
kx=b
kx=b的式子,然后逐一回代求解
x
x
x向量(用初等行变换把增广矩阵转换为行阶梯阵,然后回代求出方程的解)
将方程组的增广矩阵
A
′
=
(
A
∣
b
)
A'=(A|b)
A′=(A∣b)实施初等行变化为行的最简行,此时该最简行作为增广矩阵对应的方程组与原方程组同解,这样通过解简化的阶梯形矩阵所对应的方程组就求出原方程组的解。
编程步骤
- 消元,将矩阵化简为上三角形式
- 判断解的情况
- 回代
解的判断
-
判断一
无解: 当消元完毕后,发现一行系数都为0,但是常数项不为0,此时无解
多解: 当消元完毕后,发现有多行系数,常数项均为0,此时多解,有几行全为0,就有几个自由元,即变量的值可以任取,有无数种情况可以满足给出的方程组 -
判断二
当 n n n元一次方程组的系数矩阵 A A A可逆时,可求出方程组解 X = A − 1 b X=A^{-1}b X=A−1b,实际上这也是方程组的唯一解。如果方程组系数矩阵 A A A不可逆或者 A A A不是方阵时,
齐次线性方程组的解只有两种情况:唯一解(零解)和无穷多解(非零解)
算法实现
function elimination(mat) {
for (let i = 1,len = mat.length; i < len; i++) {
for (let j = 0; j < i; j++) {
// TODO: 被除数不能为0, 每次运算时,必须保证对角线上的元素不为0(即运算中的分母不为0)
const coe = mat[i][j] / mat[j][j] // 系数
mat[i].forEach((val, idx) => {
mat[i][idx] -= mat[j][idx] * coe
})
}
}
console.log("elimination -> mat\n", mat)
}
/**
* 返回当前方程组的解的状态
* 解的状态:无解,多解,唯一解
* @param {Matrix} mat 处理后的行阶梯阵
*/
function getSolutionStatus(mat) {
let status = 'unique' // 解的状态:无解,多解,唯一解
if (mat.some(row => row.slice(0, -1).every(val => val === 0) && row[row.length - 1] !== 0)) status = 'no'
if (mat.some(row => row.every(val => val === 0))) status = 'multiple'
return status
}
/**
* 回代
* @param {Matrix} mat 处理后的矩阵
*/
function back_substitution(mat) {
const rowLen = mat.length
const columnLen = mat[0].length
const result = Array(rowLen).fill()
for (let t = rowLen - 1; t >= 0; t--) {
const divisor = result.reduceRight((pre, cur, idx) => {
if (cur === undefined) {
return pre
} else {
return pre - mat[t][idx] * cur
}
}, mat[t][columnLen - 1])
result[t] = divisor / mat[t][t]
}
console.log("functiongaussian_elimination -> result\n", result)
return result
}
/**
* 高斯消元法
* @param {Array[[]]} mat 方程组的增广矩阵,默认方程组的系数矩阵都是合法数值
*/
function gaussian_elimination(mat) {
// mat 为增广矩阵
// 校验参数
const rowLen = mat.length
const columnLen = mat[0].length
if (columnLen - rowLen !== 1) {
throw new Error('mat required is the augmented matrix')
}
// 消元
elimination(mat)
// 解的判断
const solutionStatus = getSolutionStatus(mat) // 解的状态:无解,多解,唯一解
if (solutionStatus === 'no') {
return -1
} else if (solutionStatus === 'multiple') {
return []
} else {
// 回代
return back_substitution(mat)
}
}
// 测试数据
const arr = [
[2, 1, -1, 8],
[-3, -1, 2, -11],
[-2, 1, 2, -3]
] // output: [2, 3, -1]
const arr2 = [
[1, 1, 1, 1, 1, 1],
[16, 8, 4, 2, 1, 3],
[81, 27, 9, 3, 1, 5],
[256, 64, 16, 4, 1, 7],
[625, 125, 25, 5, 1, 2020]
] // output: [ 83.79166666667061,-837.9166666667011,2932.7083333334544,-4187.5833333335095,2010.0000000000857 ]
const arr3 = [
[3, 2, 1, 6],
[2, 2, 2, 4],
[4, -2, -2, 2]
] // output: [ 0.9999999999999997, 2.0000000000000004, -1.0000000000000002 ]
const arr4 = [
[47, 28, 19],
[89, 53, 36]
] // output: [ 1, -1 ]
gaussian_elimination(arr5)
顺序消去法虽然编程操作简单,但是存在以下两方面限制:
(1). 每次运算时,必须保证对角线上的元素不为0(即运算中的分母不为0),否则算法无法继续进行。
(2). 即使的值不为零,但如果绝对值很小,由于第
k
k
k次运算中在分母位置,因此作除数会引起很大的误差,从而影响算法的稳定性。
为了减少计算过程中舍入误差对方程组求解的影响,因此是否可以选择绝对值尽可能大的主元作为除数,基于这种思想就有了高斯消去法的改进型:列主元消去法和全主元消去法。
列主元消去法
在第 k k k步消元前,先找出 k k k行下所有第 k k k列元素最大的非零元素 a r , k \large a_{r,k} ar,k,将第 r r r行与第 k k k行进行整行交换,这样既不影响原方程的解,也可以将绝对值最大的 a r , k \large a_{r,k} ar,k作为主元,放在除数的位置上,尽可能减小引入误差。
全主元消去法与列主元消去法类似,只不过列主元消去法是从第
k
k
k列中选取一个最大的元素,与第
k
k
k行进行交换。而全主元消去法是从第k行第k列开始的右下角矩阵中所有元素中选取一个最大的元素作为主元,同时交换
r
r
r行与
c
c
c列,从而保证稳定性。
之前将算法分解为三个步骤:消元,判断解,回代。这样我们只需要写列主元消去法,后面的步骤可以直接调用之前写好的函数,列主元消去法代码如下:
/**
* 列主元消去法
* @param {Matrix} mat 方程组矩阵
*/
function columnElimination(mat) {
const swap = (a, b) => {
const temp = mat[a]
mat[a] = mat[b]
mat[b] = temp
}
const getMaxColumnIndex = (ci) => {
let result = 0
for (let i = 1, len = mat.length; i < len; i++) {
if (Math.abs(mat[i][ci]) > Math.abs(mat[result][ci])) result = i
}
return result
}
for (let i = 1,len = mat.length; i < len; i++) {
for (let j = 0; j < i; j++) {
const maxCI = getMaxColumnIndex(j)
// 矩阵交换两行(列)不需要变号,这和行列式不同。因为:行列式是值,而矩阵是表,每一行可以看做一个记录
if (maxCI !== j) swap(maxCI, j)
const coe = mat[i][j] / mat[j][j] // 系数
mat[i].forEach((val, idx) => {
mat[i][idx] -= mat[j][idx] * coe
})
}
}
console.log("elimination -> mat\n", mat)
}
数列找规律要点
f ( x ) f(x) f(x) 函数中使用的是多项式,如果你想弄得复杂点, f ( x ) f(x) f(x)中也可以使用指数函数、三角函数等等。当然,弄得太复杂的话可能会很难解出方程中常数的值来,如果你只是想凑一个公式,让数列的下一项得到任意指定的值,那么多项式就足够了。
另外,任意长度的数列找规律都可以使用这个方法,原因上面说了,这本质上就是一个给定若干散点求拟合函数的问题,使用拉格朗日插值法等方法,总是能找到符合要求的多项式函数。