matlab列主元高斯,科学网—能显示计算过程的列主元高斯消元法的MATLAB程序 - 王福昌的博文...

为了和教材上列主元高斯消元法的求解过程一致，编写了MATLAB程序，供大家参考，互相交流学习。

先写出函数

function x = gsxylzy(A,b)

% format rat;clc

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

% b = [4; 7 ;6];

M = [A,b]

n = length(b);

%%% 消元过程

for i = 1:n-1

%交换列主元

for k = i:n-1

ak = max(abs(M(k:n,k))); index = find(M(:,k)==ak);

if isempty(index)

index = find(M(:,k)==-ak);

end

temp = M(index,:);  M(index,:) = M(k,:); M(k,:) = temp;

end

fprintf(['第 ',num2str(i),' 次交换列主元\n'])

M

% 消元

for k = i:n-1

M(k+1,:) = M(k+1,:)-M(i,:)*M(k+1,i)/M(i,i);

end

fprintf(['第 ',num2str(i),' 列消元\n'])

M

end

%%% 回代过程

x = zeros(n,1);

x(n) = M(n,n+1)/M(n,n);

fprintf(['第 ',num2str(n),' 个变量 x(',num2str(n),')=%s\n'],num2str(x(n)));

for k = (n-1):-1:1

x(k) = (M(k,n+1)-M(k,k+1:n)*x(k+1:n))/M(k,k);

fprintf(['第 ',num2str(k),' 个变量 x(',num2str(k),')=% s\n'],num2str(x(k)));

end

主程序测试应用：

format rat;clc

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

b = [4; 7 ;6];

x = gsxylzy(A,b);

fprintf('方程组的解为\n')

x

运行结果为

M =

-3              2              6              4

10             -7              0              7

5             -1              5              6

第 1 次交换列主元

M =

10             -7              0              7

-3              2              6              4

5             -1              5              6

第 1 列消元

M =

10             -7              0              7

0             -1/10           6             61/10

0              5/2            5              5/2

第 2 次交换列主元

M =

10             -7              0              7

0              5/2            5              5/2

0             -1/10           6             61/10

第 2 列消元

M =

10             -7              0              7

0              5/2            5              5/2

0              0             31/5           31/5

第 3 个变量 x(3)=1

第 2 个变量 x(2)=-1

第 1 个变量 x(1)=2.6645e-16

方程组的解为

x =

1/3752999689475414 %这里严格说是0

-1

1

>>

如果不想出现这种分数近似，也可以用符号形式，则不会出现近似值.

format rat;clc

A = sym([-3 2 6;10 -7 0; 5 -1 5]);

b = sym([4; 7 ;6]);

x = gsxylzy(A,b);

fprintf('方程组的解为\n')

x

则运行结果为

M =

[ -3,  2, 6, 4]

[ 10, -7, 0, 7]

[  5, -1, 5, 6]

第 1 次交换列主元

M =

[ 10, -7, 0, 7]

[ -3,  2, 6, 4]

[  5, -1, 5, 6]

第 1 列消元

M =

[ 10,    -7, 0,     7]

[  0, -1/10, 6, 61/10]

[  0,   5/2, 5,   5/2]

第 2 次交换列主元

M =

[ 10,    -7, 0,     7]

[  0,   5/2, 5,   5/2]

[  0, -1/10, 6, 61/10]

第 2 列消元

M =

[ 10,  -7,    0,    7]

[  0, 5/2,    5,  5/2]

[  0,   0, 31/5, 31/5]

第 3 个变量 x(3)=1

第 2 个变量 x(2)=-1

第 1 个变量 x(1)=0

方程组的解为

x =

0

-1

1

>>

这里使用了符号变量命令 sym，结果和教材完全一样。

转载本文请联系原作者获取授权，同时请注明本文来自王福昌科学网博客。

链接地址：http://blog.sciencenet.cn/blog-292361-1142815.html

上一篇：LATEX中chemarrow包突然工作不正常的解决方法——只对本案例

下一篇：红玉兰