matlab elimination tree,Gauss Elimination算法分析与实现

高斯消去法分为两个过程：第一步是前向消元(forward elimination)，也就是将系数矩阵转化成上三角矩阵的过程；第二步是回代(back substitution)过程，自底向上求解方程组的过程。

选择主元(pivoting)，主元(pivot)一定不能为0.且主元应当尽可能大。所以一般情况下，我们采取部分选主元(partial pivoting)。选择主元为所在各列中绝对值最大的元素(或者非0亦可)。然后将对应的该行与主元行进行交换，然后消元。

function x = gausselimination(A,b)

ab = [A,b];

[r,c] = size(ab);

for k=1:r-1

[rm,im]=max(abs(ab(k:r,k)));//im是最大元素对应的索引

im = im +k-1;

if(ab(im,k)~=0)

if(im~=k)

ab([k im],:)=ab([im k],:);

end

end

for i=k+1:r

ab(i,k:c)=ab(i,k:c)-ab(i,k)/ab(k,k)*ab(k,k:c);

end

end

x=zeros(r,1);

x(r)=ab(r,c)/ab(r,r);

for i=r-1:-1:1

x(i)=(ab(i,c)-ab(i,i+1:r)*x(i+1:r))/ab(i,i);

end

以上是matlab代码实现(保存为gausselimination.m)。以下是运行结果：

>> A = [10 -7 0;-3 2 6;5 -1 5]

A =

10    -7     0

-3     2     6

5    -1     5

>> b = [7;4;6]

b =

7

4

6

>> gausselimination(A,b)

ans =

0.0000

-1.0000

1.0000

In linear algebra, Gaussian elimination (also known as row reduction) is an algorithm for solving systems of linear equations. It is usually understood as a sequence of operations performed on the associated matrix of coefficients. This method can also be used to find the rank of a matrix, to calculate the determinant of a matrix, and to calculate the inverse of an invertible square matrix. The method is named after Carl Friedrich Gauss (1777–1855), although it was known to Chinese mathematicians as early as 179 CE (see History section).

The process of row reduction makes use of elementary row operations, and can be divided into two parts. The first part (sometimes called Forward Elimination) reduces a given system to row echelon form, from which one can tell whether there are no solutions, a unique solution, or infinitely many solutions. The second part (sometimes called back substitution) continues to use row operations until the solution is found; in other words, it puts the matrix into reduced row echelon form.

算法解释如下：

接下来一起看一个小程序，它的目的不过是为了求解线性方程组Ax+b=x,化简之后，可以得到(I -A)x=b。请注意在MATLAB中的语法是怎样的。为此，首先我们初始化矩阵A和列向量b

>> A=[0.85 0.04;-0.04 0.85]

A =

0.8500    0.0400

-0.0400    0.8500

>> b=[0;1.6]

b =

0

1.6000

注意观察接下来的操作

>> I = eye(2,2)

I =

1     0

0     1

>> x=(I-A)\b

x =

2.6556

9.9585

这样就解得了这个方程的解。最关键的是这个符号“\”，那么我们不禁好奇，这个符号的功能是怎样实现的。

\ --“It’s amazing how much numerical linear algebra is contained in that single character”

为此介绍一个函数bslashtx

function x = bslashtx(A,b)

% BSLASHTX  Solve linear system (backslash)

% x = bslashtx(A,b) solves A*x = b

[n,n] = size(A);

if isequal(triu(A,1),zeros(n,n))

% Lower triangular

x = forward(A,b);

return

elseif isequal(tril(A,-1),zeros(n,n))

% Upper triangular

x = backsubs(A,b);

return

elseif isequal(A,A')

[R,fail] = chol(A);

if ~fail

% Positive definite

y = forward(R',b);

x = backsubs(R,y);

return

end

end

% Triangular factorization

[L,U,p] = lutx(A);

% Permutation and forward elimination

y = forward(L,b(p));

% Back substitution

x = backsubs(U,y);

% ------------------------------

function x = forward(L,x)

% FORWARD. Forward elimination.

% For lower triangular L, x = forward(L,b) solves L*x = b.

[n,n] = size(L);

x(1) = x(1)/L(1,1);

for k = 2:n

j = 1:k-1;

x(k) = (x(k) - L(k,j)*x(j))/L(k,k);

end

% ------------------------------

function x = backsubs(U,x)

% BACKSUBS.  Back substitution.

% For upper triangular U, x = backsubs(U,b) solves U*x = b.

[n,n] = size(U);

x(n) = x(n)/U(n,n);

for k = n-1:-1:1

j = k+1:n;

x(k) = (x(k) - U(k,j)*x(j))/U(k,k);

end