1. 程式人生 > >數值分析 Gauss-Seidel迭代法求解線性方程組 MATLAB程式實現

數值分析 Gauss-Seidel迭代法求解線性方程組 MATLAB程式實現

Gauss-Seidel迭代法 參考數值分析第四版 顏慶津著 P39

執行輸入為:


執行結果為:

以下是函式內容(儲存為gauss.m檔案,在MATLAB中執行):

%function [G,d,x,N]=gauss(A,b) %Gauss-Seidel迭代法求解線性方程組。迭代公式x(k+1)=Gx(k)+d,(k=0,1,…). %返回迭代矩陣G,d和方程組的解x,以及迭代次數N %2015.11.6  密密編寫  (*^__^*) 嘻嘻…… function [G,d,x2,N]=gauss(A,b) n=size(A); if n(1)~=n(2)     error('矩陣A不是方陣'); end %初始化 N=1;%迭代次數 L=zeros(n);%分解A=D+L+U,D是對角陣,L是下三角陣,U是上三角陣 U=zeros(n); D=zeros(n); G=zeros(n);%G=-inv(D+L)*U d=ones(n,1);%d=inv(D+L)*b x=ones(n,1); for i=1:n%初始化L和U     for j=1:n         if i<j             L(i,j)=A(i,j);         end         if i>j             U(i,j)=A(i,j);         end     end end for i=1:n%初始化D     D(i,i)=A(i,i); end G=-inv(D+L)*U;%初始化G d=(D+L)\b;%初始化d if vrho(G)>=1     error('Gauss-Seidel迭代法不收斂!'); end %迭代開始 x1=x; x2=G*x+d; while norm(x2-x1,inf)>10^(-4)     x1=x2;     x2=G*x2+d;     N=N+1; end x=x2; end