迭代法

 

时间函数有问题 ,当然代码也不是最好的。

 

这几天学习了三种迭代法:雅克比迭代法,高斯-赛德尔迭代法,超松弛迭代法;对方程组求解。

例如:试分别用雅克比迭代法,高斯-赛德尔迭代法,超松弛迭代法(取ω=1.15)解线性方程组

 

 当max ¦xi(k+1)-Xi(k)¦<10-5时迭代终止。方程组的精确解为X*=(1,-2,-1,3)T

一.雅克比迭代法的公式:Xi(k+1)=Xi(k)+(bi-∑aijXj(k))/aii。1≤j≤n,n为X的列数,k为迭代次数,a[4][4],b[4],X[4],三个矩阵。

二.高斯-赛德尔迭代法:1.X1(k+1)=(b1-a12X2(k)-a13X3(k)-......a1nXn(k))/a11

                               2.X2(k+1)=(b2-a21X1(k+1)-a22X2(k)-a23X3(k)-......-a2nXn(k))/a22;

                               .     .     .     .      .     .     .     .     .     .     .     .     .      . 

                               n.Xn(k+1)=(bn-an1X1(k+1)-an2X2(k+1)-an3X3(k+1)-......-annXn(k+1))/ann

三.超松弛迭代法:       Xi(k+1)=Xi+ω(bi-∑aijXj(k+1)-∑Xij(k))/aii

