前言

解基础线性方程组并不困难，基础方法在小学时就以鸡兔同笼问题呈现过了。但当线性方程组的量级扩大到上万倍时，就需要计算机参与运算了。本文介绍常见的三种方法：高斯消去法，列主元消去法，直接三角分解法

方法

高斯消去法

求解公式

此方法与小学时学习的消元法无异，核心思想就是消元

$$ \left \{ \begin{array}{c} a_{11}x_1+a_{12}x_2\cdots+a_{1n}x_n=b_1 \\ a_{21}x_1+a_{22}x_2\cdots+a_{2n}x_n=b_2 \\ \vdots \\ a_{n1}x_1+a_{n2}x_2\cdots+a_{nn}x_n=b_n \end{array} \right. $$

可以写为矩阵形式

$$ \begin{pmatrix} a_{11} & \cdots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \cdots & a_{nn} \end{pmatrix} \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} =\begin{pmatrix} b_{1} \\ b_{2} \\ \vdots \\ b_{n} \end{pmatrix} $$

简记为$Ax=b$

经过n次消元可得

$$ \begin{pmatrix} a_{11}^{(1)} & \cdots & a_{1n}^{(1)} \\ & \ddots & \vdots \\ & & a_{nn}^{(n)} \end{pmatrix} \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} =\begin{pmatrix} b_{1}^{(1)} \\ b_{2}^{(2)} \\ \vdots \\ b_{n}^{(n)} \end{pmatrix} $$

得求解公式

