主题:如何用列主元高斯消去法解线形方程?
yhfsmile
[专家分:0] 发布于 2007-10-28 19:29:00
如何用列主元高斯消去法解线形方程?
回复列表 (共4个回复)
沙发
mathsdzf [专家分:20] 发布于 2007-10-30 08:34:00
function x=GEPivShow(A,b,ptol)
%GEPivShow Show steps in Gauss elimination with partial pivoting and back
%substitution
%Snopsis: x=GEPivShow(A,b)
% x=GEPivShow(A,b,ptol)
%input A,b the coefficient matrix and right hand side vector
%ptol=(optional) tolerance for detection of zero pivot
%Defaut :ptol=50*eps
%output x=solution vector ,if solution exists
if nargin<3 ,ptol=50*eps;end
[m,n]=size(A);
if m~=n,error('A matrix needs to be square');end
nb=n+1;Ab=[A b]; %augmented system
fprintf('\n Begin elimination with augmented system:\n');disp(Ab);
% elimination
for i=1:n-1; %loop over pivot row
[pivot p]=max(abs(Ab(i:n,i)));%value and index of largest %%%%%pivot
ip=p+i-1; %p is index in subvector i:n
if ip~=i %ip is true row index of desired pivot
fprintf('\nSwap rows %d and %d new pivot=%f\n',i,ip,Ab(ip,i));
Ab([i ip],:)=Ab([ip i],:);
disp(Ab);% perform the swap
end
pivot=Ab(i,i);
if abs(pivot)<ptol,error ('zero pivot encountered after row exchange');end
for k=i+1:n %k=index of next row to be eliminated
Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);
end
fprintf('\n after elimination in column %d with pivot =%f \n',i,pivot);
disp(Ab);
end
%back substitution
x=zeros(n,1); %preallocate memory for and initialize x
x(n)=Ab(n,nb)/Ab(n,n);
for i=n-1:-1:1;
x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);
end
板凳
Guassfans [专家分:4090] 发布于 2007-10-31 12:25:00
function x=detGauss(a,b)
% Gauss列主元消去法
[n,m]=size(a);
nb=length(b);
det=1;%存储行列式值
x=zeros(n,1);
for k=1:n-1
amax=0;% 选主元
for i=k:n
if abs(a(i,k))>amax
amax=abs(a(i,k));r=i;
end
end
if amax<1e-10
return;
end
if r>k %交换两行
for j=k:n
z=a(k,j);a(k,j)=a(r,j);a(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;det=-det;
end
for i=k+1:n %进行消元
m=a(i,k)/a(k,k);
for j=k+1:n
a(i,j)=a(i,j)-m*a(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*a(k,k);
end
det=det*a(n,n);
for k=n:-1:1 %回代
for j=k+1:n
b(k)=b(k)-a(k,j)*x(j);
end
x(k)=b(k)/a(k,k);
end
3 楼
yhfsmile [专家分:0] 发布于 2007-10-31 18:15:00
谢谢!但是在运行的时候回出现问题,说没定义detGauss,怎么回事?
4 楼
Guassfans [专家分:4090] 发布于 2007-11-01 08:30:00
新建M文件,以detGauss为名,将内容复制于其中!
我来回复