这学期开了一门数值计算方法的课,要求MATLAB编程。
因为第一次学着用MATLAB,不少矩阵操作都不知道,代码写起来磕磕绊绊的,遂简单记录一下,注释应该能够让没用过MATLAB的同学也能明白每一行的作用。
过程完全按照课本所讲,教科书式代码
%Jacobi迭代法 形如 x=Mx+g
%系数矩阵A 对角矩阵D 常数项b 迭代矩阵M 常数矩阵g 初始条件x0
A=input('');
b=input('');
%保留初始x0为x00
x0=input('');
%很不错的diag用法
D=diag(diag(A));
%取合适的单位矩阵
[m,n]=size(A);
B=eye(m,n)-inv(D)*A;
g=inv(D)*b;
for i=1:10000
x1=B*x0+g;
%inf意味着所有元素中最大的那个,根据迭代终止条件定义应如此
n=norm(x1-x0,inf);
x0=x1;
%不妨取误差为0.0001
if(n<0.0001)
break
end
end
i
x1
n
%Guass-Seidel迭代法 形如 x=(D-L)'*Ux+(D-L)'b
%系数矩阵A 对角矩阵D 常数项b 下三角矩阵L 上三角矩阵U 迭代矩阵M 常数矩阵g 初始条件x0
A=input('');
b=input('');
%保留初始x0为x00
x0=input('');
%很不错的diag用法
D=diag(diag(A));
%tril triu分别为矩阵A的取上三角下三角函数,后面参数-1 ,1意思是从第-1/1个对角线开始取上下三角
L=tril(A,-1);
U=triu(A,1);
g=inv(D-L)*b;
for i=1:10000
x1=inv(D-L)*U*x0+g;
%inf意味着所有元素中最大的那个,根据迭代终止条件定义应如此
n=norm(x1-x0,inf);
x0=x1;
%不妨取误差为0.0001
if(n<0.0001)
break
end
end
i
x1
n