Lezione 5

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 separate.

Formule magiche

Finora abbiamo usato il makefile nel modo piu` esplicito possibile, indicando esattamente cosa bisognava fare per ogni mattone da costruire. L'utilita` del makefile si comprende di piu` se si tiene conto che questo puo` definire delle regole generiche ed inoltre puo` utilizzare delle variabili predefinite molto utili.

Per esempio la seguente scrittura puo` risultare molto utile:

%.o : %.cxx
       g++ -Wall -c $<

Si noti il significato speciale del carattere %, che indica "un qualunque target e la variabile $< che serve ad indicare la prima dipendenza."

Altre variabili molto utili sono $@, che indica il nome del target, e $^, che indica la lista completa dei prerequisiti. Ad esempio, ritornando al primo esercizio, per compilare un programma di ordinamento che usava il metodo del quicksort, si sarebbe potuto scrivere:

quicksort : quicksort.o ordinamenti.o scambia.o
	g++ -o $@ $^

Le funzioni integrande

• funzioni.cxx
• funzioni.h

Se si necessita di ulteriori funzioni, sara` sufficiente aggiungerle a questo file.

Il programma principale

Per variare un po' rispetto ai programmi passati, proviamo ad inserire nel programma principale un menu` 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;
    }
  }

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

Metodi di integrazione

• integrando una o piu` funzioni di cui l'integrale e' noto;
• verificare come la precisione del metodo dipende dal passo di integrazione (comparare i vari metodi a parita` di passo e lo stesso metodo variando il passo)
• negli esempi che seguono, useremo float invece di double, questo permette di incorrere piu` rapidamente in problemi di precisione numeriche. Questo ha una funzione puramente didattica.

Metodo dei trapezoidi

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;
}

Dopo che siete riusciti ad ottenere un programma compilabile, provate a vedere se fornisce il risultato corretto su alcune funzioni di prova (ad esempio sin(x) o la radice quadrata). e' istruttivo prendere un intervallo in cui la funzione varia molto e, partendo con pochi intervalli (N=2,4...), vedere quando si ottiene un risultato a vostro avviso accettabile.

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 convessita' della funzione integranda. La sua applicazione richiede che N sia pari e la formula da usare e':


La complicazione rispetto al programma precedente e' che nella lettura da terminale di N, 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 sara' 0 se N e' pari o 1 se N e' 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);
realizzanto all'interno della funzione un controllo sulla parita` del numero di intervalli, eventualmente riconducendolo al numero pari immediatamente superiore.

Si verifichi su alcune funzioni di prova che la formula di Simpson e' infatti piu' precisa di quella dei trapezoidi.

Periodo del pendolo

Uno degli esercizi della lezione 2 chiedeva di trovare la dipendenza del periodo di un pendolo reale dall'ampiezza di oscillazione integrando l'equazione differenziale che descrive il sistema.

Un altro approccio per risolvere il problema e' di sfruttare la conservazione dell'energia che, nel caso del pendolo, permette di collegare la velocita` in un punto qualunque della traiettoria con la l'ampiezza massima di oscillazione:

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

Si puo` 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 singolarita` negli estremi di integrazione, quindi non si puo` integrare fino agli estremi, ma bisogna un attimo prima;
• l'ampiezza θ_max compare sia negli estremi di integrazione che nella funzione integranda, sara` 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:

Formule magiche 2

Parte opzionale

Integrazione a precisione prefissata

Per migliorare la precisione dell'integrazione, una possibile soluzione e' 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 e' necessario ricalcolare tutto l'integrale, ma basta aggiungere alla sommatoria dei valori della funzione la valutazione nei nuovi punti aggiunti a meta' tra i punti precedentemente calcolati.

Questa possibilita' 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;

}

IMPORTANTE: 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, e' opportuno far stampare il numero di passi necessari per raggungere la precisione voluta.

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:

• campomagnetico.cxx
• campomagnetico.h

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:

Relazione