雅可比(Jacobi)迭代算法的C++实现

admin

1. /*
   
2. 　　-----------------------------------------------
   
3. 　　假设有如下方程组：
   
4. 　　Ax=b
   
5. 　　用Jacobi迭代法求解方程组的解
   
6. 　　方法：将A分裂为A=D-L-U，等价的迭代方程组为x=Bx+f。
   
7. 　　有关算法的详细说明，参看http://www.loujing.com/mywork/c++/project/Jacobi.pdf
   
8. 　　-----------------------------------------------
   
9. 　　*/
   
10. 　　#include <iostream.h>
    
11. 　　#include <stdlib.h>
    
12. 　　#include <math.h>
    
13. 　　double* allocMem(int ); //分配内存空间函数
    
14. 　　void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值
    
15. 　　void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比迭代公式求解x的值
    
16. 　　void main()
    
17. 　　{
    
18. 　　short matrixNum; //矩阵的行数(列数)
    
19. 　　double *matrixA; //矩阵A，初始系数矩阵
    
20. 　　double *matrixD; //矩阵D为A中的主对角阵
    
21. 　　double *matrixL; //矩阵L为A中的下三角阵
    
22. 　　double *matrixU; //矩阵U为A中的上三角阵
    
23. 　　double *B; //矩阵B为雅可比方法迭代矩阵
    
24. 　　double *f; //矩阵f为中间的过渡的矩阵
    
25. 　　double *x; //x为一维数组，存放结果
    
26. 　　double *xk; //xk为一维数组，用来在迭代中使用
    
27. 　　double *b; //b为一维数组，存放方程组右边系数
    
28. 　　int i,j,k;
    
29. 　　cout<<"<<请输入矩阵的行数（列数与行数一致）>>：";
    
30. 　　cin>>matrixNum;
    
31. 　　//分别为A、D、L、U、B、f、x、b分配内存空间
    
32. 　　matrixA=allocMem(matrixNum*matrixNum);
    
33. 　　matrixD=allocMem(matrixNum*matrixNum);
    
34. 　　matrixL=allocMem(matrixNum*matrixNum);
    
35. 　　matrixU=allocMem(matrixNum*matrixNum);
    
36. 　　B=allocMem(matrixNum*matrixNum);
    
37. 　　f=allocMem(matrixNum);
    
38. 　　x=allocMem(matrixNum);
    
39. 　　xk=allocMem(matrixNum);
    
40. 　　b=allocMem(matrixNum);
    
41. 　　//输入系数矩阵各元素值
    
42. 　　cout<<endl<<endl<<endl<<"<<请输入矩阵中各元素值（为 "<<matrixNum<<"*"
    
43. 　　 <<matrixNum<<"，共计 "<<matrixNum*matrixNum<<" 个元素）"<<">>："<<endl<<endl;
    
44. 　　for(i=0;i<matrixNum;i++)
    
45. 　　{
    
46. 　　 cout<<"请输入矩阵中第 "<<i+1<<" 行的 "<<matrixNum<<" 个元素：";
    
47. 　　 for(j=0;j<matrixNum;j++)
    
48. 　　 cin>>*(matrixA+i*matrixNum+j);
    
49. 　　}
    
50. 　　//输入方程组右边系数b的各元素值
    
51. 　　cout<<endl<<endl<<endl<<"<<请输入方程组右边系数各元素值，共计 "<<matrixNum<<
    
52. 　　 " 个"<<">>："<<endl<<endl;
    
53. 　　for(i=0;i<matrixNum;i++)
    
54. 　　 cin>>*(b+i);
    
55. 　　/* 下面将A分裂为A=D-L-U */
    
56. 　　//首先将D、L、U做初始化工作
    
57. 　　for(i=0;i<matrixNum;i++)
    
58. 　　 for(j=0;j<matrixNum;j++)
    
59. 　　 *(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;
    
60. 　　//D、L、U分别得到A的主对角线、下三角和上三角；其中D取逆矩阵、L和U各元素取相反数
    
61. 　　for(i=0;i<matrixNum;i++)
    
62. 　　 for(j=0;j<matrixNum;j++)
    
63. 　　 if(i==j&&*(matrixA+i*matrixNum+j)) *(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));
    
64. 　　 else if(i>j) *(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
    
65. 　　 else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
    
66. 　　//求B矩阵中的元素
    
67. 　　for(i=0;i<matrixNum;i++)
    
68. 　　 for(j=0;j<matrixNum;j++)
    
69. 　　 {
    
70. 　　 double temp=0;
    
71. 　　 for(k=0;k<matrixNum;k++)
    
72. 　　 temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
    
73. 　　 *(B+i*matrixNum+j)=temp;
    
74. 　　 }
    
75. 　　//求f中的元素
    
76. 　　for(i=0;i<matrixNum;i++)
    
77. 　　{
    
78. 　　 double temp=0;
    
79. 　　 for(j=0;j<matrixNum;j++)
    
80. 　　 temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
    
81. 　　 *(f+i)=temp;
    
82. 　　}
    
83. 　　/* 计算x的初始向量值 */
    
84. 　　GaussLineMain(matrixA,x,b,matrixNum);
    
85. 　　/* 利用雅可比迭代公式求解xk的值 */
    
86. 　　int JacobiTime;
    
87. 　　cout<<endl<<endl<<endl<<"<<雅可比迭代开始，请输入希望迭代的次数>>：";
    
88. 　　cin>>JacobiTime;
    
89. 　　while(JacobiTime<=0)
    
90. 　　{
    
91. 　　 cout<<"迭代次数必须大于0，请重新输入：";
    
92. 　　 cin>>JacobiTime;
    
93. 　　}
    
94. 　　Jacobi(x,xk,B,f,matrixNum,JacobiTime);
    
95. 　　//输出线性方程组的解 */
    
96. 　　cout<<endl<<endl<<endl<<"<<方程组运算结果如下>>"<<endl;
    
97. 　　cout.precision(20); //设置输出精度，以此比较不同迭代次数的结果
    
98. 　　for(i=0;i<matrixNum;i++)
    
99. 　　 cout<<"x"<<i+1<<" = "<<*(xk+i)<<endl;
    
100. 　　cout<<endl<<endl<<"求解过程结束..."<<endl<<endl;
     
101. 　　//释放掉所有动态分配的内存
     
102. 　　delete [] matrixA;
     
103. 　　delete [] matrixD;
     
104. 　　delete [] matrixL;
     
105. 　　delete [] matrixU;
     
106. 　　delete [] B;
     
107. 　　delete [] f;
     
108. 　　delete [] x;
     
109. 　　delete [] xk;
     
110. 　　delete [] b;
     
111. 　　}
     
112. 　　/*
     
113. 　　----------------------
     
114. 　　 分配内存空间函数
     
115. 　　----------------------
     
116. 　　*/
     
117. 　　double* allocMem(int num)
     
118. 　　{
     
119. 　　double *head;
     
120. 　　if((head=new double[num])==NULL)
     
121. 　　{
     
122. 　　 cout<<"内存空间分配失败，程序终止运行！"<<endl;
     
123. 　　 exit(0);
     
124. 　　}
     
125. 　　return head;
     
126. 　　}
     
127. 　　/*
     
128. 　　-----------------------------------------------
     
129. 　　 计算Ax=b中x的初始向量值，采用高斯列主元素消去法
     
130. 　　-----------------------------------------------
     
131. 　　*/
     
132. 　　void GaussLineMain(double* A,double* x,double* b,int num)
     
133. 　　{
     
134. 　　int i,j,k;
     
135. 　　//共处理num-1行
     
136. 　　for(i=0;i<num-1;i++)
     
137. 　　{
     
138. 　　 //首先每列选主元，即最大的一个
     
139. 　　 double lineMax=fabs(*(A+i*num+i));
     
140. 　　 int lineI=i;
     
141. 　　 for(j=i;j<num;j++)
     
142. 　　 if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;
     
143. 　　 //主元所在行和当前处理行做行交换，右系数b也随之交换
     
144. 　　 for(j=i;j<num;j++)
     
145. 　　 {
     
146. 　　 //A做交换
     
147. 　　 lineMax=*(A+i*num+j);
     
148. 　　 *(A+i*num+j)=*(A+lineI*num+j);
     
149. 　　 *(A+lineI*num+j)=lineMax;
     
150. 　　 //b中对应元素做交换
     
151. 　　 lineMax=*(b+i);
     
152. 　　 *(b+i)=*(b+lineI);
     
153. 　　 *(b+lineI)=lineMax;
     
154. 　　 }
     
155. 　　 if(*(A+i*num+i)==0) continue; //如果当前主元为0，本次循环结束
     
156. 　　 //将A变为上三角矩阵，同样b也随之变换
     
157. 　　 for(j=i+1;j<num;j++)
     
158. 　　 {
     
159. 　　 double temp=-*(A+j*num+i)/(*(A+i*num+i));
     
160. 　　 for(k=i;k<num;k++)
     
161. 　　 {
     
162. 　　 *(A+j*num+k)+=temp*(*(A+i*num+k));
     
163. 　　 }
     
164. 　　 *(b+j)+=temp*(*(b+i));
     
165. 　　 }
     
166. 　　}
     
167. 　　/* 验证Ax=b是否有唯一解，就是验证A的行列式是否为0；
     
168. 　　 如果|A|!=0，说明有唯一解
     
169. 　　*/
     
170. 　　double determinantA=1;
     
171. 　　for(i=0;i<num;i++)
     
172. 　　 determinantA*=*(A+i*num+i);
     
173. 　　if(determinantA==0)
     
174. 　　{
     
175. 　　 cout<<endl<<endl<<"通过计算,矩阵A的行列式为|A|=0，即A没有唯一解。\n程序退出..."<<endl<<endl;
     
176. 　　 exit(0);
     
177. 　　}
     
178. 　　/* 从最后一行开始，回代求解x的初始向量值 */
     
179. 　　for(i=num-1;i>=0;i--)
     
180. 　　{
     
181. 　　 for(j=num-1;j>i;j--)
     
182. 　　 *(b+i)-=*(A+i*num+j)*(*(x+j));
     
183. 　　 *(x+i)=*(b+i)/(*(A+i*num+i));
     
184. 　　}
     
185. 　　}
     
186. 　　/*
     
187. 　　--------------------------------------
     
188. 　　利用雅可比迭代公式求解x的精确值
     
189. 　　迭代公式为：xk=Bx+f
     
190. 　　--------------------------------------
     
191. 　　*/
     
192. 　　void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
     
193. 　　{
     
194. 　　int t=1,i,j;
     
195. 　　while(t<=time)
     
196. 　　{
     
197. 　　 for(i=0;i<num;i++)
     
198. 　　 {
     
199. 　　 double temp=0;
     
200. 　　 for(j=0;j<num;j++)
     
201. 　　 temp+=*(B+i*num+j)*(*(x+j));
     
202. 　　 *(xk+i)=temp+*(f+i);
     
203. 　　 }
     
204. 　　 //将xk赋值给x，准备下一次迭代
     
205. 　　 for(i=0;i<num;i++)
     
206. 　　 *(x+i)=*(xk+i);
     
207. 　　 t++;
     
208. 　　}
     
209. 　　}
     
210. 　　//程序编写结束

复制代码