只是作业，不考虑太多复杂情况。
线性方程组
基本概念
线性方程组： A X = b AX=b AX=b，
其中， A = [ a i i ] n × n = [ a 11 a 12 . . . a 1 n a 21 a 22 . . . a 2 n . . . . . . . . . . . . a n 1 a n 2 . . . a n n ] A=[a_{ii}]_{n×n}=\begin{bmatrix} a_{11}&a_{12}&...&a_{1n}\\ a_{21}&a_{22}&...&a_{2n}\\...&...&...&...\\a_{n1}&a_{n2}&...&a_{nn}\end{bmatrix} A=[aii]n×n=⎣⎢⎢⎡a11a21...an1a12a22...an2............a1na2n...ann⎦⎥⎥⎤， X = ( x 1 , x 2 , . . . , x n ) X=(x_1,x_2,...,x_n) X=(x1,x2,...,xn)， b = ( b 1 , b 2 , . . . , b n ) b=(b_1,b_2,...,b_n) b=(b1,b2,...,bn)。
3种特殊解
①
A
=
[
a
11
0
.
.
.
0
0
a
22
.
.
.
0
.
.
.
.
.
.
.
.
.
.
.
.
0
0
.
.
.
a
n
n
]
A=\begin{bmatrix} a_{11}&0&...&0\\ 0&a_{22}&...&0\\...&...&...&...\\0&0&...&a_{nn}\end{bmatrix}
A=⎣⎢⎢⎡a110...00a22...0............00...ann⎦⎥⎥⎤时，
x
i
=
b
i
a
i
i
x_i=\displaystyle\frac{b_i}{a_{ii}}
xi=aiibi。②下三角方程组，略，类似上三角。
③上三角方程组，
A
=
[
a
11
a
12
.
.
.
a
1
n
0
a
22
.
.
.
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
0
0
.
.
.
a
n
n
]
A=\begin{bmatrix} a_{11}&a_{12}&...&a_{1n}\\ 0&a_{22}&...&a_{2n}\\...&...&...&...\\0&0&...&a_{nn}\end{bmatrix}
A=⎣⎢⎢⎡a110...0a12a22...0............a1na2n...ann⎦⎥⎥⎤，其中
a
i
i
≠
0
a_{ii}≠0
aii=0，
x
i
=
b
i
−
∑
j
=
i
+
1
n
a
i
j
x
j
a
i
i
x_i=\displaystyle\frac{b_i-\displaystyle\sum_{j=i+1}^{n}{a_{ij}x_j}}{a_{ii}}
xi=aiibi−j=i+1∑naijxj，
i
=
n
,
.
.
.
1
i=n,...1
i=n,...1
上三角方程组的计算：
//利用上三角阵公式进行计算。n从1开始，表示行数（=列数），下同
void CalcUpperTraiangularMatrix(vector<vector<double>>& a, int n, vector<double>& b, vector<double> &x) {
int i, j;
double sum;
//由于索引从0开始，所以i的范围为[n-1,0]。公式中索引从1开始，故其i的范围为[n,1]
//向前不断计算x：
//也就是计算x[n-1]时用到x[n](不存在，也就用不到);
//计算x[n-2]时会用到x[n-1];
//计算x[n-3]时会用到x[n-1]和x[n-2]，不断往前
for (i = n-1; i > -1; i--)
{
sum = b[i];
for (j = i+1; j < n; j++)
{
sum -= a[i][j] * x[j];
}
x[i] = sum / a[i][i];
}
}
高斯消元法
对矩阵A和矩阵b的增广矩阵，转化为上三角阵： ( a 11 a 12 . . . a 1 n b 1 0 a 22 . . . a 2 n b 2 . . . . . . . . . . . . . . . 0 0 . . . a n n b n ) \begin{pmatrix} a_{11}&a_{12}&...&a_{1n}&b_1\\ 0&a_{22}&...&a_{2n}&b_2\\...&...&...&...&...\\0&0&...&a_{nn}&b_{n}\end{pmatrix} ⎝⎜⎜⎛a110...0a12a22...0............a1na2n...annb1b2...bn⎠⎟⎟⎞，再代入公式求解。
1、顺序消元法
对增广矩阵，从第1行开始逐行（也就是顺序）作 行的变换，得到上三角阵。
算法步骤：①消元，转化为上三角阵。②回代，将上三角阵代入公式求解。
//利用高斯顺序消元法，得到上三角阵
void GaussElimination_Sequence(vector<vector<double>> &a, int n, vector<double> &b) {
int i,k,j;
double t;
//k表示用来消元的行号，从第0行到第n-1行
for (k = 0; k <= n - 1; k++)
{
//用第k行，对第i行进行消元，i的范围：[k+1,n-1]
for (i = k+1; i < n; i++)
{
//行变换时，每行乘的非零的系数，得到新行，再用旧行减去它
t = a[i][k] / a[k][k];
//对第i行的每一个元素进行变换。元素的下标范围：[k,n-1]，用j来遍历。
for (j = k; j < n; j++)
{
a[i][j] = a[i][j] - t * a[k][j];
}
//b矩阵在a、b的增广矩阵中，一同作上述的行变换
b[i] = b[i] - t * b[k];
}
}
}
2、列主元消元法
在顺序消元法中 a n n a_{nn} ann位于被除数，所以 a n n a_{nn} ann不能为零。
另外如果 a n n a_{nn} ann很小，计算结果可能引入较大误差。
因此，在进行第k次（从1开始）行变换前，在[k,n]行之间找到列主元绝对值最大的那一行，和第k行交换位置，这样保证 a k k a_{kk} akk始终是该列中绝对值最大的那一个。
比如第一次行变换结束后，得到： ( a 11 a 12 . . . a 1 n b 1 0 a 22 . . . a 2 n b 2 0 a 32 . . . a 3 n b 3 . . . . . . . . . . . . . . . 0 a n 2 . . . a n n b n ) \begin{pmatrix} a_{11}&a_{12}&...&a_{1n}&b_1\\ 0&a_{22}&...&a_{2n}&b_2\\ 0&a_{32}&...&a_{3n}&b_3\\...&...&...&...&...\\0&a_{n2}&...&a_{nn}&b_{n}\end{pmatrix} ⎝⎜⎜⎜⎜⎛a1100...0a12a22a32...an2...............a1na2na3n...annb1b2b3...bn⎠⎟⎟⎟⎟⎞，
此时 a 22 a_{22} a22可能很小，需要在 ( a 22 , a 32 , . . . , a n 2 ) (a_{22},a_{32},...,a_{n2}) (a22,a32,...,an2)之间找到绝对值最大的那一个，比如是 a k 2 a_{k2} ak2，将其对应的第k行和第2行交换位置。
//利用高斯列主元消去法，也就是通过交换行的方法，得到上三角阵
void GaussElimination_ChangeRow(vector<vector<double>>& a, int n, vector<double>& b) {
int i, k, j;
double t;
double maxValue;//列主元最大值
int maxRow;//列主元最大行，也就是需要交换的行
//k表示用来消元的行号，从第0行到第n-1行
for (k = 0; k <= n - 1; k++)
{
//先需要找到列主元最大的行
maxRow = k;
maxValue = abs(a[k][k]);
//遍历每一行的第k个元素，找到最大值及其所在行
for (i = k + 1; i < n; i++)
{
if (abs(a[i][k]) > maxValue) {
maxRow = i;
maxValue = abs(a[i][k]);
}
}
//交换第k行和第maxRow行
//遍历每一个元素，用i下标来遍历
for (i = k; i < n; i++)
{
t = a[k][i];
a[k][i] = a[maxRow][i];
a[maxRow][i] = t;
}
//b矩阵中的对应行也需要交换
t = b[k];
b[k] = b[maxRow];
b[maxRow] = t;
//下面就是正常的高斯顺序消元法了
//用第k行，对第i行进行消元，i的范围：[k+1,n-1]
for (i = k + 1; i < n; i++)
{
//行变换时，每行乘的非零的系数，得到新行，再用旧行减去它
t = a[i][k] / a[k][k];
//对第i行的每一个元素进行变换。元素的下标范围：[k,n-1]，用j来遍历。
for (j = k; j < n; j++)
{
a[i][j] = a[i][j] - t * a[k][j];
}
//b矩阵在a、b的增广矩阵中，一同作上述的行变换
b[i] = b[i] - t * b[k];
}
}
}
迭代法
迭代法的思路：取x=0,y=0作为初值，代入x=f(y),y=g(x)，得到新的x,y的值，再反复代入，不断得到近似值。当结果收敛时，近似值会不断逼近真值。
为了对迭代进行精度控制，引入允许误差ε。当两次迭代结果之差的绝对值小于（等于）ε时，停止迭代，获得近似解。
为了对迭代次数进行控制，避免过多迭代次数甚至死循环，可以设置最大迭代次数。超过最大迭代次数后，终止迭代，返回当前近似值。
1、雅可比迭代法
首先需要满足A矩阵非奇异矩阵（行列式|A|≠0），且 a i i ≠ 0 a_{ii}≠0 aii=0。
迭代序列的分量形式： x i k + 1 = ( b i − ∑ j = 1 , j ≠ i n a i j x j ( k ) ) / a i i ， i = 1 , 2 , . . . , n x_{i}^{k+1}=(b_i-\displaystyle\sum_{j=1,j≠i}^{n}{a_{ij}x_{j}^{(k)}})/a_{ii}，i=1,2,...,n xik+1=(bi−j=1,j=i∑naijxj(k))/aii，i=1,2,...,n
能收敛的一个充要条件是：每行非对角元素的绝对值之和小于对角元素的绝对值，即
∑
j
=
1
,
j
≠
i
n
∣
a
i
j
∣
<
∣
a
i
i
∣
\displaystyle\sum_{j=1,j≠i}^{n}{|a_{ij}|}<|a_{ii}|
j=1,j=i∑n∣aij∣<∣aii∣
//雅可比迭代法，e为允许误差，maxCount为最大迭代次数
void iteration_Jacobi(vector<vector<double>>& a, int n, vector<double>& b, vector<double>& x, double e, int maxCount) {
vector<double> y(n);//每次迭代得到的结果
int i, j, cnt = 0;
double sum, max;
//设置初始值
for (j = 0; j < n; j++)
{
x[j] = 0;
}
do
{
//最大迭代次数限制
if (++cnt > maxCount) break;
//printVector(x, n);
//迭代核心计算
for (i = 0; i < n; i++)
{
sum = 0;
for (j = 0; j < n; j++)
{
if (j!=i)
{
sum += (a[i][j] * x[j]);
}
}
y[i] = (b[i] - sum) / a[i][i];
}
//找到最大误差
max = abs(y[0] - x[0]);
for (j = 1; j < n; j++)
{
if (max < abs(abs(y[j] - x[j]))) {
max = abs(abs(y[j] - x[j]));
}
}
//把计算结果y矩阵赋值给x矩阵，以便下一次迭代用
for (j = 0; j < n; j++)
{
x[j] = y[j];
}
//如果最大误差已经在允许误差范围内，停止迭代
if (max <= e) break;
} while (true);
//cout << "共计算了" << cnt - 1 << "次" << endl;
}
2、高斯-塞德尔迭代法
对雅可比迭代法进行修正。
//高斯塞尔德迭代法
void iteration_Gauss_Seidel(vector<vector<double>>& a, int n, vector<double>& b, vector<double>& x, double e, int maxCount) {
vector<double> y(n);
int i, j, cnt = 0;
double sum1, sum2, max;
//设置初始值
for (j = 0; j < n; j++)
{
x[j] = 0;
}
do
{
//最大迭代次数限制
if (++cnt > maxCount) break;
//printVector(x, n);
//迭代核心计算（与雅可比迭代法不同之处的体现）
for (i = 0; i < n; i++)
{
sum1 = 0;
for (j = 0; j <= i - 1; j++)
{
sum1 += (a[i][j] * y[j]);
}
sum2 = 0;
for (j = i + 1; j < n; j++)
{
if (j != i)
{
sum2 += (a[i][j] * x[j]);
}
}
y[i] = (b[i] - sum1 - sum2) / a[i][i];
}
//找到最大误差
max = abs(y[0] - x[0]);
for (j = 1; j < n; j++)
{
if (max < abs(abs(y[j] - x[j]))) {
max = abs(abs(y[j] - x[j]));
}
}
//把计算结果y矩阵赋值给x矩阵，以便下一次迭代用
for (j = 0; j < n; j++)
{
x[j] = y[j];
}
//如果最大误差已经在允许误差范围内，停止迭代
if (max <= e) break;
} while (true);
//cout << "共计算了" << cnt - 1 << "次" << endl;
}
3、超松弛与亚松弛迭代法
//超松弛与亚松弛迭代法，w为松弛因子。
//当w在(0,1)时，为亚松弛；
//当w>1时，为超松弛，不过一般还要<2；
//当w=1时，为高斯塞尔德迭代
void iteration_SOR(vector<vector<double>>& a, int n, vector<double>& b, vector<double>& x, double e, int maxCount, double w) {
vector<double> y(n);
int i, j, cnt = 0;
double sum1, sum2, max;
//设置初始值
for (j = 0; j < n; j++)
{
x[j] = 0;
}
do
{
//最大迭代次数限制
if (++cnt > maxCount) break;
//printVector(x, n);
//迭代核心计算
for (i = 0; i < n; i++)
{
sum1 = 0;
for (j = 0; j <= i - 1; j++)
{
sum1 += (a[i][j] * y[j]);
}
sum2 = 0;
for (j = i + 1; j < n; j++)
{
if (j != i)
{
sum2 += (a[i][j] * x[j]);
}
}
//与高斯-塞尔德迭代不同之处在于增加了w
y[i] = (1 - w) * x[i] + w * (b[i] - sum1 - sum2) / a[i][i];
}
//找到最大误差
max = abs(y[0] - x[0]);
for (j = 1; j < n; j++)
{
if (max < abs(abs(y[j] - x[j]))) {
max = abs(abs(y[j] - x[j]));
}
}
//把计算结果y矩阵赋值给x矩阵，以便下一次迭代用
for (j = 0; j < n; j++)
{
x[j] = y[j];
}
//如果最大误差已经在允许误差范围内，停止迭代
if (max <= e) break;
} while (true);
//cout << "共计算了" << cnt -1 << "次" << endl;
}
2033
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
