Project1

标题: 【C语言发布】共轭梯度法解非齐次线性方程组（中） [打印本页]
作者: zaiy2863    时间: 2016-11-25 19:51
标题: 【C语言发布】共轭梯度法解非齐次线性方程组（中）
本帖最后由 zaiy2863 于 2016-11-25 23:34 编辑

还是@RyanBern 教我的方法。谢谢rb的说。
但是直到最后也没有解决的问题。
另外感谢在格式上的一些指导。
测试函数暂时先使用Hilbert矩阵。目前没有写读入函数。但是计算精度很不好。
共轭梯度法只限于正定对称的矩阵，但是可以很好地解决一些gauss法解决不了的问题。
最后 @⑨姐姐 我爱你。

C 代码复制
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4.  
  5. void Hilbrt(int n, double *a, double *b);
  6. void Triangle(int n, double *a, double *b);//三对角矩阵 
  7. void Triangle2(int n, int i, double *a, double b);//改良的三对角矩阵 
  8. void Read(); //Read函数没有写好。先以Hilbrt代替 
  9. void Print(int n, double *x);
  10. void Initialize(int n, double *x);
  11. void Grad(int n, double *g, double *A, double *b, double *x); 
  12. int JudgeG(int n, double *g); 
  13. void Dirct(int n,double a, double *d, double *g); 
  14. double PramterT(int n, double dAd, double *d, double *g);
  15. double CoefA(int n, double dAd, double *d, double *A, double *g);
  16. double DTAd(int n, double *d, double *A); 
  17. void NextX(int n, double *d, double *x, double t);
  18.  
  19. int main(){
  20.         int i,n;
  21.         int flag;
  22.         scanf("%d", &n);
  23.         double A[n*n], b[n];//存储录入的方程
  24.         //double A[n],b;
  25.         double x[n], d[n], g[n], t, a, dAd;//存储解和中间向量 
  26.         Hilbrt(n, A, b);//录入测试矩阵 
  27.         Initialize(n, x);//先给出一组x0，代入x向量中 
  28.         Grad(n, g, A, b, x);//先把x 代入g求出g0 
  29.         flag = JudgeG(n, g);//判断G的距离是不是等于0了 
  30.         //初始化d向量 得到d0 
  31.         for(i = 0; i<n; i++){
  32.                 d[i] = -g[i];
  33.         } 
  34.         a = 1;
  35.  
  36.         while(flag){
  37.                 dAd = DTAd(n, d, A); //计算准备工作 
  38.                 if(fabs(a) > 10E-15){
  39.                         t = PramterT(n, dAd, d, g);//计算t 
  40.                         NextX(n, d, x, t);//通过t计算得到下一组x
  41.                         Grad(n, g, A, b, x);//计算g值 现在是g_k+1 
  42.                         flag = JudgeG(n, g);//判断g是否为零 
  43.                         a = CoefA(n, dAd, d, A, g);//计算系数a 
  44.                         Dirct(n, a, d, g);//根据一组d_k算出d_k+1 
  45.                 }else{
  46.                         flag = 0;
  47.                 }
  48.         }
  49.         t = PramterT(n, dAd, d, g);//计算t 
  50.         NextX(n, d, x, t);
  51.         Print(n, x);//如果g为零的话，输出答案 
  52.         return 0;
  53. } 
  54.  
  55. void Initialize(int n, double *x){
  56.         int i;
  57.         x[0] = 1.0;
  58.         for(i = 1;i<n; i++){
  59.                 x[i] = 0.0;
  60.         }
  61. } 
  62.  
  63. void Grad(int n, double *g, double *A, double *b, double *x){
  64.         int i,j;
  65.         for(i = 0; i<n; i++){
  66.                 //Triangle2(n,i,a,b);
  67.                 g[i] = 0; 
  68.                 for(j = 0; j<n; j++){
  69.                         g[i] += A[i*n + j] * x[j];
  70.                         //g[i] += A[j] * x[j];
  71.                 }
  72.                 g[i] -= b[i];
  73.         }
  74. }
  75.  
  76. int JudgeG(int n, double *g){//判断g是不是0 
  77.         int i = 0;
  78.         double sum = 0.0;
  79.         for(i = 0; i<n; i++){
  80.                 sum += g[i]*g[i];
  81.         }
  82.         //printf("\n",sum);
  83.         if(fabs(sum) > 1.0E-15){
  84.                 return 1;
  85.         }else{
  86.                 return 0;
  87.         }
  88. }
  89.  
  90. double DTAd(int n, double *d, double *A){
  91.         double dAd = 0.0;
  92.         double dA[n]; 
  93.         int i,j; 
  94.         for(i = 0; i< n; i++){
  95.                 dA[i] = 0;
  96.                 for(j = 0; j<n; j++){
  97.                         dA[i] += d[j] * A[i*n + j];
  98.                 }
  99.         } 
  100.         for(i = 0; i<n; i++){
  101.                 dAd += dA[i] * d[i]; 
  102.         }
  103.         return dAd;
  104. }
  105.  
  106. double PramterT(int n, double dAd, double *d, double *g){
  107.         double t = 0.0;
  108.         int i,j;
  109.         for(i = 0; i< n; i++){
  110.                 t -= g[i] * d[i];
  111.         }
  112.         t /= dAd;
  113.         return t;
  114. }
  115.  
  116. double CoefA(int n, double dAd, double *d, double *A, double *g){
  117.         double a = 0.0;
  118.         double dA[n]; 
  119.         int i,j; 
  120.         for(i = 0; i< n; i++){
  121.                 dA[i] = 0;
  122.                 for(j = 0; j<n; j++){
  123.                         dA[i] += d[j] * A[i*n + j];
  124.                 }
  125.         } 
  126.         for(i = 0; i<n; i++){
  127.                 a += dA[i] * g[i]; 
  128.         }
  129.         a /= dAd;
  130.         return a;
  131. }
  132.  
  133. void NextX(int n, double *d, double *x, double t){
  134.         int i;
  135.         for(i = 0; i<n; i++){
  136.                 x[i] += t *d[i];
  137.         }
  138. }
  139.  
  140. void Dirct(int n,double a, double *d, double *g){
  141.         int i;
  142.         for(i = 0; i<n; i++){
  143.                 d[i] = -g[i] + a* d[i];
  144.         }//完成方向d的更新 
  145. }
  146.  
  147. void Print(int n, double *x){
  148.         int i;
  149.         for(i = 0; i<n; i++){
  150.                 printf("%.4lf\n", x[i]);
  151.         }
  152. }
  153.  
  154. void Hilbrt(int n, double *a, double *b){
  155.         //测试矩阵：录入n阶Hilbrt矩阵 
  156.         int i,j;
  157.         for(i = 0; i<n; i++){
  158.                         b[i] = 0.0;
  159.                 }
  160.         for (i = 0 ;i < n; i++){//第i行 第j列的矩阵 
  161.                 for(j = 0; j<n; j++){
  162.                         a[i*n + j] = 1.0 / (i+j+1);
  163.                         b[i] += a[i*n + j];
  164.                 }
  165.         } 
  166. }
  167.  
  168. void Triangle(int n, double *a, double *b){
  169.         //测试矩阵：录入三对角矩阵 
  170.         int i = 0;
  171.         for(i = 0; i< n*n; i++){
  172.             a[i] = 0.0;
  173.         }
  174.         for(i = 0; i< n; i++){
  175.             a[i*n + i] = 10.0;
  176.             if(i>0){
  177.                         a[i*n + i - 1] = 1.0;
  178.                 }
  179.             if(i < n-1){
  180.                     a[i*n + i + 1] = 1.0;
  181.                 }
  182.         }
  183.         b[0] = 11;
  184.         for(i = 1; i<n-1; i++){
  185.                 b[i] = 12;
  186.         } 
  187.         b[n-1] = 11;
  188. } 
  189. void Triangle2(int n, int i, double *a, double b){
  190.         a[i] = 10;
  191.         b[i] = 12;
  192.         if(i > 0){
  193.                 a[i - 1] = 1;
  194.         }
  195.         if(i < n-1){
  196.                 a[i + 1] = 1;
  197.         }
  198.         if(i == 0 || i == n-1){
  199.                 b = 11;
  200.         }
  201. }//改良的三对角矩阵
