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++ - Formule di quadratura interpolatorie
Forum - C/C++ - Formule di quadratura interpolatorie

Avatar
Puffetta (Normal User)
Rookie


Messaggi: 21
Iscritto: 29/11/2009

Segnala al moderatore
Postato alle 19:54
Sabato, 07/04/2012
Ciao a tutti!
ho scritto un programma che deve calcolare il valore dell'integrale definito di una certa funzione scelta dall'utente utilizzando uno dei tre metodi: rettangoli, trapezi e Simpson. c'è però un errore nel calcolo finale e io credo sia dovuto alla cattiva approssimazione. potreste aiutarmi a capire dov'è l'errore?

grazie mille!
vi posto il codice che ho scritto...

Codice sorgente - presumibilmente C

  1. //foglio di es. 4
  2.  
  3. #include<stdio.h>
  4. #include<stdlib.h>
  5. #include<math.h>
  6.  
  7. typedef double vettore [100];
  8.  
  9. int sceglif();
  10. int sceglimetodo();
  11. double funzione(int t, double d);
  12. void intervallo(double &a, double &b, int t);
  13. void nodi(double a, double b, vettore X, int m, int t);
  14. void puntimedi(vettore X, int m, vettore medio);
  15. double rettangoli(int t, int m, vettore X);
  16. double trapezi(int t, int m, vettore X);
  17. double simpson(int t, int m, vettore X);
  18.  
  19. int main()
  20. {
  21.    int t, q, m;
  22.    double a, b;
  23.    vettore X;
  24.    
  25.    t=sceglif();
  26.    
  27.    intervallo(a, b, t);
  28.    
  29.    q=sceglimetodo();
  30.    
  31.    printf("Inserire il numero di nodi per partizione uniforme dell'intervallo m=");
  32.    scanf("%d", &m);
  33.    
  34.    nodi(a, b, X, m, t);
  35.    
  36.    if(q==1)
  37.        printf("Il valore dell'integrale con il metodo dei rettangoli e' %.16lf\n", rettangoli(t, m,  X));
  38.    else if(q==2)      
  39.        printf("Il valore dell'integrale con il metodo dei trapezi e' %.16lf\n", trapezi(t, m, X));
  40.    
  41.    else if(q==3)
  42.        printf("Il valore dell'integrale con il metodo di C-S e' %.16lf\n", simpson(t, m, X));
  43.    
  44.    system("PAUSE");
  45.    return 0;  
  46. }
  47.  
  48.  
  49. int sceglif()
  50. {
  51.   int t=0;
  52.  
  53.   printf("\n Scegli la funzione da approssimare con il polinomio di Lagrange:\n\n");
  54.  
  55.   printf(" 1) integrale di sen(x) per x in [0;2pi]\n\n");
  56.   printf(" 2) integrale di x*e^(-x)*cos(2x) per x in [0;2pi]\n\n");
  57.   printf(" 3) integrale di 1/(1+x^2) per x in [-5;5]\n\n");
  58.   printf(" 4) integrale di x^5/2 per x in [0;1]\n\n");
  59.   printf(" 5) integrale di x*sin(x) per x in [-2;2]\n\n");
  60.   printf(" 6) integrale di mod(x) per x in [-1;0] oppure in [-1;1]\n\n");
  61.   printf(" 7) integrale di mod(sin(x)) per x in [0;2pi]\n\n");
  62.   printf(" 8) integrale di e(-x^2) per x in [-1;1]\n\n");
  63.   printf(" 9) integrale di sgn(x-pi) per x in [-3;5]\n\n");
  64.   printf(" 10) integrale di f(x) per x in [-2;3] dove f(x)=-x per x<=0 e x^2+1 per x>0\n\n");
  65.  
  66.   printf("Quale funzione devo usare? t=");
  67.   scanf("%d", &t);
  68.  
  69.   while(t>10)
  70.         {
  71.        printf("\n Il valore inserito non è corretto!!! Inserire un valore corretto di t t=");
  72.        scanf("%d", &t);
  73.         }
  74.  
  75.   return t;
  76. }    
  77.  
  78. int sceglimetodo()
  79. {
  80.     int q;
  81.    
  82.     printf("Quale metodo vuoi utilizzare?\n1. Metodo dei rettangoli\n2.Metodo dei trapezi\n3.Metodo di Cavalieri-Simpson\n");
  83.     scanf("%d", &q);    
  84.    
  85.     return q;
  86. }
  87.    
  88. double funzione(int t, double d)
  89. {
  90.      switch(t)
  91.        {
  92.           case 1:  
  93.                  return sin(d);
  94.                  break;                    
  95.                
  96.           case 2:
  97.                  return d*(exp(-d))*(cos(2*d));
  98.                  break;                    
  99.                
  100.           case 3:      
  101.                  return 1./(1+pow(d,2));
  102.                  break;                    
  103.          
  104.           case 4:      
  105.                  return pow(d, 2.5);
  106.                  break;                    
  107.              
  108.           case 5:      
  109.                  return d*sin(d);
  110.                  break;                    
  111.              
  112.           case 6:      
  113.                  return fabs(d);
  114.                  break;                    
  115.          
  116.           case 7:      
  117.                  return fabs(sin(d));
  118.                  break;                  
  119.                
  120.           case 8:      
  121.                  return exp(-(d*d));
  122.                  break;                    
  123.            
  124.           case 9:            
  125.                 if((d-M_PI)>=0)
  126.                                return 1;
  127.                 else
  128.                        return -1;
  129.                 break;                  
  130.          
  131.                  case 10:      
  132.                 if(d<=0)
  133.                                return -d;
  134.                         else
  135.                                return d*d+1;
  136.                 break;                                
  137.    }
  138. }
  139.  
  140.  
  141. void intervallo(double &a, double &b, int t)
  142. {
  143.      switch(t)
  144.        {
  145.           case 1:  
  146.             { a=0; b=2*(M_PI);}
  147.             break;                    
  148.                
  149.           case 2:
  150.             { a=0; b=2*(M_PI);}
  151.             break;
  152.                  
  153.           case 3:      
  154.             { a=-5; b=5;}
  155.             break;
  156.  
  157.           case 4:      
  158.             { a=0; b=1;}
  159.             break;
  160.  
  161.           case 5:      
  162.             { a=-2; b=2;}
  163.             break;
  164.  
  165.           case 6:      
  166.             {
  167.                printf("Inserire l'estremo a=");
  168.                scanf("%lf", &a);
  169.                
  170.                printf("Inserire l'estremo b=");
  171.                scanf("%lf", &b);          
  172.             }
  173.             break;
  174.  
  175.           case 7:      
  176.             { a=0; b=2*(M_PI);}
  177.             break;
  178.  
  179.           case 8:      
  180.             { a=-1; b=1;}
  181.             break;
  182.  
  183.           case 9:      
  184.             { a=-3; b=5;}
  185.             break;
  186.  
  187.             case 10:      
  188.             { a=-2; b=3;}
  189.             break;                                        
  190.        }    
  191. }
  192.  
  193.  
  194. void nodi(double a, double b, vettore X, int m, int t)
  195. {
  196.   double h=(b-a)/m;
  197.        
  198.   for(int i=0; i<=m; i++)
  199.       X[i]=a+(i*h);
  200. }
  201.  
  202.  
  203. void puntimedi(vettore X, int m, vettore medio)
  204. {
  205.    for(int i=0; i<m; i++)
  206.         medio[i]=(X[i+1]+X[i])/2;          
  207. }
  208.  
  209. double rettangoli(int t, int m, vettore X)
  210. {
  211.    double r=0, h=(X[m]-X[0])/m;
  212.    vettore medio;
  213.    
  214.    printf("ht=%lf", h);  
  215.    
  216.    puntimedi(X, m, medio);
  217.    
  218.    for(int i=0; i<m; i++)
  219.         {
  220.          r=r+(h*funzione(t, medio[i]));
  221.          printf("  rett=%lf ", r);
  222.         }    
  223.    
  224.    return r;
  225. }
  226.  
  227. /*Questa funzione calcola il valore dell'integrale di f(x) con il metodo dei trapezi, cioè andiamo a suddividere l'intervallo [a,b] in tanti
  228. sottointevalli equispaziati.*/
  229. double trapezi(int t, int m, vettore X)
  230. {
  231.    double trap=0, h=(X[m]-X[0])/m;
  232.    
  233.    printf("ht=%lf", h);
  234.    
  235.    for(int i=0; i<m; i++)
  236.        {
  237.           trap=trap+((h/2)*(funzione(t, X[i])+funzione(t, X[i+1])));
  238.           printf("  trap=%lf  ", trap);
  239.        }
  240.        
  241.    return trap;
  242. }
  243.  
  244.  
  245. double simpson(int t, int m, vettore X)
  246. {
  247.        double simp=0, h=(X[m]-X[0])/m;
  248.        vettore medio;
  249.        
  250.        puntimedi(X, m, medio);
  251.        
  252.        printf("hs=%lf", h);
  253.        
  254.        for(int i=0; i<m; i++)
  255.           {
  256.            simp=simp+(h/6)*(funzione(t, X[i])+4*funzione(t, medio[i])+funzione(t, X[i+1]));
  257.            printf("  simp=%lf ", simp);
  258.           }
  259.          
  260.        return simp;
  261. }


