MATLAB入门学习-#6-Jacobi、Gauss-Seidel、SOR迭代法编程练习
这三种迭代法是在数值分析课程里学到的,都是求解线性方程组用的,相关的知识可以在数值计算方法的第三章里面找,可以说已经把最核心的怎么算和为什么这么算给交代了,剩下的任务就仅仅是编程实现就完了。
毕竟自己也没怎么编过程序,更是很久以来第一次写function,所以还是比较简陋的,但是也能完成相应的操作。在此就仅仅是记录一下自己学习matlab的历程把…
1.Jacobi迭代法
function [B,g,x]=jacobi(A,b,times,x0)
%[B,g,x]=jacobi(A,b,times,x0)
%
%B为jacobi迭代矩阵
%g为jacobi变形后的b
%x为方程组的解
%A为系数矩阵
%b为矩阵右侧
%times为迭代次数
%x0为初始x
D=diag(diag(A));
c=size(A);
I=eye(c);
B=I-D\A; %使用左除或者右除要比使用inv计算快
g=D\b;
if vrho(B)>=1 %求B的谱
error('jacobi迭代法不收敛!');
end
i=1;
while i<=times
x=B*x0+g; %来自书上
x0=x;
disp(['第',num2str(i),'次']);
disp(x);
i=i+1;
end
2.Gauss-Seidel迭代法
function [M,x]=gauss_seidel(A,b,x0,times)
%[M,x]=gauss_seidel(A,b,x0,times)
%
%M为jacobi迭代矩阵
%x为方程组的解
%A为系数矩阵
%b为矩阵右侧
%times为迭代次数
%x0为初始x0
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
M=(D-L)\U;
if vrho(M)>=1 %求B的谱
error('jacobi迭代法不收敛!');
end
i=1;
while i<=times
x=((D-L)\U)*x0+(D-L)\b; %来自书上
x0=x;
disp(['第',num2str(i),'次']);
disp(x);
i=i+1;
end
3.SOR迭代法(松弛法)
function [M,x]=SOR(A,b,x0,times,w)
%[M,x]=SOR(A,b,x0,times,w)
%
%M为jacobi迭代矩阵
%x为方程组的解
%A为系数矩阵
%b为矩阵右侧
%times为迭代次数
%x0为初始x0
%w为参数因子
D=diag(diag(A));
L=D-tril(A);
U=D-triu(A);
M=(D-w.*L)\((1-w).*D+w.*U);
if vrho(M)>=1 %求B的谱
error('jacobi迭代法不收敛!');
end
i=1;
while i<=times
x=M*x0+(D-w.*L)\(w.*b); %来自书上
x0=x;
disp(['第',num2str(i),'次']);
disp(x);
i=i+1;
end