C/C++ - Formule di quadratura interpolatorie

Forum - C/C++ - Formule di quadratura interpolatorie

Puffetta (Normal User)
Rookie

Messaggi: 21
Iscritto: 29/11/2009

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

//foglio di es. 4
 
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
 
typedef double vettore [100];
 
int sceglif();
int sceglimetodo();
double funzione(int t, double d);
void intervallo(double &a, double &b, int t);
void nodi(double a, double b, vettore X, int m, int t);
void puntimedi(vettore X, int m, vettore medio);
double rettangoli(int t, int m, vettore X);
double trapezi(int t, int m, vettore X);
double simpson(int t, int m, vettore X);
 
int main()
{
   int t, q, m;
   double a, b;
   vettore X;
   
   t=sceglif();
   
   intervallo(a, b, t);
   
   q=sceglimetodo();
   
   printf("Inserire il numero di nodi per partizione uniforme dell'intervallo m=");
   scanf("%d", &m);
   
   nodi(a, b, X, m, t);
   
   if(q==1)
       printf("Il valore dell'integrale con il metodo dei rettangoli e' %.16lf\n", rettangoli(t, m,  X));
   else if(q==2)      
       printf("Il valore dell'integrale con il metodo dei trapezi e' %.16lf\n", trapezi(t, m, X));
   
   else if(q==3)
       printf("Il valore dell'integrale con il metodo di C-S e' %.16lf\n", simpson(t, m, X));
    
   system("PAUSE");
   return 0;  
}
 
 
int sceglif()
{
  int t=0;
  
  printf("\n Scegli la funzione da approssimare con il polinomio di Lagrange:\n\n");
  
  printf(" 1) integrale di sen(x) per x in [0;2pi]\n\n");
  printf(" 2) integrale di x*e^(-x)*cos(2x) per x in [0;2pi]\n\n");
  printf(" 3) integrale di 1/(1+x^2) per x in [-5;5]\n\n");
  printf(" 4) integrale di x^5/2 per x in [0;1]\n\n");
  printf(" 5) integrale di x*sin(x) per x in [-2;2]\n\n");
  printf(" 6) integrale di mod(x) per x in [-1;0] oppure in [-1;1]\n\n");
  printf(" 7) integrale di mod(sin(x)) per x in [0;2pi]\n\n");
  printf(" 8) integrale di e(-x^2) per x in [-1;1]\n\n");
  printf(" 9) integrale di sgn(x-pi) per x in [-3;5]\n\n");
  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");
  
  printf("Quale funzione devo usare? t=");
  scanf("%d", &t);
  
  while(t>10)
        {
       printf("\n Il valore inserito non è corretto!!! Inserire un valore corretto di t t=");
       scanf("%d", &t);
        }
 
  return t;
}    
 
int sceglimetodo()
{
    int q;
    
    printf("Quale metodo vuoi utilizzare?\n1. Metodo dei rettangoli\n2.Metodo dei trapezi\n3.Metodo di Cavalieri-Simpson\n");
    scanf("%d", &q);    
    
    return q;
}
    
double funzione(int t, double d)
{
     switch(t)
       {
          case 1:   
                 return sin(d);
                 break;                    
                
          case 2:
                 return d*(exp(-d))*(cos(2*d));
                 break;                    
               
          case 3:       
                 return 1./(1+pow(d,2));
                 break;                    
          
          case 4:       
                 return pow(d, 2.5);
                 break;                    
              
          case 5:       
                 return d*sin(d);
                 break;                    
              
          case 6:       
                 return fabs(d);
                 break;                     
          
          case 7:       
                 return fabs(sin(d));
                 break;                  
               
          case 8:       
                 return exp(-(d*d));
                 break;                    
           
          case 9:             
                if((d-M_PI)>=0) 
                               return 1;
                else
                       return -1;
                break;                   
          
                 case 10:       
                if(d<=0) 
                               return -d;
                        else
                               return d*d+1;
                break;                                
   }
}
 
 
void intervallo(double &a, double &b, int t)
{
     switch(t)
       {
          case 1:   
            { a=0; b=2*(M_PI);}
            break;                    
                
          case 2:
            { a=0; b=2*(M_PI);}
            break;
                 
          case 3:       
            { a=-5; b=5;}
            break;
  
          case 4:       
            { a=0; b=1;}
            break;
 
          case 5:       
            { a=-2; b=2;}
            break;
 
          case 6:       
            { 
               printf("Inserire l'estremo a=");
               scanf("%lf", &a);
               
               printf("Inserire l'estremo b=");
               scanf("%lf", &b);           
            }
            break;
 
          case 7:       
            { a=0; b=2*(M_PI);}
            break;
 
          case 8:       
            { a=-1; b=1;}
            break;
 
          case 9:       
            { a=-3; b=5;}
            break; 
 
            case 10:       
            { a=-2; b=3;}
            break;                                         
       }    
}
 
 
void nodi(double a, double b, vettore X, int m, int t)
{
  double h=(b-a)/m;
        
  for(int i=0; i<=m; i++)
      X[i]=a+(i*h);
}
 
 
void puntimedi(vettore X, int m, vettore medio)
{
   for(int i=0; i<m; i++)
        medio[i]=(X[i+1]+X[i])/2;          
}
 
double rettangoli(int t, int m, vettore X)
{
   double r=0, h=(X[m]-X[0])/m;
   vettore medio;
   
   printf("ht=%lf", h);  
   
   puntimedi(X, m, medio);
   
   for(int i=0; i<m; i++)
        {
         r=r+(h*funzione(t, medio[i]));
         printf("  rett=%lf ", r);
        }    
   
   return r;
}
 
/*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
sottointevalli equispaziati.*/
double trapezi(int t, int m, vettore X)
{
   double trap=0, h=(X[m]-X[0])/m;
   
   printf("ht=%lf", h);
   
   for(int i=0; i<m; i++)
       {
          trap=trap+((h/2)*(funzione(t, X[i])+funzione(t, X[i+1])));
          printf("  trap=%lf  ", trap);
       }
       
   return trap;
}
 
 
double simpson(int t, int m, vettore X)
{
       double simp=0, h=(X[m]-X[0])/m;
       vettore medio;
       
       puntimedi(X, m, medio);
       
       printf("hs=%lf", h);
       
       for(int i=0; i<m; i++)
          {
           simp=simp+(h/6)*(funzione(t, X[i])+4*funzione(t, medio[i])+funzione(t, X[i+1]));
           printf("  simp=%lf ", simp);
          }
          
       return simp;
}

nessuno (Normal User)
Guru^2

Messaggi: 6402
Iscritto: 03/01/2010

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 ...
---
Il grande studioso italiano Bruno de Finetti ( uno dei padri fondatori del moderno Calcolo delle probabilità ) chiamava il gioco del Lotto Tassa sulla stupidità.

arack95 (Member)
Pro

Messaggi: 144
Iscritto: 15/11/2010

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

Puffetta (Normal User)
Rookie

Messaggi: 21
Iscritto: 29/11/2009

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

arack95 (Member)
Pro

Messaggi: 144
Iscritto: 15/11/2010

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

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

Puffetta (Normal User)
Rookie

Messaggi: 21
Iscritto: 29/11/2009

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