Questo sito utilizza cookies, anche di terze parti, per mostrare pubblicità e servizi in linea con il tuo account. Leggi l'informativa sui cookies.
Username: Password: oppure
C/C++ - Risoluzione sistemi lineari: metodi diretti
Forum - C/C++ - Risoluzione sistemi lineari: metodi diretti

Avatar
Puffetta (Normal User)
Rookie


Messaggi: 21
Iscritto: 29/11/2009

Segnala al moderatore
Postato alle 17:30
Venerdì, 25/05/2012
Ciao a tuttti!
ho scritto questo programma che deve risolvere un sistema lineare del tipo Ax=b con il metodo di gauss o di cholesky(metodo scelto dall'utente) ho però un grosso problema: il metodo di gauss funziona perfettamente(a quanto pare), mentre se uso il metodo di cholesky il programma va in loop. mi potreste aiutare a trovare l'errore? io non riesco a vederlo:-(

grazie a tutti!!!!

Codice sorgente - presumibilmente C

  1. /*Risoluzione di sistemi lineari usando metodi diretti quali metodo di eliminazione di gauss, cholesky e thomas*/
  2.  
  3. #include<stdio.h>
  4. #include<stdlib.h>
  5. #include<math.h>
  6. #include<time.h>
  7.  
  8. typedef double matrice[300][300];
  9. typedef double vettore[300];
  10.  
  11.  
  12. int sceglimatrice();
  13. void creamatrice(int mat, matrice A, vettore b, int n);
  14. double U();
  15. void prodotto(matrice A, matrice B, int n);
  16. int sceglimetodo(int mat);
  17. void gauss(matrice A, vettore b, vettore x, int n);
  18. void riduci(matrice, vettore, int);
  19. int cercapivot(vettore, int, int);
  20. void scambiarighe(matrice, int, int, int);
  21. void cholesky(matrice A, vettore b, vettore x, int n);
  22. void risolvi_sup(matrice A, vettore b,vettore x, int n);
  23. void risolvi_inf(matrice A, vettore b,vettore x, int n);
  24.  
  25.  
  26. int main()
  27. {
  28.       int n, metodo, mat;
  29.       matrice A;
  30.       vettore x, b;
  31.      
  32.       printf("Risolviamo il sistema lineare Ax=b.\nInseriamo tutti i dati riguardanti la matrice A\n");
  33.      
  34.       mat=sceglimatrice();
  35.  
  36.       printf("\nQuante righe e colonne ha la matrice A? n=");
  37.       scanf("%d", &n);
  38.      
  39.       creamatrice(mat, A, b, n);
  40.      
  41.       metodo=sceglimetodo(mat);
  42.       /*
  43.       for(int i=1; i<=n; i++)
  44.          {
  45.             for(int j=1; j<=n; j++)
  46.             printf("  %lf  ", A[i][j]);
  47.  
  48.                 printf("\n");
  49.           }
  50.  
  51.      
  52.       for(int i=1; i<=n; i++) printf("\t%lf", b[i]);
  53.      
  54.       printf("\nmetodo=%d\n", metodo);
  55.       */
  56.      
  57.       switch(metodo)
  58.           {
  59.             case 1:  gauss(A, b, x, n);
  60.                      break;
  61.              
  62.             case 2:  cholesky(A, b, x, n);
  63.                      break;
  64.           }
  65.      
  66.      printf("\n\n Le soluzioni del sistema sono:");
  67.      
  68.      for(int i=1; i<=n; i++) printf("\t%lf", x[i]);
  69.  
  70.      system("PAUSE");
  71.      return 0;
  72. }
  73.  
  74.  
  75. double U()
  76. {
  77.        return (double)rand()/RAND_MAX;
  78. }
  79.  
  80.  
  81. int sceglimatrice()
  82. {
  83.      int mat=0;
  84.  
  85.      printf("Scegliere il tipo da matrice da utilizzare:\n");
  86.  
  87.      printf("\nMATRICE RANDOM\n\n");
  88.      printf(" 1) Matrice randomica generica\n\n");
  89.      printf(" 2) Matrice randomica simmetrica e definita positiva\n\n");
  90.      printf(" 3) Matrice randomica tridiagonale\n\n");
  91.      printf(" 4) Matrice randomica diagonale dominante\n");
  92.      printf("\nMATRICE INSERITA DA TASTIERA\n\n");
  93.      printf(" 5) Matrice generica\n\n");
  94.      printf(" 6) Matrice simmetrica e definita positiva\n\n");
  95.      printf(" 7) Matrice tridiagonale\n\n");
  96.      printf(" 8) Matrice diagonale dominante\n\n");
  97.      
  98.      printf("Quale tipo di matrice devo usare? ");
  99.      scanf("%d", &mat);    
  100.  
  101.      while((mat<=0)||(mat>8))
  102.        {
  103.           printf("\nIl valore inserito non e' corretto. Inserire nuovamente il valore  ");
  104.           scanf("%d", &mat);
  105.        }
  106.      
  107.      return mat;
  108. }
  109.  
  110. //Questa funzione crea la matrice in base alle esigenze dell'utente
  111. void creamatrice(int mat, matrice A, vettore b, int n)
  112. {
  113.    srand(time(0));
  114.    
  115.    switch(mat)    
  116.      {
  117.         //Matrice randomica generica e vettore randomico
  118.         case 1:
  119.              {                    
  120.                        for(int i=1; i<=n; i++)
  121.                            {
  122.                       for(int j=1; j<=n; j++)
  123.                                                A[i][j]=U();//matrice
  124.                                            
  125.                                   b[i]=U();//vettore
  126.                    }
  127.              }
  128.              break;
  129.        
  130.         //Matrice randomica simmetrica e definita positiva e vettore randomico  
  131.         case 2:
  132.              {
  133.                 matrice B;
  134.                
  135.                 for(int i=1; i<=n; i++)
  136.                                 {
  137.                       for(int j=i; j<=n; j++)
  138.                                                 {
  139.                           B[i][j]=U();
  140.                           B[j][i]=B[i][j];
  141.                                                 }
  142.                      
  143.                              
  144.                       b[i]=U();
  145.                      }
  146.                
  147.                 //Uso la matrice B per creare la matrice A richiesta facendo A=(Bt)*B con B simmetrica
  148.                 prodotto(A, B, n);  
  149.              }
  150.              break;
  151.        
  152.         //Matrice randomica tridiagonale e vettore randomico            
  153.         case 3:
  154.              {      
  155.                for(int i=1; i<=n; i++)
  156.                       {
  157.                                           A[i][i]=U();
  158.                        
  159.                       if(i!=n)
  160.                            A[i+1][i]=U();
  161.                          
  162.                       if(i!=1)
  163.                            A[i-1][i]=U();
  164.                          
  165.                       b[i]=U();
  166.                   }
  167.                         }
  168.                         break;
  169.            
  170.         case 4:
  171.              {
  172.                   //diagonale dominante!!!!!!
  173.              }
  174.              break;
  175.      }
  176.    
  177.    //Matrici inserite da tastiera. uso dunque lo scanf per acquisire le componenti
  178.    if((mat==5)||(mat==6)||(mat==7)||(mat==8))
  179.       {
  180.          printf("\nInserire i coefficienti della matrice A:");
  181.          for(int i=1; i<=n; i++)
  182.                 for(int j=1; j<=n; j++)
  183.                         scanf("%lf", &(A[i][j]));
  184.          
  185.          printf("\nInserire il vettore dei termini noti");
  186.          for(int i=1; i<=n; i++)
  187.              {
  188.                 printf("Inserire il termine noto b[%d]=", i);
  189.                 scanf("%lf", &b[i]);
  190.              }                                                        
  191.       }
  192. }
  193.  
  194.  
  195. int sceglimetodo(int mat)
  196. {
  197.   int metodo=0;
  198.    
  199.   printf("\nQuale metodo devo usare per risolvere il sistema lineare quadrato Ax=b?\n1.Gauss   2.Cholesky   3.Thomas\n");
  200.   scanf("%d", &metodo);
  201.  
  202.   while((metodo<=0)||(metodo>4))
  203.         {
  204.        printf("\n Il valore inserito non e' corretto!!! Scegliere di nuovo il metodo ");
  205.        scanf("%d", &metodo);
  206.         }
  207.  
  208.   if((metodo==2)&&((mat==2)||(mat==6)))
  209.                                    return metodo;
  210.          
  211.   if((metodo==3)&&((mat==3)||(mat==7)))
  212.                                    return metodo;
  213.  
  214.   while((metodo!=1)&&(mat!=1)&&(mat!=5))
  215.        {
  216.           printf("\nIl metodo scelto non puo' essere usato sulla matrice A.\nscegliere un altro metodo  ");
  217.           scanf("%d", &metodo);
  218.          
  219.           if((metodo==2)&&((mat==2)||(mat==6)))
  220.                                          break;
  221.          
  222.           if((metodo==3)&&((mat==3)||(mat==7)))
  223.                                           break;                                                                
  224.       }
  225.      
  226.   return metodo;
  227. }
  228.  
  229.  
  230. void gauss(matrice A, vettore b, vettore x, int n)
  231. {
  232.   riduci(A, b, n);
  233.      
  234.   risolvi_sup(A, b, x, n);
  235. }
  236.  
  237.  
  238. void riduci(matrice A, vettore b, int n)
  239. {
  240.      double mik; vettore kcolonna; int kmax;
  241.      
  242.      for(int k=1; k<=n-1; k++)
  243.              {
  244.                   for(int h=k; h<=n; h++)
  245.                           {
  246.                                kcolonna[h]=A[h][k];
  247.                                
  248.                                kmax=cercapivot(kcolonna, k, n);
  249.                                
  250.                                if(kmax!=k)    scambiarighe( A, k, kmax, n);                              
  251.                           }
  252.              
  253.                   for(int i=k+1; i<=n; i++)
  254.                      {
  255.                             mik=A[i][k]/A[k][k];
  256.                            
  257.                             for(int j=k+1; j<=n; j++)       A[i][j]=A[i][j]-mik*A[k][j];
  258.                            
  259.                             b[i]=b[i]-mik*b[k];
  260.                      }
  261.              }
  262.      return;    
  263. }
  264.  
  265. int cercapivot(vettore kcolonna, int k, int n)
  266. {
  267.     double maxmodulo=fabs(kcolonna[k]);
  268.    
  269.     int imax=k;
  270.    
  271.     for(int i=k+1; i<=n; i++)
  272.                    if(!maxmodulo)        exit (1);
  273.            
  274.     return imax;
  275. }
  276.  
  277.  
  278. void scambiarighe(matrice A, int k, int kmax, int n)
  279. {
  280.      double aux;
  281.      
  282.      for(int j=k; j<=n; j++)
  283.              {
  284.                   aux=A[j][k];
  285.      
  286.                   A[j][k]=A[j][kmax];
  287.      
  288.                   A[j][kmax]=aux;
  289.              }
  290.      return;
  291. }
  292.  
  293.  
  294. void cholesky(matrice A, vettore b, vettore x, int n)
  295. {
  296.      printf(" ********** ");
  297.      matrice H, HT;
  298.      vettore y;
  299.      
  300.      printf(" ++++++++++  ");
  301.      
  302.      //creo la matrice H
  303.      for(int i=1; i<=n; i++)
  304.        for(int j=1; j<=n; j++)
  305.            {
  306.               if(i==j==1)
  307.                  H[i][j]=sqrt(A[i][i]);
  308.                  
  309.               if(j<i)
  310.                   H[i][j]=0;
  311.                  
  312.               if(j>i)
  313.                 {
  314.                   double somma=0;
  315.                  
  316.                   for(int k=1; k<=j-1; k++)
  317.                                 somma+=H[i][k]*H[j][k];
  318.                                
  319.                   H[i][j]=(A[i][j]-somma)/H[j][j];    
  320.                 }  
  321.              
  322.               if(j==i)
  323.                 {
  324.                   double somma=0;
  325.                  
  326.                   for(int k=1; k<=i-1; k++)
  327.                                 somma+=pow(H[i][k], 2);
  328.                                
  329.                   H[i][j]=A[i][i]-somma;    
  330.                 }    
  331.            }    
  332.      
  333.      for(int i=1; i<=n; i++)
  334.         for(int j=1; j<=n; j++)
  335.                            HT[i][j]=H[j][i];
  336.                            
  337.      risolvi_inf(HT, b, y, n);
  338.      
  339.      risolvi_sup(H, y, x, n);
  340. }
  341.  
  342. //Questa funzione mi risolve un sistema la cui matrice e' triangolare inferiore
  343. void risolvi_inf(matrice A, vettore b,vettore x, int n)
  344. {    
  345.      x[1]=b[1]/A[1][1];
  346.      
  347.      for(int i=2; i<=n; i++)
  348.              {
  349.                 for(int j=1; j<=i; j++)        b[i]=b[i]-A[i][j]*x[j];
  350.                    
  351.                 x[i]=b[i]/A[i][i];
  352.              }
  353.      return;  
  354. }
  355.  
  356.  
  357. //Questa funzione mi risolve un sistema la cui matrice e' triangolare superiore
  358. void risolvi_sup(matrice A, vettore b,vettore x, int n)
  359. {    
  360.      x[n]=b[n]/A[n][n];
  361.      
  362.      for(int i=n-1; i>=1; i--)
  363.              {      
  364.                     for(int j=i+1; j<=n; j++)        b[i]=b[i]-A[i][j]*x[j];
  365.                    
  366.                     x[i]=b[i]/A[i][i];
  367.              }
  368.      return;  
  369. }
  370.  
  371.  
  372. //Questa matrice mi fa il prodotto (Bt)*B e pone il risultato della nuova matrice A
  373. void prodotto(matrice A, matrice B, int n)
  374. {
  375.    for(int i=1; i<=n; i++)
  376.        for(int j=1; j<=n; j++)
  377.              {
  378.                A[i][j]=0;
  379.                
  380.                for(int k=1; k<=n; k++)
  381.                        A[i][j]=A[i][j]+B[i][k]*B[j][k];
  382.              }
  383.    return;
  384. }


PM Quote