软件工程

本类阅读TOP10

·PHP4 + MYSQL + APACHE 在 WIN 系统下的安装、配置
·Linux 入门常用命令(1)
·Linux 入门常用命令(2)
·使用 DCPROMO/FORCEREMOVAL 命令强制将 Active Directory 域控制器降级
·DirectShow学习(八): CBaseRender类及相应Pin类的源代码分析
·基于ICE方式SIP信令穿透Symmetric NAT技术研究
·Windows 2003网络负载均衡的实现
·一网打尽Win十四种系统故障解决方法
·数百种 Windows 软件的免费替代品列表
·收藏---行百里半九十

分类导航
VC语言Delphi
VB语言ASP
PerlJava
Script数据库
其他语言游戏开发
文件格式网站制作
软件工程.NET开发
全主元Gauss-Jordan消元法(Blitz++库)

作者:未知 来源:月光软件站 加入时间:2005-2-28 月光软件站

全主元Gauss-Jordan消元法

 

 

 

 

Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法,该方法的优点是稳定,而全主元[1]法可能更加稳定,同时,它也存在着弱点——求解速度可能较慢。

 

Gauss-Jordan消元法主要特点是通过交换任意的行列,把矩阵A约化为单位矩阵,约化完成后,方程组右端项向量b即为解向量。我们知道,选择绝对值最大的元素作为主元是很好的办法,所以,全主元法目的就是在一个大的范围里面寻找主元,以达到比较高的精度。

 

下面,我使用Blitz++Array库,编写全主元Gauss-Jordan消元法。

 

Code

 

#include <blitz/array.h>

#include <cstdlib>

#include <algorithm>

#include <vector>

using namespace blitz;

 

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b)

{

     int n = A.rows(), m = b.cols();

     int irow, icol;

     vector<int> indexcol(n), indexrow(n), piv(n);

 

     for (int j=0; j<n; ++j)

         piv.at(j) = 0;

    

     //寻找绝对值最大的元素作为主元

     for (int i=0; i<n; ++i) {

         double big = 0.0;

 

         for (int j=0; j<n; ++j)

              if (piv.at(j) != 1)

                   for (int k=0; k<n; ++k) {

                       if (piv.at(k) == 0) {

                            if (abs(A(j, k)) >= big) {

                                 big = abs(A(j, k));

                                 irow = j;

                                 icol = k;

                                 if (irow == icol) break;

                            }

                       }

                   }

 

         ++piv.at(icol);

        

         //进行行交换,把主元放在对角线位置上,列进行假交换,

         //使用向量indexrowindexcol记录主元位置,    

         //这样就可以得到最终次序是正确的解向量。

         if (irow != icol) {

              for (int l=0; l<n; ++l)

                   swap(A(irow, l), A(icol, l));

 

              for (int l=0; l<m; ++l)

                   swap(b(irow, l), b(icol, l));

         }

 

         indexrow.at(i) = irow;

         indexcol.at(i) = icol;

        

         try {

              double pivinv = 1.0 / A(icol, icol);

 

              for (int l=0; l<n; ++l)

                   A(icol, l) *= pivinv;

              for (int l=0; l<m; ++l)

                   b(icol, l) *= pivinv;

 

              //进行行约化

              for (int ll=0; ll<n; ++ll)

                   if (ll != icol) {

                       double dum = A(ll, icol);

 

                       for (int l=0; l<n; ++l)

                            A(ll, l) -= A(icol, l)*dum;

                       for (int l=0; l<m; ++l)

                            b(ll, l) -= b(icol, l)*dum;

                   }

         }

         catch (...) {

              cerr << "Singular Matrix";

         }

     }

}

 

int main()

{

     //测试矩阵

     Array<double, 2> A(3,3), b(3,1);

     A =  10,-19,-2,

         -20, 40, 1,

           1,  4, 5;

 

     b = 3,

         4,

         5;

    

     Gauss_Jordan(A, b);

    

     cout << "Solution = " << b <<endl;

}

 

Result

 

Solution = 3 x 1

[   4.41637

    2.35231

   -1.76512 ]

 

 

从代码的过程可以看出,矩阵A的逆在A中逐步构造,最终矩阵A演变成单位矩阵,解向量X也逐步替代右端项向量,且使用同一存储空间。

 

 

注释:[1]主元,又叫主元素,指用作除数的元素

 




相关文章

相关软件