jacobi.m

来自「数值分析常用的一些公式」· M 代码 · 共 20 行

 function s=jacobi(a,b,x0)
 %a 为系数矩阵，b为右端向量，x0为初始值
 esp=1.0e-4;
 D=diag(diag(a));
 D=inv(D);
 L=tril(a,-1);%求出严格下三角矩阵
 U=triu(a,1); %求出严格上三角矩阵
 B=-D*(L+U);
 f=D*b;
 s=B*x0+f;
 while norm(s-x0)>=esp
     x0=s;
     s=B*x0+f;
 end
 return
 
 %a=[10 -2 -1;-2 10 -1;-1 -2 5];
 %b=[3 15 10]';
 %x0=[0 0 0]';
 %s=jacobi(a,b,x0)