一、消元过程
设线性方程组为:
式(1.1)
消元目标为一个上三角形方程组:
式(1.2)
第1步,设,记,用-li1乘以方程组式(1.1)的第1个方程,加到第i个方程(i=2,3,....,n),消去方程组式(1.1)的第2个直到第n个方程中含有x1的项,得到一个与式(1.1)等价的方程组式(1.3):
式(1.3)
其中有
第2步,设,记,用-li2乘以方程组式(1.3)的第2个方程,加到第i个方程(i=3,4,....,n),消去方程组式(1.,3)的第n-2个方程中含有x2的项,得到一个与式(1.3)等价的方程组式(1.4):
式(1.4)
其中有
类似的,若(i=3,4,...,n-1),完成第n-1步后,则可得一个等价的上三角形方程组
式(1.5)
由方程式(1.1)转化为方程组式(1.5)是Gauss消去法的消元过程,该过程的计算公式为:
k = 1,2,3,...,n-1,计算:
式(1.6)
二、回代过程
通过消元过程,将原来的方程组式(1.1)转化为上三角形方程组式(1.5),当时,该上三角形方程组可以通过如下公式求解:
式(1.7)
即回代过程
三、代码实现
根据式(1.6)和(1.7)的迭代过程,该算法实现的代码如下:
%% Gauss消元法的Matlab实现
function [r_matrix] = GaussElimination(Coefficient_matrix,Load_matrix)
% 2017-11-05 xh_scu
% inputs:
% Coefficient_matrix:系数矩阵,为n*n维方阵
% Load_matrix :载荷矩阵,为n*1维矩阵
% outputs:
% r_matrix:计算结果向量,为n*1为矩阵
% 判断输入矩阵维度是否满足要求
[row_coeff,col_coeff] = size(Coefficient_matrix);
[row_load,~] = size(Load_matrix');
% 初始化r_matrix矩阵
r_matrix = zeros(row_load,1);
% 判断输入的维度是否满足要求
if (row_coeff ~= col_coeff) || (row_coeff ~= row_load)
% 不满足则输出错误提示
print('输入错误!');
else
% 满足则进行下一步运算
for k = 1:row_coeff-1
% 检查主对角元素第i行的第i个元素是否为0
if Coefficient_matrix(k,k) == 0
print('主对角元素错误!');
else
% 循环计算第k+1行到最后一行
for i = k+1:row_coeff
L_ik = Coefficient_matrix(i,k) / Coefficient_matrix(k,k); %更新L_ik
% 更新每一行从第i个元素开始后的所有元素
for j = k+1:col_coeff
Coefficient_matrix(i,j) = Coefficient_matrix(i,j) - ...
L_ik*Coefficient_matrix(k,j); % 更新a(i,j)
end
Load_matrix(i) = Load_matrix(i) - Load_matrix(k)*L_ik; % 更新b(i)
Coefficient_matrix(i,1) = 0;
end %for循环结束
end % 条件2结束
end % for循环结束
% 回代过程
r_matrix(row_coeff) = Load_matrix(row_coeff)/Coefficient_matrix(row_coeff,col_coeff);
for k = row_coeff-1:-1:1
sum_temp = 0;
for j = k+1:col_coeff
sum_temp = sum_temp + Coefficient_matrix(k,j)*r_matrix(j);
end
r_matrix(k) = (Load_matrix(k) - sum_temp)/Coefficient_matrix(k,k);
end
end % 条件判断结束
end % 函数结束
四、测试分析
设输入a = [1,2;5,6], b = [12,24],则计算结果为result=[-6,9]。通过手算验证该结果为正确结果。
下面对算法的时间代价进行测试,输入维度分别取100,200,300,400,500,600,700,800,900,1000,输入元素采用1-100的整数,计算时间代价如下表所示。
维度 | 100 | 200 | 300 | 400 | 500 | 600 | 700 | 800 | 900 | 1000 |
时间 | 0.0280 | 0.2230 | 0.7520 | 1.7890 | 3.4770 | 6.1610 | 10.2220 | 15.0800 | 21.7500 | 31.4570 |
计算时间示意图