线性方程组迭代方法
用Jacobi迭代计算如下方程组
推导过程
编程思路
在求解线性方程组时, 程序中输入的是线性方程组的增广矩阵,并设定送
代结果的精度要求, 在开始迭代之前则可以进行以下判断:
- 检测该用户输入矩阵是否有误,并判断该线性方程组是否有唯一解;
- 判断系数矩阵的对角元是否有零元素。
由于迭代过程的本质是矩阵乘法,因此程序中直接调用 MATLAB 的矩阵乘
法来完成送代,选代均以零向量为初始向量,并且在送代之前判断该送代过程是
否收敛。 倘若收敛,则输出结果并打印迭代次数。
MATLAB代码
function函数:
function solution = Jacobi(Ab,epsilon)
%% Ab为输入的增广矩阵
%% epsilon为输入的精度要求
%% 输入参数检查
if nargin ==1
disp('请输入精度要求epsilon')
return
end
row = size(Ab,1);
col = size(Ab,2);
if ndims(Ab)~= 2 | col - row ~= 1
disp('矩阵的大小有误,不能采用Jacobi迭代法')
return
end
A = Ab(:,1:row);
b = Ab(:,col);
ddet = abs(det(A));
ddiag = abs(det(diag(diag(A))));
if ddet < eps | ddiag < eps
disp('该方程的系数矩阵行列式不为0,方程组无解或有无穷解,或系数矩阵的对角线有零元,不能采用Jacobi迭代法')
return
end
%%提取上下三角矩阵及对角矩阵
U = -triu(A,1); %提取对角元素为0的上三角矩阵
L = -tril(A,-1); %提取对角元素为0的下三角矩阵
Dinv = diag(diag(A).^(-1));
B = Dinv * (L + U);
f = Dinv * b;
if max(abs(eig(B))) > 1
disp('迭代法不收敛!!!')
return
end
%% 迭代过程
error = 10;
n = 0;
start = zeros(row,1);
xk = start;
xknext = B * xk +f;
error = norm(xk - xknext);
while error > epsilon
xk = xknext;
xknext = B * xk + f;
error = norm(xk - xknext);
n = n+1;
end
fprintf('Jacobi迭代次数:%d',n)
solution = xknext;
主函数:
Ab = [10 -2 -1 3;-2 10 -1 15;-1 -2 5 10]
Jacobi(Ab,1e-8)
**注意:**function函数和主函数是在同一个文件夹路径下的两个文件!