原理
高斯消去法
的基本原理就是用初等变换将用行的,逐次消去未知数的方法,把原来的方程组,化为与其等价的上三角方程组。
设有线性方程
A
x
=
B
\boldsymbol{A}\boldsymbol{x}=\boldsymbol{B}
Ax=B,写成矩阵如下
[
a
11
a
12
⋯
a
1
n
a
21
a
22
⋯
a
2
n
⋮
⋮
⋱
⋮
a
m
1
a
m
2
⋯
a
m
n
]
[
x
1
x
2
⋮
x
n
]
=
[
b
1
b
2
⋮
b
n
]
\left[ \begin{matrix} a_{11}& a_{12}& \cdots& a_{1n}\\ a_{21}& a_{22}& \cdots& a_{2n}\\ \vdots& \vdots& \ddots& \vdots\\ a_{m1}& a_{m2}& \cdots& a_{mn}\\ \end{matrix} \right]\left[ \begin{matrix} x_1\\ x_2\\ \vdots\\ x_n\\ \end{matrix} \right]= \left[ \begin{matrix} b_1\\ b_2\\ \vdots\\ b_n\\ \end{matrix} \right]
⎣⎢⎢⎢⎡a11a21⋮am1a12a22⋮am2⋯⋯⋱⋯a1na2n⋮amn⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤
该线性方程组的初始增广矩阵为
[
a
11
(
0
)
a
12
(
0
)
⋯
a
1
n
(
0
)
b
1
(
0
)
a
21
(
0
)
a
22
(
0
)
⋯
a
2
n
(
0
)
b
2
(
0
)
⋮
⋮
⋱
⋮
⋮
a
n
1
(
0
)
a
n
2
(
0
)
⋯
a
m
n
(
0
)
b
n
(
0
)
]
\left[\begin{matrix} a_{11}^{(0)} & a_{12}^{(0)} & \cdots & a_{1 n}^{(0)} & b_{1}^{(0)} \\ a_{21}^{(0)} & a_{22}^{(0)} & \cdots & a_{2 n}^{(0)} & b_{2}^{(0)} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n 1}^{(0)} & a_{n 2}^{(0)} & \cdots & a_{m n}^{(0)} & b_{n}^{(0)} \end{matrix}\right]
⎣⎢⎢⎢⎢⎡a11(0)a21(0)⋮an1(0)a12(0)a22(0)⋮an2(0)⋯⋯⋱⋯a1n(0)a2n(0)⋮amn(0)b1(0)b2(0)⋮bn(0)⎦⎥⎥⎥⎥⎤
利用初等行变换的方法,即消元,将其化为如下形式
[
1
a
12
(
n
)
⋯
a
1
n
(
n
)
b
1
(
n
)
0
1
⋯
a
2
n
(
n
)
b
2
(
n
)
⋮
⋮
⋱
⋮
⋮
0
0
⋯
1
b
n
(
n
)
]
\left[ \begin{matrix} 1& a_{12}^{\left( n \right)}& \cdots& a_{1n}^{\left( n \right)}& b_{1}^{\left( n \right)}\\ 0& 1& \cdots& a_{2n}^{\left( n \right)}& b_{2}^{\left( n \right)}\\ \vdots& \vdots& \ddots& \vdots& \vdots\\ 0& 0& \cdots& 1& b_{n}^{\left( n \right)}\\ \end{matrix} \right]
⎣⎢⎢⎢⎢⎡10⋮0a12(n)1⋮0⋯⋯⋱⋯a1n(n)a2n(n)⋮1b1(n)b2(n)⋮bn(n)⎦⎥⎥⎥⎥⎤
其上标
n
n
n表示进行第
n
n
n消元后的值
列主元高斯消去法
列主元高斯消去法较高斯消去法多了一步选取该列主元,即选取该列上绝对值最大的元素作为主元的步骤,交换该行与主元所在的行。目的是使得高斯消去法更稳定。举例说明,比如说现在有一个线性方程组
A
x
=
B
\boldsymbol{A}\boldsymbol{x}=\boldsymbol{B}
Ax=B如下
2
x
1
−
x
2
+
3
x
3
=
1
4
x
1
+
2
x
2
+
5
x
3
=
4
x
1
+
2
x
2
=
7
\begin{array}{l} 2 x_{1}-x_{2}+3 x_{3}=1 \\ 4 x_{1}+2 x_{2}+5 x_{3}=4 \\ x_{1}+2 x_{2}=7 \end{array}
2x1−x2+3x3=14x1+2x2+5x3=4x1+2x2=7
其增广矩阵为
[
2
−
1
3
1
4
2
5
4
1
2
0
7
]
\left[\begin{matrix} 2&-1&3&1\\ 4&2&5&4\\ 1&2&0&7\\ \end{matrix}\right]
⎣⎡241−122350147⎦⎤
在进行第一步消元时,先选取第一列中绝对值最大的元素4作为主元,交换第一行和第二行,得到
[
4
2
5
4
2
−
1
3
1
1
2
0
7
]
\left[\begin{matrix} 4&2&5&4\\ 2&-1&3&1\\ 1&2&0&7\\ \end{matrix}\right]
⎣⎡4212−12530417⎦⎤
再进行消元,得到
[
4
2
5
4
0
−
2
0.5
−
3
0
1.5
−
1.25
6
]
\left[\begin{matrix} 4&2&5&4\\ 0&-2&0.5&-3\\ 0&1.5&-1.25&6\\ \end{matrix}\right]
⎣⎡4002−21.550.5−1.254−36⎦⎤
算法
-
∃ l ∈ { k , k + 1 , … , n } a l k = max k ≤ i ≤ n a i k , a l j ⇔ a k j , j = k , k + 1 , … , n \exists l \in \{k,k+1,\ldots,n\}\\a_{lk}=\max\limits_{k\le i \le n}a_{ik},a_{lj}\Leftrightarrow a_{kj},j=k,k+1,\ldots,n ∃l∈{k,k+1,…,n}alk=k≤i≤nmaxaik,alj⇔akj,j=k,k+1,…,n
-
a k j ( k ) = a k j ( k − 1 ) / a k k ( k − 1 ) a i j ( k ) = a i j ( k − 1 ) − a i k ( k − 1 ) a k j ( k ) , i , j = k , k + 1 , … , n a_{k j}^{(k)}=a_{k j}^{(k-1)} / a_{k k}^{(k-1)}\\a_{i j}^{(k)}=a_{i j}^{(k-1)}-a_{i k}^{(k-1)} a_{k j}^{(k)},i,j=k,k+1,\ldots,n akj(k)=akj(k−1)/akk(k−1)aij(k)=aij(k−1)−aik(k−1)akj(k),i,j=k,k+1,…,n
-
x n = a n , n + 1 ( n ) , i = n x i = a i , n + 1 ( n ) − ∑ j = i + 1 n a i j ( n ) x j , i = n − 1 , n − 2 , ⋯ , 1 \begin{array}{l} x_{n}=a_{n, n+1}^{(n)},i=n \\ x_{i}=a_{i, n+1}^{(n)}-\sum_{j=i+1}^{n} a_{i j}^{(n)} x_{j},i=n-1, n-2, \cdots, 1 \\ \end{array} xn=an,n+1(n),i=nxi=ai,n+1(n)−∑j=i+1naij(n)xj,i=n−1,n−2,⋯,1
c程序
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void gauss(double[], double[], double[], int);
int main(int argc, char* argv[])
{
int n = 4;
double a[4 * 4] = { 1, 1, 0, 3,
2, 1, -1, 1,
3, -1, -1, 3,
-1, 2, 3, -1 };
double b[4] = { 4, 1, -3, 4 };
double x[4] = {};
printf("增广矩阵\n");
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
printf("%f\t", a[i * n + j]);
}
printf("%f\n", b[i]);
}
printf("解向量\n");
gauss(a, b, x, n);
for (int i = 0; i < n; i++)
printf("%f\t", x[i]);
printf("\n");
return 0;
}
void gauss(double a[], double b[], double x[], int n)//列主元高斯消去法,a[]系数矩阵,b[]结果向量,x[]解向量,n阶数
{
int i, j, k, exchangeline, exchangeflag = 0;
double temp, max;
for (k = 0; k < n - 1; k++) { //k迭代次数
max = a[k * n + k];
for (i = k + 1; i < n; i++) { //寻找主元,i行号
if (fabs(max) < fabs(a[n * i + k])) {
max = a[n * i + k];
exchangeflag = 1; //交换标志
exchangeline = n * i; //记录需要交换的行号
}
}
if (exchangeflag) { //换行,j列号
for (j = 0; j < n; j++) {
temp = a[exchangeline + j]; //对系数矩阵操作
a[exchangeline + j] = a[k * n + j];
a[k * n + j] = temp;
}
temp = b[exchangeline / n]; //对结果向量操作
b[exchangeline / n] = b[k];
b[k] = temp;
exchangeflag = 0; //清除交换标志
}
for (i = k + 1; i < n; i++) { //消元
temp = a[i * n + k] / a[k * n + k];
b[i] = b[i] - b[k] * temp; //对结果向量操作
for (j = k; j < n; j++)
a[i * n + j] = a[i * n + j] - a[k * n + j] * temp; //对系数矩阵操作
}
}
for (i = n - 1; i > -1; i--) { //回代
temp = b[i];
for (j = n - 1; j > i; j--)
temp = temp - a[n * i + j] * x[j];
x[i] = temp / a[i * n + i];
}
}