matlab 列主元素高斯消去法,高斯消元法和列主元素法

高斯消元法

设线性方程组

A

x

=

b

Ax=b

Ax=b，其中

(

A

(

1

)

,

b

(

1

)

)

=

(

a

11

(

1

)

a

12

(

1

)

a

13

(

1

)

…

a

1

n

(

1

)

b

1

(

1

)

a

21

(

1

)

a

22

(

1

)

a

23

(

1

)

…

a

2

n

(

1

)

b

2

(

1

)

a

31

(

1

)

a

32

(

1

)

a

33

(

1

)

…

a

3

n

(

1

)

b

3

(

1

)

⋮

⋮

⋮

⋱

⋮

⋮

a

n

1

(

1

)

a

n

2

(

1

)

a

n

3

(

1

)

…

a

n

n

(

1

)

b

n

(

1

)

)

(A^{(1)},b^{(1)})= \left(\begin{array}{cccccc} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ a_{21}^{(1)} & a_{22}^{(1)} & a_{23}^{(1)} & \ldots & a_{2 n}^{(1)} & b_{2}^{(1)} \\ a_{31}^{(1)} & a_{32}^{(1)} & a_{33}^{(1)} & \ldots & a_{3 n}^{(1)} & b_{3}^{(1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n 1}^{(1)} & a_{n 2}^{(1)} & a_{n 3}^{(1)} & \ldots & a_{n n}^{(1)} & b_{n}^{(1)} \end{array}\right)

(A(1),b(1))=⎝⎜⎜⎜⎜⎜⎜⎛​a11(1)​a21(1)​a31(1)​⋮an1(1)​​a12(1)​a22(1)​a32(1)​⋮an2(1)​​a13(1)​a23(1)​a33(1)​⋮an3(1)​​………⋱…​a1n(1)​a2n(1)​a3n(1)​⋮ann(1)​​b1(1)​b2(1)​b3(1)​⋮bn(1)​​⎠⎟⎟⎟⎟⎟⎟⎞​

用高斯消元法时，步骤如下：

计算

−

l

i

1

=

−

a

i

1

(

1

)

a

11

(

1

)

,

(

i

=

2

,

3

,

.

.

.

,

n

)

-l_{i1}=-\frac{a^{(1)}_{i1}}{a^{(1)}_{11}},(i=2,3,...,n)

−li1​=−a11(1)​ai1(1)​​,(i=2,3,...,n)，然后将矩阵的第一行乘这个系数分别加到每一行，这样，第一列的对角元一下就全是0。

(

A

(

2

)

,

b

(

2

)

)

=

(

a

11

(

1

)

a

12

(

1

)

a

13

(

1

)

…

a

1

n

(

1

)

b

1

(

1

)

0

a

22

(

2

)

a

23

(

2

)

…

a

2

n

(

2

)

b

2

(

2

)

0

a

32

(

2

)

a

33

(

2

)

…

a

3

n

(

2

)

b

3

(

2

)

⋮

⋮

⋮

⋱

⋮

⋮

0

a

n

2

(

2

)

a

n

3

(

2

)

…

a

n

n

(

2

)

b

n

(

2

)

]

\left(A^{(2)}, b^{(2)}\right)=\left(\begin{array}{cccccc} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ 0 & a_{22}^{(2)} & a_{23}^{(2)} & \ldots & a_{2 n}^{(2)} & b_{2}^{(2)} \\ 0 & a_{32}^{(2)} & a_{33}^{(2)} & \ldots & a_{3 n}^{(2)} & b_{3}^{(2)} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & a_{n 2}^{(2)} & a_{n 3}^{(2)} & \ldots & a_{n n}^{(2)} & b_{n}^{(2)} \end{array}\right]

(A(2),b(2))=⎝⎜⎜⎜⎜⎜⎜⎛​a11(1)​00⋮0​a12(1)​a22(2)​a32(2)​⋮an2(2)​​a13(1)​a23(2)​a33(2)​⋮an3(2)​​………⋱…​a1n(1)​a2n(2)​a3n(2)​⋮ann(2)​​b1(1)​b2(2)​b3(2)​⋮bn(2)​​⎦⎥⎥⎥⎥⎥⎥⎤​

然后用除去第一行和第一列的矩阵重复上述步骤，n-1步之后，将会得到

(

A

(

n

)

,

b

(

n

)

)

=

(

a

11

(

1

)

a

12

(

1

)

a

13

(

1

)

…

a

1

n

(

1

)

b

1

(

1

)

a

22

(

2

)

a

23

(

2

)

…

a

2

n

(

2

)

b

2

(

2

)

a

33

(

3

)

…

a

3

n

(

3

)

b

3

(

3

)

⋱

⋮

⋮

a

n

n

(

n

)

b

n

(

n

)

)

\left(A^{(n)}, b^{(n)}\right)=\left(\begin{array}{llllll} a_{11}^{(1)} & a_{12}^{(1)} & a_{13}^{(1)} & \ldots & a_{1 n}^{(1)} & b_{1}^{(1)} \\ & a_{22}^{(2)} & a_{23}^{(2)} & \ldots & a_{2 n}^{(2)} & b_{2}^{(2)} \\ & & a_{33}^{(3)} & \ldots & a_{3 n}^{(3)} & b_{3}^{(3)} \\ & & & \ddots & \vdots & \vdots \\ & & & & a_{n n}^{(n)} & b_{n}^{(n)} \end{array}\right)

(A(n),b(n))=⎝⎜⎜⎜⎜⎜⎜⎛​a11(1)​​a12(1)​a22(2)​​a13(1)​a23(2)​a33(3)​​………⋱​a1n(1)​a2n(2)​a3n(3)​⋮ann(n)​​b1(1)​b2(2)​b3(3)​⋮bn(n)​​⎠⎟⎟⎟⎟⎟⎟⎞​

这个时候消元完成，然后进行求解，如下：

{

x

n

=

b

n

(

n

)

/

a

n

n

(

n

)

x

i

=

(

b

i

(

i

)

−

∑

j

=

i

+

1

n

a

i

j

(

i

)

x

j

)

/

a

i

i

(

i

)

,

i

=

n

−

1

,

n

−

2

,

⋯

,

1

\left\{\begin{array}{l} x_{n}=b_{n}^{(n)} / a_{n n}^{(n)} \\ x_{i}=\left(b_{i}^{(i)}-\sum_{j=i+1}^{n} a_{i j}^{(i)} x_{j}\right) / a_{i i}^{(i)}, \quad i=n-1, n-2, \cdots, 1 \end{array}\right.

{xn​=bn(n)​/ann(n)​xi​=(bi(i)​−∑j=i+1n​aij(i)​xj​)/aii(i)​,i=n−1,n−2,⋯,1​

matlab代码如下：

function x=gaosi(A,b)

[rows,cols]=size(A);

if rows~=cols

error("rows != cols")

end

n=rows;

Ab=[A,b];

%消元

for i=1:n-1

l= Ab((i+1):end,i)/Ab(i,i);

rep_Ab_row = repmat(Ab(i,:),n-i,1);

Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;

end

A=Ab(:,1:n);

b=Ab(:,n+1);

x=zeros(n,1);

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);

end

关于代码中的l.*rep_Ab_row，l是一个m行1列的系数矩阵，rep_Ab_row是一个m行n+1列的矩阵，这里有点像numpy里的广播机制，前者每一行乘到后者每一行的每个元素上。

比如

>> a

a =

1

2

3

>> b

b =

1 2

3 4

5 6

>> a.*b

ans =

1 2

6 8

15 18

列主元素法

高斯消元法虽然简单，但是要确保消元过程中

a

i

,

i

(

i

)

≠

0

a^{(i)}_{i,i}\neq0

ai,i(i)​​=0，因为要作为除数，而且当它比较小的时候，作为除数会导致结果特别大，使得计算过程中的误差变大。

列主元素法就是在每一次消元之前，先对A进行初等行变换，即在这一列中，选取绝对值最大元素所在的行，然后和当前的第一行交换，可以使得误差较小。

%列主元素法

function x=lzys(A,b)

[rows,cols]=size(A);

if rows~=cols

error("rows != cols")

end

n=rows;

Ab=[A,b];

%消元

for i=1:n-1

%先找列主元素

[M,idx] = max(abs(Ab(i:end,i)));

%交换

idx = idx + i-1;%max得到的idx=1时，指的是第i行

Ab([i,idx],:)=Ab([idx,i],:);%去掉这个分号，可以选择列主元素的结果

l= Ab((i+1):end,i)/Ab(i,i);

rep_Ab_row = repmat(Ab(i,:),n-i,1);

Ab((i+1):end,:) = Ab((i+1):end,:)-l.*rep_Ab_row;

end

A=Ab(:,1:n);

b=Ab(:,n+1);

x=zeros(n,1);

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

x(i)=(b(i)-sum(A(i,i+1:end).*(x(i+1:end)')))/A(i,i);

end

标签：Ab,end,rows,高斯消,array,列主,ldots,元法,vdots

来源： https://blog.csdn.net/gamebyc/article/details/114777230