本系列内容主要针对数值计算方法
第一部分 引言(误差、有效数字、减小误差的方法)
第二部分 非线性方程解法(二分法、牛顿迭代法等)
第三部分 线性方程解法(直接法和迭代法)
第四部分 差值与曲线拟合(拉格朗日、牛顿插值法、最小二乘法等)
第五部分 数值积分与数值微分(牛顿-柯特斯公式,复化积分法等)
第六部分 常微分方程的初值问题的解法
线性方程解法一般分为直接法和迭代法,这次内容我们主要讲述线性方程的直接法,本文主要是简述多种直接法求解,内容包含理论、代码、以及例题分析.
一、高斯Gauss消去法
在线性代数中我们曾经学习过克莱姆Cramer法则,这是一种不实用的直接法,本章介绍几种实用的直接法。高斯Gauss消去法是一种规则化的加减消元法,其基本思想是通过逐次消元计算,把一般线性方程组的求解转化为等价的上三角形方程组的求解。
1.顺序高斯消去法
先来一道简单的三元一次线性方程组:
那么我们应该如何用代码实现呢?对于这个题目呢,我们首先要录入线性方程的增广矩阵
[2 4 -2 | 2]
[1 -3 -3 | -1]
[4 2 2 | 3]
首先我们需要将列出的矩阵导入到代码当中
double A[N][N+1] = {{2, 4, -2, 2}, {1, -3, -3, -1}, {4, 2, 2, 3}};
然后构造函数计算方法来对矩阵进行处理,具体的完整代码如下:
#include <stdio.h>
#define N 3 // 线性方程组的阶数
void gaussian_elimination(double A[N][N+1]) {
int i, j, k;
double factor, temp;
// 列主元消去
for (i = 0; i < N - 1; i++) {
for (j = i + 1; j < N; j++) {
factor = A[j][i] / A[i][i];
for (k = i; k <= N; k++) {
A[j][k] -= factor * A[i][k];
}
}
}
// 回代求解
for (i = N - 1; i >= 0; i--) {
temp = A[i][N];
for (j = i + 1; j < N; j++) {
temp -= A[i][j] * A[j][N];
}
A[i][N] = temp / A[i][i];
}
}
int main() {
double A[N][N+1] = {{2, 4, -2, 2}, {1, -3, -3, -1}, {4, 2, 2, 3}};
gaussian_elimination(A);
printf("Solution:\n");
for (int i = 0; i < N; i++) {
printf("x%d = %.2f\n", i+1, A[i][N]);
}
return 0;
}
运行上面代码,你会得到Solution:x1 = 0.50 x2 = 0.33 x3 = 0.17
即(x1 = 1/2 x2 = 1/3 x3 = 1/6)
2.列主元素高斯消元法
在高斯消去法中,消元过程中可能出现某个元素接近无穷小的情况,这时消去法将无法进行;即使主元素≠0但很小,用其作除数,也会导致其他元素数量级的严重增长和舍入误差的扩散([数值计算方法] 零基础的学习资料和总结 ( c语言描述版本))关于误差我在这篇文章中有详细的描述。为避免此种情况的发生,可通过交换方程的次序,选取绝对值大的元素作为主元。基于这种想法导出了主元素法。
步骤:
- 构造增广矩阵,即系数矩阵A和常数向量b组合在一起。
- 对于矩阵的第k列,找到绝对值最大的元素ar,k,并将该行与第k行交换。
- 使用这个主元进行消元,使得第k列的第k+1行到底部的元素都变为零。
- 重复步骤2和3,直到所有的主对角线上的元素都成为主元,且其下方的元素都是零,此时矩阵成为上三角形式。
- 使用回代算法求解上三角形式,得到方程组的解.
例题: 现在让我们看这么一道题
2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3
使用列主元素高斯消元法,我们首先找到第一列绝对值最大的元素,即-3
,然后交换第一行和第二行。接下来,我们对第二列进行消元,使得第二列的第二行和第三行元素变为零。重复这个过程,直到所有列都被消元完毕。最后,通过回代求解出每个变量的值。如下就是代码的实现方式,根据大家的需求进行修改就可以解不同的题目啦!
#include <stdio.h>
#include <math.h>
#define N 3 // 定义矩阵的维度
void pivotGaussianElimination(double A[N][N+1], double x[N]) {
int i, j, k, r;
double max, temp, sum;
// 遍历每一列,进行列主元素高斯消元法
for (k = 0; k < N - 1; k++) {
max = fabs(A[k][k]); // 初始化主元素的绝对值为当前列的第一个元素的绝对值
r = k; // 记录当前主元素所在的行
// 在当前列中找到绝对值最大的元素,确定主元素
for (i = k + 1; i < N; i++) {
if (fabs(A[i][k]) > max) {
max = fabs(A[i][k]);
r = i;
}
}
// 将主元素所在的行与当前行交换
for (j = k; j <= N; j++) {
temp = A[k][j];
A[k][j] = A[r][j];
A[r][j] = temp;
}
// 将主元素下方的元素消为零
for (i = k + 1; i < N; i++) {
A[i][k] /= A[k][k];
for (j = k + 1; j <= N; j++) {
A[i][j] -= A[k][j] * A[i][k];
}
}
}
// 回代求解未知变量
x[N - 1] = A[N - 1][N] / A[N - 1][N - 1];
for (i = N - 2; i >= 0; i--) {
sum = 0.0;
for (j = i + 1; j < N; j++) {
sum += A[i][j] * x[j];
}
x[i] = (A[i][N] - sum) / A[i][i];
}
}
int main() {
double A[N][N+1] = {
{2, 1, -1, 8},
{-3, -1, 2, -11},
{-2, 1, 2, -3}
};
double x[N];
pivotGaussianElimination(A, x);
printf("Solution:\n");
for (int i = 0; i < N; i++) {
printf("x%d = %f\n", i + 1, x[i]);
}
return 0;
}
列主元Gauss消去法是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素。而全主元Gauss消去法是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增,实际应用中并不经常使用。
3.高斯-若尔当消去法
高斯-若尔当消去法是高斯消去法的扩展,它不仅将线性方程组转换为上三角形式,而且进一步通过行操作将其转换为简化行阶梯形(Reduced Row Echelon Form, RREF),即对角线上的元素都是1,而对角线以上的元素都是0的矩阵。这种方法可以精确地解决线性方程组,而不受舍入误差的影响。
步骤:
- 构造增广矩阵,将系数矩阵A和常数向量b组合在一起。
- 对于矩阵的第k列,找到绝对值最大的元素ar,k,并将该行与第k行交换。
- 使用这个主元进行消元,不仅使第k列的第k+1行到底部的元素都变为零,而且还要使第k行的主元变为1。
- 对于第k行,除以主元使它变为1,然后使用这一行去消去第k+1行到第n行的第k列元素。
- 重复步骤2到4,直到所有的对角线上的元素都成为1,且对角线以上的元素都是零。
例题:
2x - 3y + z = 4
-6x + 9y - 3z = -12
x - y + 0z = 1
首先,我们构造增广矩阵:
| 2 -3 1 | 4 |
| -6 9 -3 | -12|
| 1 -1 0 | 1 |
通过行交换和消元操作,我们逐步将矩阵转换为简化行阶梯形:
| 1 -1 0 | 1 |
| 0 3 -1 | 2 |
| 0 0 1 | 4 |
观察上面矩阵的第三行可以容易的得出Z=4,带入第二行得Y=2,最后计算第一行得出X=1。
下面我们给出具体代码(我已经在关键的步骤上写了中文注释):
#include <stdio.h>
#include <math.h>
#define MAX_SIZE 100
void gaussianJordanElimination(double A[][MAX_SIZE], int n, double x[]) {
int i, j, k, max_row;
double max_val, temp;
// 遍历每一列进行消元
for (k = 0; k < n - 1; k++) {
// 找到主元素所在的行并交换到当前行
max_val = fabs(A[k][k]);
max_row = k;
for (i = k + 1; i < n; i++) {
if (fabs(A[i][k]) > max_val) {
max_val = fabs(A[i][k]);
max_row = i;
}
}
// 交换行 k 和 max_row
for (j = 0; j <= n; j++) {
temp = A[k][j];
A[k][j] = A[max_row][j];
A[max_row][j] = temp;
}
// 将对角线元素化为1,将对角线上方和下方的元素化为0
for (i = 0; i < n; i++) {
if (i != k) {
A[i][k] -= A[k][k] * A[i][k];
}
}
A[k][k] = 1.0 / A[k][k]; // 对角线元素归一化为1
for (j = k + 1; j <= n; j--) {
A[k][j] *= A[k][k];
}
}
// 最后一行自动成为单位行,直接将解复制过来
for (i = 0; i < n; i++) {
x[i] = A[n - 1][i];
}
}
int main() {
int n;
printf("输入线性方程组的阶数: ");
scanf("%d", &n);
double A[MAX_SIZE][MAX_SIZE], x[n];
// 输入增广矩阵
printf("输入增广矩阵:\n");
for (int i = 0; i < n; i++) {
for (int j = 0; j <= n; j++) {
scanf("%lf", &A[i][j]);
}
}
gaussianJordanElimination(A, n, x);
// 输出解向量
printf("解向量为:\n");
for (int i = 0; i < n; i++) {
printf("x%d = %f\n", i + 1, x[i]);
}
return 0;
}
二、矩阵三角分解法
1.Gauss消去法的矩阵表示
步骤:
- 初始化一个单位矩阵作为下三角矩阵的初始状态。
- 对于每一步骤的高斯消去,确定需要的行变换(例如,交换行、行乘以标量、行相加)。
- 构建一个初等矩阵来表示该步的行变换。
- 将这个初等矩阵乘以当前的下三角矩阵,得到新的下三角矩阵。
- 应用相同的行变换到原始矩阵上,得到上三角矩阵。
- 不断循环步骤2到5,直到原始矩阵完全转换为上三角形式。
例题: 接下来我们增加一下例题难度,来看下面例题:
3x + 4y - 5z + 2w = 1
x - 2y + 3z - 4w = 2
-6x + 8y - 3z + 7w = 3
9x - 4y + z - 2w = 4
让我们分步看看这道题:首先依旧是把系数提取出来构造增广矩阵
int main() {
double A[N][N+1] = {
{3, 4, -5, 2, 1},
{1, -2, 3, -4, 2},
{-6, 8, -3, 7, 3},
{9, -4, 1, -2, 4}
};
接下来构造函数,在每一步中,我们寻找当前列绝对值最大的元素作为主元,并将其所在的行与该步的初始行交换,以确保每一步的消元过程都是有效的。利用主元,我们通过行变换将当前列的其他行中的对应元素消为零。当矩阵转换为行阶梯形后,我们从最后一行开始,逐步回代求解每个变量的值。最终我们可以构造出下面完整代码。
#include <stdio.h>
#define N 4
// 打印矩阵
void printMatrix(double A[N][N+1]) {
for (int i = 0; i < N; i++) {
for (int j = 0; j < N+1; j++) {
printf("%f ", A[i][j]);
}
printf("\n");
}
}
// 高斯消去法求解线性方程组
void gaussianElimination(double A[N][N+1]) {
int i, j, k, max_row;
double max_val, multiplier;
// 进行消元操作,将矩阵转化为上三角形式
for (k = 0; k < N - 1; k++) {
max_val = fabs(A[k][k]);
max_row = k;
// 找到最大元素所在的行,进行主元素交换
for (i = k + 1; i < N; i++) {
if (fabs(A[i][k]) > max_val) {
max_val = fabs(A[i][k]);
max_row = i;
}
}
// 交换最大元素所在的行和当前行
for (j = k; j <= N; j++) {
double temp = A[max_row][j];
A[max_row][j] = A[k][j];
A[k][j] = temp;
}
// 使用当前行将下方元素消为0
for (i = k + 1; i < N; i++) {
multiplier = A[i][k] / A[k][k];
// 执行行操作:A[i] = A[i] - multiplier * A[k]
for (j = k; j <= N; j++) {
A[i][j] -= multiplier * A[k][j];
}
}
}
// 矩阵现在已经处于行阶梯形式。进行回代以找到解
for (i = N - 1; i >= 0; i--) {
for (j = i + 1; j < N; j++) {
A[i][N] -= A[i][j] * A[j][N];
}
A[i][N] /= A[i][i];
}
}
int main() {
double A[N][N+1] = {
{3, 4, -5, 2, 1},
{1, -2, 3, -4, 2},
{-6, 8, -3, 7, 3},
{9, -4, 1, -2, 4}
};
printf("原始增广矩阵:\n");
printMatrix(A);
gaussianElimination(A);
printf("行阶梯形式矩阵:\n");
printMatrix(A);
return 0;
}
2.杜里特尔 Doolittle分解法
步骤:
- 首先,我们定义并输入系数矩阵A。
- 通过Doolittle分解法,我们将矩阵A分解为下三角矩阵L和上三角矩阵U
- 将得到的下三角矩阵L和上三角矩阵U回代求解
例题:依旧是一小道
x + 3y + 7z = 6
4x + 2y + 6z = 13
5x + 9y + z = 14
根据我给出的步骤也不难得出下面代码:
#include <stdio.h>
#define N 3
// Doolittle分解法进行LU分解
void DoolittleDecomposition(double A[N][N], double L[N][N], double U[N][N]) {
// 初始化L矩阵为对角线为1,其余为0
for (int i = 0; i < N; i++) {
// 对角线元素为1
for (int j = 0; j < (i + 1); j++) {
L[i][j] = 0;
}
L[i][i] = 1; // 对角线元素为1
// 初始化U矩阵为0
for (int j = i; j < N; j++) {
U[i][j] = 0;
}
}
// 执行分解
for (int j = 0; j < N; j++) {
for (int i = j; i < N; i++) {
// 计算U矩阵的元素
if (j == i) {
U[i][j] = A[i][j]; // 主对角线元素直接赋值
} else {
double sum = 0;
// 计算累加和
for (int k = 0; k < j; k++) {
sum += L[i][k] * U[k][j];
}
U[i][j] = A[i][j] - sum; // 计算U矩阵的非对角线元素
}
}
}
}
// 使用Doolittle分解得到的L和U矩阵进行回代求解方程组
void solveEquations(double L[N][N], double U[N][N], double b[N], double x[N]) {
double y[N];
// 前向替换求解Ly = b
for (int i = 0; i < N; i++) {
double sum = 0;
for (int j = 0; j < i; j++) {
sum += L[i][j] * y[j];
}
y[i] = b[i] - sum;
}
// 回代求解Ux = y
for (int i = N - 1; i >= 0; i--) {
double sum = 0;
for (int j = i + 1; j < N; j++) {
sum += U[i][j] * x[j];
}
x[i] = (y[i] - sum) / U[i][i];
}
}
int main() {
double A[N][N] = {
{1, 3, 7},
{4, 2, 6},
{5, 9, 1}
};
double L[N][N], U[N][N], b[N], x[N];
// 设置右端向量b
b[0] = 1;
b[1] = 2;
b[2] = 3;
// 进行Doolittle分解
DoolittleDecomposition(A, L, U);
// 求解方程组Ax = b
solveEquations(L, U, b, x);
// 打印解向量x
printf("解向量x:\n");
for (int i = 0; i < N; i++) {
printf("%f\n", x[i]);
}
return 0;
}
以上就是本次的全部内容了,希望给大家带来指导和帮助。
版权声明:本文为博主作者:OnlyOswald原创文章,版权归属原作者,如果侵权,请联系我们删除!
原文链接:https://blog.csdn.net/m0_74431077/article/details/138195989