将矩阵 A 通过一系列的线性变换转换为三角矩阵(上三角:upper triangular,下三角:lower triangular),然后便可轻松的求解联立线性方程组。这也是诸如LU分解存在的理由。
0、 思路
这里,我们将会看到所谓思路(人脑的思维方式)和计算机的执行方式,并不完全一致。
我们对
我们以如下的方程组为例:
其对应的计算流程应当是:
- x[4] = b[4]/A[4][4]
- x[3] = (b[3] - A[3][4]*x[4])/A[3][3]
- x[2] = (b[2] - A[2][4]*x[4] - A[2][3]*x[3])/A[2][2]
- x[1] = (b[1] - A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]) / A[1][1]
一、 matlab解法
clear all; close all; clc
A = [1, 2, 3, 4; 0, 5, 6, 7; 0, 0, 8, 9; 0, 0, 0, 10];
b = [10, 18, 17, 10]';
x = zeros(size(A, 2), 1);
for i = size(A, 1):-1:1,
x(i) =b(i)/A(i, i);
j = 1:(i-1);
b(j) = b(j) - A(j, i)*x(i);
end
x
二、 C++实现
2.1 数据结构的定义:
typedef vector<int> Vec;
typedef vector<Vec > Mat;
这里需要注意Mat
矩阵的初始化不同于Vec
向量的初始化,Vec
可以直接动态的增加元素,Mat
却需要首先添加Vec
作为其每一行,然后再在Vec
中增加元素。
对向量Vec
和矩阵Mat
的初始化分别为:
int arr[] = { 10, 18, 17, 10 };
Vec b(arr, arr + 4);
Mat A;
for (int i = 0; i < 4; ++i)
A.push_back(Vec(4));
int count = 1;
for (int i = 0; i < 4; ++i)
for (int j = i; j < 4; ++j)
A[i][j] = count++;
2.2 函数接口设计
数学形式:
为了模拟数学的感觉,变量名我们沿用数学的记号,在不引起混淆的前提下;
void solve(const Mat& A, const Vec& b, Vec& x);
2.3 implementation
void solve(const Mat& A, const Vec& b, Vec& x)
{
size_t row = A.size();
Vec b2 = b; // 为了保持输入的原始数据不被修改
for (int i = row - 1; i >= 0; --i)
{
x[i] = b2[i] / A[i][i];
for (int j = 0; j < i; ++j)
b2[j] -= A[j][i] * x[i];
}
}
这里需要注意两个问题:
unsigned类型变量为负时:
for (size_t i = row-1; i >= 0; --i)
// 或者unsigned i;
// typedef _W64 unsigned int size_t;
当i
减少到0,再执行减一的动作时,不会变成-1,而是4294967295
,也即(
232−1
)。此时循环不会退出,但如果试图索引x[4294967295]
,则会报错。
返回局部类对象的引用:
切勿返回局部类对象的引用(如下代码所示),因为函数调用完毕,编译器执行清栈操作,会调用局部对象的析构函数。对象如果不存在的话,何谈其引用呢。
Vec& solve(const Mat& A, const Vec& b)
{
size_t row = A.size();
Vec x(row);
...
return x;
}
Vec x = solve(A, b); // x.size() == 0
对问题2的改造:
如2.3节所示:
void solve(const Mat& A, const Vec& b, Vec& x);
// 在客户端完成对x的初始化
或者:
Vec solve(const Mat& A, const Vec& b);
// 返回的是类对象,而不是类对象的引用,这样在函数退出时,会调用Vec类的拷贝构造函数,创建并赋值一个临时(未名)对象
三、 实现方法之二
我们再来观察当 A <script type="math/tex" id="MathJax-Element-119">A</script>是四阶矩阵的计算流程:
- x[4] = b[4]/A[4][4]
- x[3] = (b[3] - A[3][4]*x[4])/A[3][3]
- x[2] = (b[2] - A[2][4]*x[4] - A[2][3]*x[3])/A[2][2]
- x[1] = (b[1] - A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]) / A[1][1]
A[1][4]*x[4] - A[1][3]*x[3] - A[1][2]*x[2]
这种A
的行和x
的乘积再加和的形式对应着线性代数的一个经典的概念就叫内积。我们便由此可以利用A
的行和x
向量已解部分的内积来进行相关的计算。
3.1 matlab
clear all; close all; clc
A = [1, 2, 3, 4; 0, 5, 6, 7; 0, 0, 8, 9; 0, 0, 0, 10];
b = [10, 18, 17, 10]';
x = zeros(size(A, 2), 1);
r = size(A, 1);
for i = r:-1:1,
j = i+1:r;
x(i) = (b(i) - A(i, j)*x(j)) / A(i, i);
end
x
因为未解部分已被初始化为0,可进一步将程序简化为:
for i = r:-1:1,
x(i) = (b(i) - A(i, :)*x) / A(i, i);
end
3.2 C++
void solve2(const Mat& A, const Vec& b, Vec& x)
{
size_t r = A.size();
for (int i = r - 1; i >= 0; --i)
{
int tmp = b[i];
for (int j = 0; j < r; ++j)
tmp -= A[i][j] * x[j];
x[i] = tmp / A[i][i];
}
}