防灾数学分享 http://blog.sciencenet.cn/u/fzmath 防灾科技学院数学教研室

博文

能显示计算过程的列主元高斯消元法的MATLAB程序

已有 5584 次阅读 2018-10-25 23:17 |系统分类:教学心得

 为了和教材上列主元高斯消元法的求解过程一致,编写了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,结果和教材完全一样。



https://blog.sciencenet.cn/blog-292361-1142815.html

上一篇:LATEX中chemarrow包突然工作不正常的解决方法——只对本案例
下一篇:红玉兰
收藏 IP: 124.238.133.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 16:42

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部