假设矩阵A是一个n阶非奇异方阵,那么对于任意一个n维向量b, 线性方程组Ax=b有唯一的解.
考虑如下方程组
我们通常为了求解它,我们把-4乘以第一行加到第二行,-7乘以第一行加到第三行,然后在进行求解.这个过程就是高斯消元法.下面我们看看通过matlab怎么实现.
%第一次消元
for i = 2:n % 矩阵的维数
mul = A(i,1)/ A(1,1);%求出2到n行的第一元素与第一行第一个元素的倍数
A(i,:)= A(i,:)-mul*A(1,:);%2到n行每个元素都减去与第一行的倍数
b(i)=b(i)-mul*b(i);%对列向量b进行同样操作
end
%以上过程可以得到矩阵A第二行至n行第一列所有元素为0
下面需要重复执行这个过程,是的第二列,第三列,第n-1列元素对角线一下都是0
%第2-n-1次消元
for j = 1:n-1
for i = j+1:n % 矩阵的维数
mul = A(i,j)/ A(j,j);%求出2到n行的第一元素与第一行第一个元素的倍数
A(i,:)= A(i,:)-mul*A(j,:);%2到n行每个元素都减去与第一行的倍数
b(i)=b(i)-mul*b(i);%对列向量b进行同样操作
end
end
%以上过程可以得到矩阵A第二行至n行第一列所有元素为0
以上过程可能会出现问题.因为对角元有可能非常小,甚至为0. 下面需要考虑这种情况,即选主元素消元.
%第2-n-1次消元
for j = 1:n-1
[pp,k]= max(abs(j:n,j));%找到第j列最大的元素所在位置
tempA = A(j,:);
A(j,:)= A(k+j-1,:);
A(k+j-1,:)=tempA;
tempb = b(k);
b(j,:)= b(k+j-1,:);
b(k+j-1,:)=temp;%换主元
for i = j+1:n % 矩阵的维数
mul = A(i,j)/ A(j,j);%求出2到n行的第一元素与第一行第一个元素的倍数
A(i,:)= A(i,:)-mul*A(j,:);%2到n行每个元素都减去与第一行的倍数
b(i)=b(i)-mul*b(i);%对列向量b进行同样操作
end
end
%以上过程可以得到矩阵A第二行至n行第一列所有元素为0
以上就是高斯消元法,下面是完整代码
clear all;
clc;
A = [1 2 3;
3 2 4;
1 3 2];
b = [1;1;3];
Ab
x= GaussianSolver0(A,b)
function x=GaussianSolver0(A,b)
[n,~] = size(A);
x=zeros(n,1);
for j = 1:n-1
for i = j+1:n
mul = A(i,j)/A(j,j);
A(i,:)= A(i,:)-mul*A(j,:);
b(i)= b(i)-mul*b(j);
end
end
for i=n:-1:1
sum=0;
for j=n:-1:i+1
sum=sum+x(j)*A(i,j);
end
x(i)=(b(i)-sum)/A(i,i);
end
end
