jacobi_gs_sor.m

来自「数值分析, 同济大学教材<<现代数值数学和计算>>数值分」· M 代码 · 共 81 行

%gauss迭代
%author:skysmiler
%create date:2007-12-1
%调用格式：gauss(a,b,x0,eps）
%a为系数矩阵,b为方程组常数向量，x0为初值,eps为精度。
function res=gauss(a,b,x0,eps) 
D=diag(diag(a));%求矩阵对角阵
L=-tril(a,-1);%求严格下三角矩阵
U=-triu(a,1);%求严格上三角矩阵
C=inv(D-L);
B=C*U;
if max(abs(eig(B)))>=1
    fprintf('B的谱半径>=1.');
    max(abs(eig(B)))
    return; 
end 
g=C*b;
res=B*x0+g;
n=1;
while norm(res-x0)>=eps 
    x0=res;
    res=B*x0+g;
    n=n+1;
end 
n

%jacobi迭代
%author:skysmiler
%create date:2007-12-1
%调用格式：jacobi(a,b,x0,eps）
%a为系数矩阵,b为方程组常数向量，x0为初值,eps为精度。
function res=jacobi(a,b,x0,eps) 
D=diag(diag(a));%求矩阵对角阵
D=inv(D);
L=-tril(a,-1);%求严格下三角矩阵
U=-triu(a,1);%求严格上三角矩阵
B=D*(L+U); 
if max(abs(eig(B)))>=1
    fprintf('B的谱半径>=1.');
    max(abs(eig(B)))
    return; 
end 
g=D*b; 
res=B*x0+g;
n=1;
while norm(res-x0)>=eps 
    x0=res; 
    res=B*x0+g; 
    n=n+1;
end  
n

%sor迭代
%author:skysmiler
%create date:2007-12-1
%调用格式：sor(a,b,x0,eps）
%a为系数矩阵,b为方程组常数向量，x0为初值,w为松弛因子，eps为精度。
function res=sor(a,b,x0,w,eps) 
D=diag(diag(a));%求矩阵对角阵
L=-tril(a,-1);%求严格下三角矩阵
U=-triu(a,1);%求严格上三角矩阵
C=inv(D-w.*L);
B=C*((1-w).*D+w.*U);
if max(abs(eig(B)))>=1
    fprintf('B的谱半径>=1.');
    max(abs(eig(B)))
    return; 
end 
g=w.*C*b;
res=B*x0+g;
n=1;
while norm(res-x0)>=eps 
    x0=res;
    res=B*x0+g;
    n=n+1;
end 
n