Lezione 3

Introduzione

Scopo di questa sessione di laboratorio e' di implementare diversi metodi di integrazione numerica su di una variabile, in modo da poterne verificare la precisione.

Verranno poi forniti degli esempi di applicazione di questi metodi alla risoluzione di problemi fisici.

Anche questa volta utilizzeremo la solita struttura di avere una funzione principale che permette di selezionare alcune integrande predefinite, chiede gli estremi di integrazione ed effettua il calcolo.

Le funzioni integrande saranno definite in un file a parte e dovranno essere funzioni che ricevono come argomento un double e restiruiscono un double.

I diversi metodi di integrazione, saranno invece create in funzioni che converrà tenere separate, per poterle riustilizzare in seguito.

Le funzioni integrande

Le funzioni integrande che realizzeremo sono della forma
double funzione(double); La maggior parte delle funzioni delle librerie matematiche, come sin, cos, sono già definite con questa convenzione. Alcune funzioni che useremo sono definite in

funzioni.cxx
funzioni.h

Se si necessita di ulteriori funzioni, sarà sufficiente aggiungerle a questo file.

In funzioni.h si noti l'utilizzo delle direttive del preprocessore per evitare multiple inclusioni di un header file: questa è una buoan abitudine per evitare catene di file che si includono reciprocamente.

La struttura è del tipo:

#ifndef VARIABILE_CHE_IDENTIFICA_IL_FILE
#define VARIABILE_CHE_IDENTIFICA_IL_FILE
.
.
.
#endif

dove la prima riga controlla se è già stata definita la variabile indicata (usare un nome che identifica in maniera ragionevolmente univoca il file da usare). Se questa variabile è già stata definita, significa che il file è già stato letto d il resto del file, fino alla direttiva #endif viene saltato.

Se invece la variabile non è stata ancora definita, la parte di header fino all'#endif viene letta, ma la prima operazione che essa fa è di definire la variabile, in modo che si sia certi che essa verrà trovata definita al momento della prossima inclusione.

Il programma principale

Per variare un po' rispetto ai programmi passati, proviamo ad inserire nel programma principale un menù per selezionare le funzioni da integrare.

Il frammento di codice che implementa questa funzione potrebbe essere:


  int id;
  double (*fun)(double)=0; // definisco fun come un puntatore a funzione
                           // il valore iniziale della funzione e` nullo
  cerr << "Seleziona la funziona da integrare" << endl;
  cerr << " 1 - seno" << endl;
  cerr << " 2 - coseno" << endl;
  cerr << " 3 - radice quadrata" << endl;
  cerr << " 4 - quadrato" << endl;
  while ( fun==0 ) {       // ripeto il ciclo fino a quando non riesco ad 
                           // inizializzare correttamente fun
    cin >> id;
    switch (id) {
      case 1:
        fun=sin;           // sin e` gia` definito in math.h
        break;
      case 2:
        fun=cos;           // cos e` gia` definito in math.h
        break;
      case 3:
        fun=radice;        // radice e` definito in funzioni.h
        break;
      case 4:
        fun=quadrato;      // quadrato e` definito in funzioni.h
        break;
      default:
        cerr << "La funzione selezionata: id=" << id << " non esiste!" << endl;
        cerr << "Seleziona la funzione da integrare: ";
        break;
    }
  }

Una struttura simile si potrà utilizzare anceh per scrivere un menù che permetta di scegliere uno dei diversi metodi di integrazione illustrati nella prossima sezione.

Piu` avanti nel programma, ci sarà da chiamare il metodo che effettua l'integrazione, passando come argomento la funzione integranda selezionata. Ad esempio,

cerr << "Estremo inferiore: ";
cin >> xmin;
cerr << "Estremo superiore: ";
cin >> xmax;
cerr << "Numero di intervalli: ";
cin >> npassi;
double integrale = trapezoidi(fun,xmin,xmax,npassi)

Metodi di integrazione

Il problema che vogliamo affrontare inizialmente è quello di calcolare un integrale del tipo:

suddividendo l'intervallo iniziale [a,b] in N intervalli di lunghezza

ed applicando su ciascuno di essi un approssimazione della funzione. Adesso provare ad implementare i diversi metodi di integrazione spiegati a lezione. Per tutti, bisogna controllare che l'implementazione sia corretta:

integrando una o più funzioni di cui l'integrale è noto;
verificare come la precisione del metodo dipende dal passo di integrazione (comparare i vari metodi a parità di passo e lo stesso metodo variando il passo)

Metodo del midpoint

In questo metodo l'approssimazione dell'integrale è data dalla formula:

che si può dimostrare avere un errore sul calcolo dell'integrale O(h²).

Per verificarlo, fare un grafico dell'errore numerico in funzione di h² e verificareche i punti stanno su di una retta.

Nota bene: questo è l'unico metodo che vedremo che non necessita di valutare la funzione negli estremi dell'intervallo [a,b].

Metodo dei trapezoidi

Il metodo dei trapezoidi è dato dalla formula:
Formula dei trapezoidi

Quindi un approccio per la funzione potrebbe essere:


double trapezoidi(double (*funzione)(double), float xmin, float xmax, int npassi) {
  float integrale=0;
  float h=(xmax-xmin)/npassi;
  integrale = // inizializzazione usando i valori della funzione agli estremi
  for ( int i=1; i<npassi; i++ ) {
    float x = // calcolare il valore di x
    integrale += h*(*funzione)(x);
  }
  return integrale;
}

Anche questo metodo ha un errore proporzionale ad h².

Integrazione alla Simpson

Mantenendo lo stesso numero di valutazioni della funzione sull'intervallo, la formula di Simpson migliora la precisione dell'integrale tenendo in conto anche la convessità della funzione integranda e fornendo un errore di troncamento pari a h⁴. La sua applicazione richiede che N sia pari e la formula da usare è:

La complicazione rispetto al programma precedente è che nella lettura da terminale di npassi, bisogna controllare che questo sia pari. Inoltre all'interno del ciclo for bisogna usare pesi differenti per i termini pari (2/3) e quelli dispari (4/3) della sommatoria.

A questo proposito, il C++ fornisce un utile operatore: l'operatore % "resto della divisione per". Ad esempio il valore dell'espressione N%2 sarà 0 se N è pari o 1 se N è dispari.

Si costruisca una funzione che implementi il metodo di Simpson e soddisfi il prototipo
double simpson(double (*funzione)(double), float xmin, float xmax, int npassi);
realizzando all'interno della funzione un controllo sulla parità del numero di intervalli, eventualmente riconducendolo al numero pari immediatamente superiore.

Si verifichi su alcune funzioni di prova che la formula di Simpson è infatti più precisa di quella dei trapezoidi o del midpoint.

Uso di TGraph

Per verificare la dipendenza funzionale degli errori dei diversi metodi in funzione del passo, può essere utile stampare a schermo le coppie di valori:
h² <integrale calcolato - valore vero >
e eventualmente mandarli su di un file usando la ridirezione dello standard output.

All'interno di ROOT, si possono rileggere i dati usando un costruttore molto utile degli oggetti TGraph: si puo` costruire grafico contenente le coppie semplicemente usando il costruttore:

TGraph grafico("nomefile")

Una volta costruito il grafico, si può verificarne la dipendenza funzionale facendone un fit polinomiale:

grafico.Fit("pol2");

Integrazione a precisione prefissata

Per migliorare la precisione dell'integrazione, una possibile soluzione è aumentare il numero di punti di integrazione. La formula dei trapezoidi ha una caratteristica attraente per fare questa operazione: se si raddoppia il numero dei sottointervalli, non è necessario ricalcolare tutto l'integrale, ma basta aggiungere alla sommatoria dei valori della funzione la valutazione nei nuovi punti aggiunti a metà tra i punti precedentemente calcolati.

Questa possibilità offre lo spunto per un programma che abbia in input solo gli estremi dell'intervallo e poi proceda a valutare l'integrale raddoppiando ogni volta il numero di intervalli fin tanto che l'errore relativo sull'integrale non raggiunge la precisione desiderata. Per valutare l'errore, possiamo considerare la differenza tra due valutazioni successive dell'integrale.

La funzione portebbe avere un prototipo del tipo:
double iterativo(double (*funzione)(double), float xmin, float xmax, float precisione);

Siccome il numero di passi non e' noto a priori, concettualmente il processo di raddoppio degli intervalli puo' stare bene all'interno di un ciclo while:

while ( fabs( (newint-oldint)/oldint ) > precisione ) {
  oldint=newint;
  npassi*=2;
  h/=2;
  for ( /* costruite il for in modo che si cicli solo sui nuovi punti */ ) {
    x=...
    sum+=(*funzione)(x);
}
newint = h*sum;
}

La variabile sum contiene la somma delle valutazioni delle funzioni (la parte tra parentesi nella formula dei trapezoidi) e si e' preferito separarla dalla variabile contentente l'integrale perche' le due sono legate dal fattore h che varia variando il numero di intervalli.

Nota Bene: in questo programma bisogna fare molta attenzione all'inizializzazione corretta di tutte le variabili!

Al termine del programma, oltre a farsi stampare il valore dell'integrale, è opportuno far stampare il numero di passi necessari per raggungere la precisione voluta.

Midpoint iterativa

Le stesse considerazione fatte per il metodo dei trapezoidi, si possono applicare a quello del midpoint. Solo che in questo caso bisogna triplicare i punti tra un'iterazione e l'altra.

È opportuno implementare anche tale metodo perché, non richiedendo di valutare la funzione sugli estremi dell'intervallo, risulta utilizzabile anche in caso di singolarità sugli estremi di integrazione, come nel caso che segue.

Approfondimenti

Prima di proseguire con gli esercizi, consultate i consigli su come trattare integrali impropri.

Periodo del pendolo

Un problema tipico in meccanica è trovare la dipendenza del periodo di un sistema fisico reale dall'ampiezza di oscillazione.

Un possibile approccio per risolvere il problema è di sfruttare la conservazione dell'energia che, nel caso di moto in un potenziale, permette di collegare la velocità in un punto qualunque della traiettoria, basandosi sul calcolo dell'energia cinetica come differenza tra l'energia totale e quella potenziale.

Come esempio, consideriamo il caso di un pendolo reale, ovvero una particella che si muova in un potenziale -mlgcosθ. La conservazione dell'energia si può scrivere come:

Separando il t e θ ed integrando, si ottiene la seguente relazione per il periodo di un'oscillazione:

Si pò` quindi determinare il periodo del pendolo integrando l'equazione precedente.

Fare un grafico del periodo del pendolo in funzione dell'ampiezza di oscillazione θ_max, per valori di θ_max tra 0.01 ed 1.5 rad. Come valore dei parametri, utilizzate g=9.8 m/s² e l=30 cm.

Ci sono un paio di accortezze da avere:

la funzione integranda presenta una singolarità negli estremi di integrazion;
l'ampiezza θ_max compare sia negli estremi di integrazione che nella funzione integranda, sarà quindi utile utilizzare una variabile condivisa tra la funzione principale e quella che chiama il metodo di integrazione.

Alla fine dovreste ottenere un grafico simile al seguente:

Parte opzionale

Calcolo del campo magnetico prodotto da bobine di Helmotz

Adesso avete in mano tutti gli strumenti per risolvere un altro problema numerico che troverete in un'esperienza del laboratorio del secondo semestre.

Siccome non aveta ancora studiato la legge di Biot-Savart, potete per il momento lasciare in sospeso questo esercizio, ma ripescarlo quando studierete tale esperienza.

Innanzitutto definiamo il problema di calcolare il campo magnetico generato da un circuito in cui passa una corrente. Questo e' dato dalla legge di Biot e Savart:

Nel caso specifico di una spira circolare di raggio R, con il centro nell'origine e orientata nel piano x-y:

le quantita` nella funzione integranda diventano:

I seguenti file:

contengono uno schema per realizzare delle funzioni che generino il campo magnetico prodotto da una singola spira in un punto generico dello spazio: dovete solo inserire correttamente il vostro algoritmo di integrazione.

L'apparato che userete in laboratorio per la determinazione del rapporto e/m dell'elettrone, consiste di una coppia di bobine, composta ciascuna di 120 spire, di raggio 15.5 cm e poste ad una distanza di 15.5 cm.

Costruite una funzione che calcoli il campo magnetico generato nel piano intermedio tra le bobine per una corrente di 1 A.

Fare un grafico del rapporto tra la componente B_z del campo ad un certo raggio dal centro delle bobine ed il suo valore al centro del sistema. La forma del grafico dovrebbe ricalcare la seguente:

Integrazione alla Romberg

Come indicato nella lezione, si può ridurre l'errore di troncamento effettuando un'interpolazione polinomiale attraverso una sequenza di integrali calcolati con diversa granularità.

Per esempio con il metodo dei trapezoidi:

calcolando gli integrali per N, 2N e 4N punti
usando i valori ottenuti per costruire un polinomio di secondo grado in h², passante per i punti con ascissa rispettivamente (b-a)²/N², (b-a)²/4N², (b-a)²/16N²,
usare come stime dell'integrale l'estrapolazione di questo polinomio per h²=0.

Implementare il metodo suddetto e verificare che è di ordine h⁶.

Laboratorio di Calcolo 2