数值计算消去法 C语言编写

时间:2012-06-20 01:31:50
【文件属性】:
文件名称:数值计算消去法 C语言编写
文件大小:5KB
文件格式:TXT
更新时间:2012-06-20 01:31:50
消去法 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;ifabs(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=0;--i) { //消第i列上三角元素 for(j=0;j=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;kmaxdx) 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;kmaxdx) 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 #include #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; }

网友评论