作者: RyanBern    时间: 2016-11-25 22:27
本帖最后由 RyanBern 于 2016-11-26 21:32 编辑

代码写得比上次有进步。但是对函数的划分怎么这么别扭呢。

@zaiy2863 明天或者后天我补一份样例吧。感兴趣的话务必看一下。

源码范例：https://github.com/RyanBernX/math/tree/master/CG.rbp/examples/C
说明书：https://ryanbernx.github.io/math/cg/index.html

最后，三对角的英文是 Tridiagonal，不是 Triangle 啦。

C 代码复制
  1. void Triangle2(int n, int i, double *a, double b){
  2.         a[i] = 10;
  3.         b[i] = 12;
  4.         if(i > 0){
  5.                 a[i - 1] = 1;
  6.         }
  7.         if(i < n-1){
  8.                 a[i + 1] = 1;
  9.         }
  10.         if(i == 0 || i == n-1){
  11.                 b = 11;
  12.         }
  13. }//改良的三对角矩阵

这一段代码有问题，请仔细看一下。
作者: 浮云半仙    时间: 2016-11-27 15:18
膜拜叶子姐姐！！
我试试去用民科算法之三分法去算x向量的极小值(雾)



欢迎光临 Project1 (https://rpg.blue/) Powered by Discuz! X3.1