(完整)高斯消去法解线性方程的Matlab程序.doc

1151091  杨晨辉 高斯消去法解线性方程的  Matlab 程序 方法一： function x = gauss(A,b) n = length(b); for k = 1 : n-1 if A(k,k)==0 fprintf( 'Error: the %dth pivot element equal to zero!\n'  ,k); return ; end index = [k+1:n]; m = -A(index,k)/A(k,k); A(index,index) = A(index,index) + m*A(k,index); b(index) = b(index) + m*b(k); end x = zeros(n,1); x(n) = b(n)/A(n,n); for i = n-1:-1:1 x(i) = ( b(i) - A(i,[i+1:n])*x([i+1:n]) )/A(i,i); end 运行结果： >> A=[1 1.355 1.4 2; 3 3.5 0.22 1; 0.5 2 2.1 3; 0.3 0.1 -0.55 2]; b=[2.00,1.00,0.55,3.00]' b = 2.0000 1.0000 0.5500 3.0000 gauss(A,b) ans = 2.5225 -2.2313 0.0177 1.2381 方法二： 矩阵求逆： function [ B ] = qiuni( A ) %UNTITLED Summary of this function goes here Detailed explanation goes here n=numel(A); r=rank(A); B=eye(r); if n==r^2 for k=1:r for i=1:r for j=1:r ii=i-1; jj=j-1; if ii==0 && jj==0; ii=r; jj=r; B(ii,jj)=1/A(1,1); elseif ii==0 && jj~=0; ii=r; B(ii,jj)=-A(1,j)/A(1,1); elseif jj==0 && ii~=0; jj=r; B(ii,jj)=A(i,1)/A(1,1); else B(ii,jj)=A(i,j)-A(i,1)*A(1,j)/A(1,1); end end end A=B; B=eye(r); end B=A; else msgbox( ' 矩阵不可逆 ' , 'message' , 'warn' ); end end 当方程不可逆时： 第二种： function [ B ] = lyxqinv( A ) %UNTITLED Summary of this function goes here % Detailed explanation goes here n=numel(A); r=rank(A); B=eye(r); if n==r^2 E=eye(r); A=[A E]; for k=1:(r-1) for i=(k+1):r for j=(k+1):2*r A(i,j) = A(i,j)-A(i,k)* A(k,j) / A(k,k); end j=k; A(i,j) = A(i,j)-A(i,k)* A(k,j) / A(k,k); end end for k=2:r for i=1:(k-1) for j=(k+1):2*r A(i,j) = A(i,j)-A(i,k)*A(k,j) / A(k,k); end j=k; A(i,j) = A(i,j)-A(i,k)* A(k,j) / A(k,k); end end for i=1:r for j=1:2*r B(i,j)=A(i,j)/A(i,i); end end B=B(:,(r+1):2*r); else msgbox( ' 矩阵不可逆 ' , 'message' , 'warn' ); end end 方程求解： function X=jie(A,b); [m,n]=size(A);B=[A,b]; RA=rank(A);RB=rank(B);format rat if RA==RB & RA==n %判断方程有唯一解 X=A\b; else if RA==RB & RA<n %判断方程有无穷解 X=A\b %求特解 C=null(A, 'R' ); %求 AX=0的基础解系 else X=' 方程无解 ' ; %判断方程无解 end end 上图的方程有两个，第一个有解，显示结果；第二个无解，显示 方程无解 当方程有无穷解时，显示其特解