迭代法
时间函数有问题 ,当然代码也不是最好的。
这几天学习了三种迭代法:雅克比迭代法,高斯-赛德尔迭代法,超松弛迭代法;对方程组求解。
例如:试分别用雅克比迭代法,高斯-赛德尔迭代法,超松弛迭代法(取ω=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 }