$$ \left \{ \begin{array}{c} x_n=b_n^{(n)}/a_{nn}^{(n)},\\ x_k=(b_k^{(k)}-\sum_{j=k+1}^na_{kj}^{(k)}x_j)/a_{kk}^{(k)}, k=n-1,n-2\cdots,1. \end{array} \right. $$

代码

main.m

clear;clc
A=input('(A,b)=');
n=size(A,1);
for i=1:n
    for j=i+1:n
        A(j,:)=-A(j,i)/A(i,i)*A(i,:)+A(j,:);
    end
    disp(['(A(' num2str(i) '),b)=']);
    disp(A);
end

x(n,1)=0;
x(n)=A(n,n+1)/A(n,n);
for i=size(A,1)-1:-1:1
    x(i)=((A(i,n+1))-f_sum(A,x,i,n))/A(i,i);
end

disp(x);

f_sum.m

function sum=f_sum(A,x,i,n)
sum=0;
for j=i+1:n
    sum=sum+A(i,j)*x(j);
end

缺点

• 该算法时间复杂度较高，不适合规模大的计算
• 在主元为绝对值较小的浮点数时结果误差较大，不适合精度要求高的运算
• 出现主元为0的情况无法计算

列主元消去法

求解公式

如上文所提到的，高斯消去法在主元为较小的浮点数时会出现较大误差
此时可以引入新方法，即列主元消去法
为避免消元过程中除去小主元所带来的误差，在每次消元前将本列绝对值最大的主元行与本行对换

代码

main.m

clear;clc
A=input('(A,b)=');
n=size(A,1);

for i=1:n
    [p,]=find(A(:,i)==max(A(i:n,i)),1);
    A([i p],:)=A([p i],:);
    for j=i+1:n
        A(j,:)=-A(j,i)/A(i,i)*A(i,:)+A(j,:);
    end
    disp(['(A(' num2str(i) '),b)=']);
    disp(A);
end

x(n,1)=0;
x(n)=A(n,n+1)/A(n,n);
for i=size(A,1)-1:-1:1
    x(i)=((A(i,n+1))-f_sum(A,x,i,n))/A(i,i);
end

disp(x);

f_sum.m

function sum=f_sum(A,x,i,n)
sum=0;
for j=i+1:n
    sum=sum+A(i,j)*x(j);
end

解析

其中，寻找绝对值最大主元代码如下：

 [p,]=find(A(:,i)==max(A(i:n,i)),1);

对于find()函数，常见用法如下：

i = find(X)
i = find(X, k) % 返回前k个不为零的元素索引值
i = find(X, k, 'first') % 返回前k个不为零的元素索引值
i = find(X, k, 'last') % 返回后k个不为零的元素索引值
[r,c] = find(X, ...) % 返回矩阵X中非零元素的行和列的索引值
[r,c,v] = find(X, ...) % 返回矩阵X中非零元素的行和列的索引值以及该不为零元素

这里需要获取最大值所在第几行的值，所以只需要取[r,c]的第一个值，第二个留空即可
如果有多个条件满足find，该函数会返回一个向量，这里只需要第一个最大值，所以在find中加了参数，即find(A,1)

这里行对换的代码如下：

A([i p],:)=A([p i],:);

列对换同理，不再赘述

直接三角分解法

求解公式

对于系数矩阵$A$，在其为非奇异矩阵的情况下，有分解式
$$A=LU$$
其中$L$为单位下下三角矩阵，$U$为上三角矩阵

$$ A=\begin{pmatrix} 1 & & & \\ l_{21} &1 & &\\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \end{pmatrix} \begin{pmatrix} u_{11} &u_{12} &\cdots & u_{1n}\\ &u_{22} &\cdots &u_{2n}\\ & & \ddots & \vdots \\ && & u_{nn} \end{pmatrix} $$

可以得到计算公式
(1)

$$ u_{1i}=a_{1i},(i=1,2\cdots ,n)， \\ l_{i1}=a_{i1}/u_{11},(i=2,3\cdots ,n). $$

计算$U$的第$r$行，$L$的第$r$列元素$(r=2,3\cdots ,n)$
(2)

$$ u_{ri}=a_{ri}-\sum_{k=1}^{r-1}l_{rk}u_k,i=r,r+1,\cdots ,n; \\ $$

(3)

$$ l_{ir}=(a_{ir}-\sum_{k=1}^{i-1}l_{ik}u_{kr}/u_{rr},i=r+1,\cdots ,n,且r\neq n. $$

求解$Ly=b$，$Ux=y$：
(4)

$$ \left \{ \begin{array}{c} y_1=b_1, \\ y_i=b_i-\sum_{k=1}^{i-1}l_{ik}y_k,i=2,3,\cdots,n; \end{array} \right. $$

(5)

$$ \left \{ \begin{array}{c} x_n=y_n/u_{nn}, \\ x_i=(y_i-\sum_{k=i+1}^ni_{ik}x_k)/i_{ii},i=n-1,n-2,\cdots ,1; \end{array} \right. $$

代码

main.m

clear;clc

%% input
A=input('A=');
b=input('b=');
n=size(A,1);

%%l
for i=1:n
    l(i,i)=1;
end

%% u1i lil
u(1,:)=A(1,:);
l(2:n,1)=A(2:n,1)./u(1,1);

%% uri lir
for r=2:n
    for i=r:n
        u(r,i)=A(r,i)-f_sum1(l,u,r,i);
        if i==r
            continue
        end
        l(i,r)=(A(i,r)-f_sum2(l,u,r,i))/u(r,r); 
    end
end

%% y
y(1)=b(1);
for i=2:n
    y(i)=b(i)-f_sum3(l,y,i);
end

%% x
x(n)=y(n)/u(n,n);
for i=n-1:-1:1
    x(i)=(y(i)-f_sum4(u,x,i,n))/u(i,i);
end

%% output
x

f_sum1.m

function sum=f_sum1(l,u,r,i)
sum=0;
for k=1:r-1
    sum=sum+l(r,k)*u(k,i);
end

f_sum2.m

function sum=f_sum2(l,u,r,i)
sum=0;
for k=1:r-1
    sum=sum+l(i,k)*u(k,r);
end

f_sum3.m

function sum=f_sum3(l,y,i)
sum=0;
for k=1:i-1
    sum=sum+l(i,k)*y(k);
end

f_sum4.m

function sum=f_sum4(u,x,i,n)
sum=0;
for k=i+1:n
    sum=sum+u(i,k)*x(k);
end

仓库地址

本文所有涉及的代码均已上传至GitHub
仓库地址