PM Quote
Avatar
nessuno (Normal User)
Guru^2


Messaggi: 5475
Iscritto: 03/01/2010

Segnala al moderatore
Postato alle 10:53
Domenica, 08/04/2012
Forse è meglio anche fare un esempio di dati da inserire, risultati attesi e risultati ottenuti ... che ne dici?


Ricorda che nessuno è obbligato a risponderti e che nessuno è perfetto ...
PM Quote
Avatar
arack95 (Member)
Pro


Messaggi: 144
Iscritto: 15/11/2010

Segnala al moderatore
Postato alle 22:01
Domenica, 08/04/2012
L'errore c'è sempre ed è normale, quelli sono metodi per calcolare gli integrali in maniera approssimativa, comunque dovresti curare il codice: nome variabili, parametri di funzioni inutilizzati(l'intero t che passi alla funzione "intervallo"), variabili non controllate(ad esempio nella funzione sceglif posso inserire un valore < 1, come 0 o -42, oppure dopo posso scegliere un numero di nodi < 0 o > 100) e altro ancora :nono:

Ultima modifica effettuata da arack95 il 08/04/2012 alle 22:04


Take life less seriously!
PM Quote
Avatar
Puffetta (Normal User)
Rookie


Messaggi: 21
Iscritto: 29/11/2009

Segnala al moderatore
Postato alle 22:24
Martedì, 10/04/2012
avete ragione entrambi... iniziamo col rispondere a nessuno: uso la funzione di Runge (caso3) e il risultato dovrebbe essere all'incirca 2,747 con 2 nodi (per tutti e tre i metodi). a me invece viene 1,3793(rettangoli), 5,1923(trapezi) e 2,6508(Simpson)

altro esempio: funzione caso2. Per un numero di nodi pari a 4 il risultato dovrebbe essere -0.122122 e invece a me viene -0,0000(rettangoli), -0.3569(trapezi) e -0.1118(Simpson)


risposta per arack95: so che il mio codice non è dei migliori, ma purtroppo non mi è stato insegnato molto bene a programmare, ho tante lacune e vorrei colmarle. dimmi altre cose che secondo te non vanno così le sistemo:-)
ho però una domanda: nella funzione intervallo mi dici che il parametro t è inutilizzato, ma se non lo passo come faccio il valore che mi interessa allo switch?

