第一大类是直接法，理论基础是高等代数里的矩阵和行列式，也是众所周知的消元法。

x+y=a
x-y=b

小学就开始会解了，所以消元法在每个人心里都很明白，计算机算法也相对简单和容易理解，有完全主元消去法，列主元消去法，高斯若当消去法（消去对角线下方元素的同时消去上方元素）等。

高斯列主元素消去法实现：

//高期列主元素消去法，ab为{A,b}增广矩阵
void gauss1(CMatrix &ab)
{
 int h,w;
 ab.size(h,w);
 if(h+1!=w)//要求n阶方阵
  return;
 int n=h;
 int i,j;
 for(i=0;i<n;++i)
 {
  //从a[i,i]到a[n,i]找出最大元素所在行
  int max=i;//max指向最大列主元素所在行
  for(j=i+1;j<h;++j)
  {
   if(fabs(ab.elem(j,i))>fabs(ab.elem(max,i)))
    max=j;
  }
  ab.swap(i,max);//交换行
  if(ab.elem(i,i)==0.0)//det=0，计算停止
   return;
  //消元
  for(j=i+1;j<h;++j)
  {
   ab.pt(i,j,-(ab.elem(j,i)/ab.elem(i,i)));
  }
  //将a[i,i]化成1.0
  ab.mul(i,1.0/ab.elem(i,i));
 }
 for(i=n-1;i>=0;--i)
 {
  //消第i列上三角元素
  for(j=0;j<i;++j)
   ab.pt(i,j,-ab.elem(j,i));
 }
}

最后还有三角分解的方法，直接三角分解（解相同系数的方程组Ax=(b1,b2,...)），平方根法（对称正定方程组），追赶法（对角占优的三对角线方程组）。

追赶法实现：

void zhuiganfa(double const *a,double const *b,double const *c, // [in] 方程的三对角
        double const *f, // [in] f向量
        double *x, // [out] 方程的解
        int n, // [in] 方阵的阶
        double *beta=NULL // [in,none]
        )
{
 int i,tag=0;
 if(!beta)
 {
  beta=new double[n-1];
  tag=1;
 }
 double *y=x;
 beta[0]=c[0]/b[0];
 for(i=1;i<n-1;++i)
  beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);
 y[0]=f[0]/b[0];
 for(i=1;i<n;++i)
  y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]); // 追
 for(i=n-2;i>=0;--i)
  x[i]=y[i]-beta[i]*x[i+1]; // 赶
 if(tag)
  delete[] beta;
}

第二大类就是迭代法，有雅可比迭代，高斯-塞德尔迭代法，超松弛迭代法等。

高斯-塞德尔迭代法实现：

int gauss_seidel(CMatrix &ab, // [in] 方程的增广矩阵 A+b
    CMatrix &x, // [in,out] 迭代方程的解
    double eps // [in] 精度控制
    )
{
 int w,h,wx,hx;
 ab.size(h,w); // 得到行数和列数
 x.size(hx,wx);
 assert(wx==1 && h==hx);
 double dx,maxdx;
 int cnt=0; // 统计迭代次数作为函数返回值
 do
 {
  maxdx=0.0;
  for(int k=0;k<h;++k)
  {
   double s=0.0;
   for(int i=0;i<h;++i)
   {
    if(k!=i)
     s+=ab.elem(k,i)*x.elem(i,0);
   }
   double rst=(ab.elem(k,h)-s)/ab.elem(k,k); // 迭代计算x(k)
   if((dx=fabs(rst-x.elem(k,0)))>maxdx)
    maxdx=dx;
   x.elem(k,0)=rst;
  }
  ++cnt;
 }
 while(maxdx>eps);
 return cnt;
}

超松弛迭代法实现：

int sor(CMatrix &ab, // [in] 方程的增广矩阵 A+b
 CMatrix &x, // [in,out] 迭代方程的解
 double omg, // [in] 松弛因子
 double eps // [in] 精度控制
 )
{
 int w,h,wx,hx;
 ab.size(h,w); // 得到行数和列数
 x.size(hx,wx);
 assert(wx==1 && h==hx);
 double dx,maxdx;
 int cnt=0; // 统计迭代次数作为函数返回值
 do
 {
  maxdx=0.0;
  for(int k=0;k<h;++k)
  {
   double s=0.0;
   for(int i=0;i<h;++i)
   {
    if(k!=i)
     s+=ab.elem(k,i)*x.elem(i,0);
   }
   double rst=(ab.elem(k,h)-s)/ab.elem(k,k); // 迭代计算x~(k)
   rst=(1-omg)*x.elem(k,0)+omg*rst; // 加权平均
   if((dx=fabs(rst-x.elem(k,0)))>maxdx)
    maxdx=dx;
   x.elem(k,0)=rst;
  }
  ++cnt;
 }
 while(maxdx>eps);
 return cnt;
}

消元法是确定的解法，但在计算机里由于浮点数存在误差，所以实际上也是近似解，对于一些病态方程仍然无从下手，其时间复杂度在O（n3），对一些特殊方程的特殊解法除外，比如追赶法就是O（n）的线性时间复杂度，而迭代法是由一组初值，进行多次迭代收敛到近似解的过程，所以有收敛性问题，在一次迭代中需要一次矩阵乘法的运算量，是O（n2），所以迭代法至少是O（n2）以上的方法，而收敛的速度取决于迭代的次数，这是收敛速度的问题。

迭代法在解一些超大型方程组时是很有效的方法，比如成百上千的未知量，你要说现实中哪有这么大的方程，有的，比如椭圆方程最简单的五点格式，网格中有多少内点就有多少未知量，一个差分网格当然是越多越精确，而这个时候O（n3）的直接法就显得像蜗牛一样，是迭代法大展拳脚的时候。

超松弛迭代法有一个松弛因子，它的研究中心就是怎么求最佳松弛因子，现实中拿到一个方程，你总不能说，我先试一下，慢慢试，试出最佳因子，你试第一个的时候，结果就已经出来了，那你还要松弛因子干嘛，你要的是方程的解。最好的情况是，在迭代过程中，松弛因子会自动适应，解在逼近方程的解，而松弛因子也在逼近最佳松弛因子。

最后一个测试的简单程序：

#include "stdafx.h"
#include <cmath>
#include <cassert>
#include "cmatrix.h"
int main(int argc, char* argv[])
{
 double a[]={0,1,1},b[]={3,3,3},c[]={2,2,0},f[]={1,2,3},x[3];
 zhuiganfa(a,b,c,f,x,3);
 printf("x={%lf,%lf,%lf}\n",x[0],x[1],x[2]);
 double data[]={3,2,0,1,
  1,3,2,2,
  0,1,3,3};
 CMatrix ab(data,3,4),mx(3,1);
 int cnt;
 cnt=gauss_seidel(ab,mx,1e-5);
 //cnt=sor(ab,mx,1.2,1e-5);
 printf("x={%lf,%lf,%lf}\ncnt=%d\n",mx.elem(0,0),mx.elem(1,0),mx.elem(2,0),cnt);
 return 0;
}

rickone 2007/5/23

评论