用matlab解线性方程组.docVIP

用matlab解线性方程组 2008-04-12 17:00 一。高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[] [n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1 a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去 end x=[0 0 0 0] %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end 2.列主高斯消去法 *由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。 直接编写命令文件 a=[] d=[] [n,n]=size(a); c=n+1 a(:,c)=d; %(增广) for k=1:n-1 [r,m]=max(abs(a(k:n,k))); %选主 m=m+k-1; %(修正操作行的值) if(a(m,k)~=0) if(m~=k) a([k m],:)=a([m k],:); %换行 end a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去 end end x=[0 0 0 0] %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end 3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果 a=[0 1 2 3;9 11 23 34;62.5 23.4 15.5 17.2;192.01 124 25.1 59.3] d=[1;1;1;1] 顺序高斯消去法：提示“Warning: Divide by zero.” x =NaN NaN NaN NaN 列主高斯消去法：x =-1.2460 2.8796 5.5206 -4.3069 由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。 4. 将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“15.5”改为“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。 x =-0.8081 1.8226 3.5568 -2.7047 很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。 这是因为系数矩阵的条件数很大：cond(a)=2.0695e+003 二。迭代法 J迭代公式 a=[5 2 1;-1 4 2;2 -3 10] d=[-12;20;3] x=[0;0;0]; %初始向量 stop=1.0e-4 %迭代的精度 L=-tril(a,-1) U=-triu(a,1) D=inv(diag(diag(a))) X=D*(L+U)*x+D*d; % J迭代公式 n=1; while norm(X-x,inf)=stop % 时迭代中止否则继续 x=X; X=D*(L+U)*x+D*d; n=n+1; end X n G-S迭代公式 a=[5 2 1;-1 4 2;2 -3 10] x=[0;0;0]; %初始向量 d=[-12;20;3] stop=1.0e-4 L=-tril(a,-1) U=-triu(a,1) D=(diag(diag(a))) X=inv(D-L)*U*x+inv(D-L)*d; % G-S迭代公式 n=1; while norm(X-