grazie a tutti, siete utilissimi:)

PM Quote
Avatar
arack95 (Member)
Pro


Messaggi: 144
Iscritto: 15/11/2010

Segnala al moderatore
Postato alle 14:04
Mercoledì, 11/04/2012
Hai ragione ho toppato io, era la funzione nodi :asd:

Per le altre cose da modificare: potresti usare le librerie standard del C++ al posto di quelle deprecate del C, potresti rendere più esplicativi i nomi delle variabili e delle funzioni...

Per gli errori, basta che aumenti il numero di nodi e i risultati saranno più precisi, del resto ce lo dice proprio la matematica :D Ho provato il 2 caso, se metti 100 nodi esce -0.121956 con i rettangoli, -0.122455 con i trapezi e -0.122122 con Simpson (valori approssimati)

Ultima modifica effettuata da arack95 il 11/04/2012 alle 14:14


Take life less seriously!
PM Quote
Avatar
Puffetta (Normal User)
Rookie


Messaggi: 21
Iscritto: 29/11/2009

Segnala al moderatore
Postato alle 15:01
Giovedì, 12/04/2012
per quanto riguarda le librerie del C purtroppo non posso toglierle visto che mi è stato imposto di usare proprio quelle, per i nomi delle variabili invece ho cercato di mettere dei nomi che si avvicinassero il più possibile alle formule teoriche studiate, in modo da non fare confuzione...

per quanto riguarda gli errori invece ho visto anche io che effettivamente sono giuste le approssimazioni quindi :)

PM Quote