下面看一下代码:

  1 #include"stdio.h"
  2 #include"iostream.h"
  3 #include"math.h"
  4 #include "time.h"
  5 double Jacobi(double a[][4],double b[])
  6 {
  7     int q=0,i;
  8     clock_t start,finish;
  9     
 10     double c[4]={0},x[4]={0},t[4]={0};
 11     double d,duration;
 12     start=clock();
 13     while(1)
 14     {
 15         q++;
 16         t[0]=x[0]+(b[0]-a[0][0]*x[0]-a[0][1]*x[1]-a[0][2]*x[2]-a[0][3]*x[3])/a[0][0];
 17         t[1]=x[1]+(b[1]-a[1][0]*x[0]-a[1][1]*x[1]-a[1][2]*x[2]-a[1][3]*x[3])/a[1][1];
 18         t[2]=x[2]+(b[2]-a[2][0]*x[0]-a[2][1]*x[1]-a[2][2]*x[2]-a[2][3]*x[3])/a[2][2];
 19         t[3]=x[3]+(b[3]-a[3][0]*x[0]-a[3][1]*x[1]-a[3][2]*x[2]-a[3][3]*x[3])/a[3][3];
 20         for( i=0;i<4;i++)
 21             c[i]=fabs(t[i]-x[i]);
 22         d=c[0];
 23         for( i=1;i<4;i++)
 24             if(d<c[i])
 25                 d=c[i];
 26         for( i=0;i<4;i++)
 27             x[i]=t[i];
 28     if(d<0.00001)
 29     {
 30         cout<<"迭代次数为"<<q<<"次后的结果为:"<<endl;
 31         for(i=0;i<4;i++)
 32             printf("x[%d]=%.7lf\n",i,x[i]);
 33             //cout<<"x["<<i<<"]="<<x[i]<<endl;
 34         finish=clock();
 35         duration=(double)(finish-start);
 36         cout<<"耗时"<<duration<<"ms"<<endl;
 37         break;
 38     }    
 39     }
 40     return 0;
 41 }
 42 double Gauss_Seidel(double a[][4],double b[])
 43 {
 44     int q=0,i;
 45     clock_t start,finish;
 46     double c[4]={0},x[4]={0},t[4]={0};
 47     double d,duration;
 48     start=clock();
 49     while(1)
 50     {
 51         q++;
 52         t[0]=x[0]+(b[0]-a[0][0]*x[0]-a[0][1]*x[1]-a[0][2]*x[2]-a[0][3]*x[3])/a[0][0];
 53         t[1]=x[1]+(b[1]-a[1][0]*t[0]-a[1][1]*x[1]-a[1][2]*x[2]-a[1][3]*x[3])/a[1][1];
 54         t[2]=x[2]+(b[2]-a[2][0]*t[0]-a[2][1]*t[1]-a[2][2]*x[2]-a[2][3]*x[3])/a[2][2];
 55         t[3]=x[3]+(b[3]-a[3][0]*t[0]-a[3][1]*t[1]-a[3][2]*t[2]-a[3][3]*x[3])/a[3][3];
 56         for( i=0;i<4;i++)
 57             c[i]=fabs(t[i]-x[i]);
 58         d=c[0];
 59         for( i=1;i<4;i++)
 60             if(d<c[i])
 61                 d=c[i];
 62         for( i=0;i<4;i++)
 63             x[i]=t[i];
 64         if(d<0.00001)
 65         {
 66             cout<<"迭代次数为"<<q<<"次后的结果为:"<<endl;
 67             for(i=0;i<4;i++)
 68                 printf("x[%d]=%.7lf\n",i,x[i]);
 69                 //cout<<"x["<<i<<"]="<<x[i]<<endl;
 70             finish=clock();
 71             duration=(double)(finish-start);
 72             //printf("%.3lf\n",finish-start);
 73             cout<<"耗时"<<duration<<"ms"<<endl;
 74             break;
 75         }
 76     }
 77     return 0;
 78 }
 79 double SOR(double a[][4],double b[])
 80 {
 81     int q=0,i;
 82     clock_t start,finish;
 83     double c[4]={0},x[4]={0},t[4]={0};
 84     double d,duration;
 85     start=clock();
 86     while(1)
 87     {
 88         q++;
 89         t[0]=x[0]+(b[0]-a[0][0]*x[0]-a[0][1]*x[1]-a[0][2]*x[2]-a[0][3]*x[3])/a[0][0]*1.15;
 90         t[1]=x[1]+(b[1]-a[1][0]*t[0]-a[1][1]*x[1]-a[1][2]*x[2]-a[1][3]*x[3])/a[1][1]*1.15;
 91         t[2]=x[2]+(b[2]-a[2][0]*t[0]-a[2][1]*t[1]-a[2][2]*x[2]-a[2][3]*x[3])/a[2][2]*1.15;
 92         t[3]=x[3]+(b[3]-a[3][0]*t[0]-a[3][1]*t[1]-a[3][2]*t[2]-a[3][3]*x[3])/a[3][3]*1.15;
 93         for( i=0;i<4;i++)
 94         {
 95             c[i]=fabs(t[i]-x[i]);
 96         }
 97         d=c[0];
 98         for( i=1;i<4;i++)
 99         {
100             if(d<c[i])
101                 d=c[i];
102         }
103         for( i=0;i<4;i++)
104             x[i]=t[i];
105         if(d<0.00001)
106         {
107             cout<<"迭代次数为"<<q<<"次后的结果为:"<<endl;
108             for(i=0;i<4;i++)
109                 printf("x[%d]=%.7lf\n",i,x[i]);
110                 //cout<<"x["<<i<<"]="<<x[i]<<endl;
111             finish=clock();
112             duration=(double)(finish-start);
113             //printf("%.3lf\n",finish-start);
114             cout<<"耗时"<<duration<<"ms"<<endl;
115             break;
116         }
117     }
118     return 0;
119 }
120 void main()
121 {
122     int n;
123     double v[4][4]={5,1,-1,-2,
124                    2,8,1,3,
125                    1,-2,-4,-1,
126                    -1,3,2,7};
127     double f[4]={-2,-6,6,12};
128     double c[4]={0},x[4]={0},m[4]={0},t[4]={0};
129     printf("\n****************************");
130     printf("\n*  1.  雅克比迭代法        *");
131     printf("\n*  2.  高斯-塞德尔迭代法   *");
132     printf("\n*  3.  超松弛迭代法        *");
133     printf("\n*  0.    GAME OVER         *");
134     printf("\n****************************");
135     cout<<endl;
136     printf("请选择一种迭代方法:");
137     cin>>n;
138     while(n!=0)
139     {
140         switch(n)
141         {
142         case 0: break;
143         case 1: Jacobi(v,f);break;
144         case 2: Gauss_Seidel(v,f);break;
145         case 3: SOR(v,f);break;
146         default: cout<<"Error"<<endl;   
147         }
148         printf("请选择一种迭代方法:");
149         cin>>n;
150     }
151             
152 }

 

posted @ 2012-10-19 14:27  Marshalkk  阅读(4215)  评论(0编辑  收藏  举报