1、函数输入形参要求如下,也可以修改得更好、更适合你自己.
function Xsolve=myJacobi(A.b.X0 mytol)
%A为系数矩阵
點为右端常向量
X0为初始层,默认取0.
%mytol表示允许误差.要求默认值为le-b.
%Xsolve表示收敛情况下满足精度的数值解
Xnargin表示输入参数个数,下面这一段是简单的判断,你可以编写的更好if nargin<4 %如果参数<4.
mytol=le-b:
if nargin<3 %如果参数3
X0=zeros(size(b));
if nargins2
error(‘输入参数不足‘);
end
end
end
2.第一步描述出矩阵DL.U
3、第二步要判断系数矩阵是否可逆,用“abs(det(A))
5、第四步用迭代法基本定理(计算谱半径)判断迭代法是否收敛,要在命令窗口输出是否收敛的结论:
6、第五步在收敛的基础上进行迭代格式的程序编写:
7、迭代格式终止循环的判断方式用后一层和前一层计算结果的差的绝对值的最大值小于某个给定的小量mytol. 比如le-b,即.max(abs(A-B))emytol也可以用其它方式终止循环:
8、收敛的情况下,和精确解做简单比较,精确解可以采用Ab的形式:
9、注意:命令窗口中显示的结果不是真实结果,只是在输入格式的条件下进行四舍五入的显示,要想看到小数点后多个位的值,可以采用format long 等形式.
下面是一个用Jacobi迭代格式解线性方程组的MATLAB函数实现,实现了上述的各个步骤:
function Xsolve = myJacobi(A, b, X0, mytol)
% A为系数矩阵,b为右端常向量,X0为初始层,默认取0.
% mytol表示允许误差。要求默认值为le-6。
% Xsolve表示收敛情况下满足精度的数值解X
% 判断输入参数
if nargin < 4 % 如果参数<4
mytol = 1e-6;
if nargin < 3 % 如果参数<3
X0 = zeros(size(b));
if nargin < 2 % 如果参数<2
error('输入参数不足');
end
end
end
% 计算矩阵DL、U
D = diag(diag(A));
L = tril(A, -1);
U = triu(A, 1);
% 判断矩阵是否可逆
if abs(det(A)) < mytol
error('系数矩阵不可逆');
end
% 判断主对角元素是否为零
if any(diag(A) == 0)
error('系数矩阵主对角元素有零元素');
end
% 计算迭代矩阵T和常数向量C
T = -D \ (L + U);
C = D \ b;
% 计算谱半径
rhoT = max(abs(eig(T)));
% 判断迭代法是否收敛
if rhoT >= 1
error('迭代法不收敛');
else
fprintf('迭代法收敛,谱半径为 %f\n', rhoT);
end
% 迭代格式
Xn = X0;
diff = inf;
while diff > mytol
Xn1 = T * Xn + C;
diff = max(abs(Xn1 - Xn));
Xn = Xn1;
end
Xsolve = Xn;
% 与精确解比较
if nargin >= 2 % 如果有精确解,则比较
Xexact = A \ b;
err = norm(Xsolve - Xexact);
fprintf('与精确解比较,误差为 %e\n', err);
end
其中,利用MATLAB自带的运算符 tril 和 triu 分别计算下三角矩阵和上三角矩阵,利用 eig 计算矩阵的特征值, norm 计算向量的范数。注意,该函数实现了默认输入参数的处理,即如果调用时不传入某些参数,则使用默认值,这样方便用户的调用。最后,函数中进行了误差比较,可以对比结果的精度。
不知道你这个问题是否已经解决, 如果还没有解决的话: