大家对高斯消去法应该比较了解了,高代和线性代数中做的已经不少了,但是计算机实现的时候还是要注意一些东西,
%列选主元的高斯消去法
function X=lufact_my(A,B)
%Inpiut A 是系数矩阵,B是右端项
%Output x是解
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
C=zeros(1,N);
R=1:N;
k=1;
while k<=N-1
%求列中最大的值赋给max
[max1,j]=max(abs(A(k:N,k)));
%交换行
C=A(k,:);%C为A的第k列的值
A(k,:)=A(j+k-1,:); %将A的第K列赋为最大
A(j+k-1,:)=C;
d=R(k);
R(k)=R(j+k-1);
R(j+k-1)=d;
%主元为0的情况
if A(k,k)==0
'A is singular. no unique solution'
break
end
%化为上三角
for m=k+1:N
mult=A(m,k)/A(k,k);
A(m,k)=mult;
A(m,k+1:N)=A(m,k+1:N)-mult*A(k,k+1:N);
m=m+1;
end
k=k+1;
end
%对右端项做处理,但要保证行的交换相同,要注意R(k)的作用
Y(1)=B(R(1));
for k=2:N
Y(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);
end
X(N)=Y(N)/A(N,N);
for k=N-1:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end
演示---
>> A=[2,-1,3;4,2,5;2,1,2];
>> B=[1,4,5];
>> lufact_my(A,B)
ans =
9
-1
-6
>>
此算法不能求主元为零的,局限比较大。下一篇写一个改进的。。。