본문 바로가기

IT/Numerical analysis

야코비(Jacobi) Method C언어로 구현하기


  1.  
  2. /* Jacobi Method */
  3.  
  4. #include <stdio.h>
  5. #include <conio.h>
  6. #include <stdlib.h>
  7. #include <math.h>
  8.  
  9. int it=0;
  10. double **a;
  11. double *b;
  12. double *s;
  13. #define         EPS             1e-5            // 최대 오차
  14. #define         ITOR    500             // 최대 반복 횟수
  15.  
  16. int n; // 행렬의 갯수
  17.  
  18. FILE *fp;
  19.  
  20.  
  21. double** matrix_malloc(int n)
  22. {
  23.  
  24.         double **a;
  25.  
  26.         int i,j;
  27.         a=(double **)malloc(n*sizeof(double));
  28.         for(i=0;i<n;i++)
  29.         a[i]=(double *)malloc(n*sizeof(double));
  30.          for(i=0;i<n;i++){
  31.                  for(j=0;j<n;j++){
  32.                          a[i][j] = 0.0;
  33.                  }
  34.          }
  35.  
  36.         return a;
  37.  
  38. }
  39. double * b_malloc(int n)
  40. {
  41.         double *b;
  42.  
  43.         b = (double*)malloc(n*sizeof(double));
  44.  
  45.         return b;
  46.  
  47. }
  48.  
  49. double * s_malloc(int n)
  50. {
  51.         double *s;
  52.  
  53.         s = (double*)malloc(n*sizeof(double));
  54.  
  55.         return s;
  56.  
  57. }
  58.  
  59.  
  60. void matrix_free(int n)
  61. {
  62.  
  63.         int i;
  64.          
  65.         for(i=0;i<n;i++)
  66.                 free(a[i]);
  67.  
  68.          free(a);
  69.  
  70. }
  71.  
  72. void b_free()
  73. {
  74.         free(b);
  75. }
  76.  
  77.  
  78.  
  79. void s_free()
  80. {
  81.         free(s);
  82. }
  83.  
  84. void Jacobi_Method(void)
  85. {
  86.         int i,j;
  87.     double a_norm, x_norm, w;
  88.     double *x, *x1;
  89.         int it =0;
  90.  
  91.         x = s_malloc(n);
  92.         x1 = s_malloc(n);
  93.  
  94.     for(i=0;i<n;i++)
  95.             x[i] = 0.0;
  96.                
  97.          printf("Iteration");
  98.         for(i=0;i<n;i++)
  99.                 printf("    x%d    ",i);
  100.         printf("\n");
  101.  
  102.  
  103.  
  104.           do{
  105.                 it++;
  106.                 for(j=0;j<n;j++)
  107.                                         x1[j]=x[j];
  108.  
  109.                 a_norm = 0.0;
  110.                 x_norm = 0.0;
  111.  
  112.              
  113.                                 for(i=0;i<n;i++) {
  114.                                         w = b[i];
  115.                                         for(j=0;j<n;j++)
  116.                                         {
  117.                                                 if(j!=i)-= a[i][j]*x1[j];
  118.                                         }
  119.  
  120.                         w /= a[i][i];
  121.                         a_norm += fabs(x[i]-w);
  122.                         x_norm += fabs(w);
  123.                         x[i] = w;
  124.                                 }
  125.                 printf("    %-2d ",it);
  126.                 for(i=0;i<n;i++)
  127.                         printf(" %9.5lf", x[i]);
  128.                 printf("\n");
  129.                 if( it > ITOR){
  130.                         printf("Iteration= %d    Check Error!!!!", it);
  131.                         exit(1);
  132.                 }
  133.         } while(a_norm/x_norm > EPS);
  134.  
  135.          free(x);
  136.          free(x1);
  137.  
  138. }
  139.  
  140.  
  141. int main(void)
  142. {
  143.         int i, j;
  144.  
  145.  
  146.        
  147.         if((fp=fopen("matrix.txt""r"))==NULL){
  148.                  printf("File error!!\n");
  149.                 exit(1);
  150.         }
  151.  
  152.         fscanf(fp, "%d"&n);
  153.  
  154.         a = matrix_malloc(n);
  155.         b = b_malloc(n);
  156.         s = s_malloc(n);
  157.  
  158. // 파일 입력
  159.         for(i=0;i<n;i++){
  160.                  for(j=0;j<n;j++){
  161.  
  162.                          fscanf(fp,"%lf",&a[i][j]);
  163.                          printf("%.2f ",a[i][j]);
  164.                  }
  165.                 fscanf(fp,"%lf",&b[i]);
  166.         printf("%.2f \n",b[i]);
  167.         }
  168.  
  169.  
  170.         Jacobi_Method();
  171.  
  172.  
  173.         matrix_free(n);
  174.  
  175.         b_free();
  176.         s_free();
  177.         fclose(fp);
  178.        
  179.         return 0;
  180. }