Lezione 6 : ricerca zeri di funzione

In questa sesta lezione affronteremo il problema della ricerca di zeri di una funzione. Per fare questo realizzeremo per prima cosa due classi astratte per rappresentare rispettivamente una generica funzione di una variabile $f(x)$ e un generico algoritmo per la ricerca di zeri. Affronteremo poi casi concreti (ricerca degli zeri per una parabola usando il metodo della bisezione) sfruttanto ereditarietà e polimorfismo.

ESERCIZIO 6.0 - Metodi virtual :

Prima di incominciare con la lezione vera e propria proviamo a riflettere sul significato del polimorfismo attraverso l'uso del qualificatore virtual. Consideriamo le classi Particella e la sue derivata Elettrone costruite nella lezione 5. Implementare ed eseguire il seguente programma:

#include "Particella.h"

int main() {

  Particella *p = new Particella(1.,2.);
  Elettrone  *e = new Elettrone();
  Particella *b = new Elettrone(); // Puntatore a Particella che punta ad un oggetto Elettrone

  p->Print(); // Metodo Print di Particella
  e->Print(); // Metodo Print di Elettrone
  b->Print(); // Quale Print ???

  delete p;
  delete e;
  delete b;

  return 0;

}

Mentre le prime due chiamate p->Print(); e e->Print(); sono immediate da capire, il caso più interessante è quello di b : b è tecnicamente di tipo puntatore a Particella ma punta ad un oggetto di tipo Elettrone. Questa sitauzione è permessa in C++ e costituisce la base del polimorfismo.

• Se utilizziamo le classi scritte nella lezione 5 senza modifiche, si dovrebbe constatare nell'istruzione b->Print() viene invocato il metodo Print() della classe madre (Particella), anche quando gli oggetti riferiti dai puntatori sono classi figlie (Elettrone in questo caso).
• Aggiungere ora negli header file delle classi Particella ed Elettrone il qualificatore virtual davanti alla dichiarazione del metodo Print(). Ricompilando e rigirando il programma si dovrebbe adesso vedere che per ciascun oggetto viene invocato il metodo corrispondente all'oggetto al quale il puntatore punta: nel caso specifico in b->Print() viene invocato il metodo Print() di Elettrone in quanto b punta ad un oggetto di tipo Elettrone.
• Una classe che implementa o eredita un metodo virtuale si dice polimorfica: l'istruzione b->Print() darà esiti diversi a seconda della classe a cui b punta nonstante tecnicamente b sia di tipo puntatore a Particella. Questa proprietà del linguaggio sarà molto utile per sviluppare il codice di ricerca degli zeri di una funzione.

Metodi virtuali e classi astratte

La parola chiave virtual informa il compilatore che il metodo indicato potrà venire sovrascritto dalle classi figlie. Se essa viene fornita, durante l'esecuzione del programma, al momento di chiamare il metodo a partire da un puntatore alla classe madre, il programma valuterà se l'oggetto indirizzato è del tipo classe madre o una delle figlie. In quest'ultimo caso, invocherà il metodo appropriato della classe figlia reale. Un metodo virtuale può anche essere posto a 0 (vedi prossimo esercizio), in tal caso è obbligatorio per le classi figlie implementarlo. In questo caso il metodo si dice virtuale puro e una classe che implementa almeno un metodo virtuale puro si dice classe astratta.

Nota sull'organizzazione del codice

In linea di principio si potrebbe pensare di costruire una coppia di files .h e .cpp per ogni classe che scriviamo. Questo potrebbe portare ad una proliferazione incontrollata del numero di files del nostro progetto. Per limitare questo effetto possiamo pensare di scrivere nello stesso file un set di classi affini tra loro: in questo caso specifico per esempio Particella ed Elettrone possono stare ragionevolemnte nello stesso file.

ESERCIZIO 6.1 - Classe astratta FunzioneBase :

In questo e nei prossimi esercizi avremo a che fare con diverse funzioni di una sola variabile e dovremo effettuare delle operazioni generiche su queste funzioni (come trovarne gli zeri o calcolarne l'integrale in un certo intervallo). Dovremo pertanto strutturare il codice in modo che ogni algoritmo possa lavorare allo stesso modo su funzioni che possono essere diverse. Per fare questo è utile definire una classe astratta che implementi le proprietà generali della classe attraverso metodi virtuali puri e poi lasciare alle classi derivate il compito di implementare concretamente tutti questi metodi ( ed eventualmente quelli aggiuntivi necessari). Notate in questo caso che è buona prassi definire il distruttore di una classe madre astratta come virtual in modo che il distruttore della classe figlia venga invocato correttamente quando un oggetto di tipo classe figlia viene distrutto. Per capire meglio la struttura del codice procediamo per passi. In questo esercizio scriviamo un set di classi per modellizzare una generica funzione di una variabile:

• definire la classe astratta FunzioneBase che implementi un metodo double Eval(double x) che modellizza l'effetto di applicare $f$ in $x$.
• implementare una classe derivata Parabola che descriva una funzione del tipo $f(x)=ax^{2}+bx+c$ (chiaramente questa classe dovrà avere dei data membri aggiuntivi per i parameteri a, b, c insieme a metodi dedicati per l'accesso in lettura e scrittura a questi parametri).
• verificare il funzionamento della classe Parabola costruendo un programmino che, dati i parametri di una parabola, ne stampi il valore della $y$ nel vertice $x_{v}=-b/2a$. A titolo di esempio utilizzare la parabola $f(x)= 3x^{2}+5x-2$.

La classe FunzioneBase

class FunzioneBase {

public:

  virtual double Eval(double x) const =0;
  virtual ~FunzioneBase() {;};

};

La classe Parabola

class Parabola : public FunzioneBase {

public:

  Parabola()  { m_a = 0.; m_b = 0.; m_c = 1.; }  
  Parabola(double a, double b, double c)  { m_a = a; m_b = b; m_c = c; }
  ~Parabola() {;} ;
  virtual double Eval(double x) const override {return m_a*x*x+m_b*x+m_c;}
  void SetA(double a) { m_a = a; }
  double GetA() const {return m_a;} 
  void SetB(double b) { m_b = b; }
  double GetB() const {return m_b;} 
  void SetC(double c) { m_c = c; }
  double GetC() const {return m_c;} 

  double GetVertex() const { return -m_b / (2*m_a) ;} ;

private:

  double m_a, m_b, m_c;

};

Notate la keyword override che ci avvisa con un errore se cerchiamo per errore di definire nella classe figlia il metodo virtuale della classe madre con parametri di input diversi.

Nota sulle classi astratte

Quando dichiariamo nullo un metodo virtual, la classe difetta dell'implementazione di tale metodo, e quindi non si possono creare oggetti di quella classe. Solo le sue classi derivate implementeranno il metodo, e possono essere implementate. Siccome non si possono costruire oggetti di tale classe, non è necessario definire dei costruttori. Se una classe astratta implementa solo metodi virtuali nulli, come in questo caso, non è neanche necessario realizzare un file di implementazione, né creare un file oggetto, dato che tutta l'informazione è contenuta nell'header.

ESERCIZIO 6.2 - Metodo della bisezione (da consegnare):

Passiamo ora alla codifica della parte relativa agli algoritmi di ricerca degli zeri di una funzione: proviamo a scrivere un programma che calcoli gli zeri della funzione $f(x) = 3x^{2}+5x-2$. Il programma dovrà leggere da riga di comando gli estremi dell'intervallo in cui cercare lo zero e la precisione richiesta. Per calcolare gli zeri, implementeremo una classe astratta Solutore ed una classe concreta che realizza il metodo della bisezione visto a lezione. Nulla vieta di implementare anche il metodo delle secanti usando lo stesso schema. Si richiede obbligatoriamente che il programma stampi l'ascissa dello zero con un numero di cifre significative corrispondente alla precisione immessa.

La funzione segno

Ci sono vari modi di implementare una funzione segno da utilizzare nella codifica dell'algoritmo di bisezione.

1. funzione libera

   int sign(double x){
   
    if(x<0)
      return -1;
    else 
      return 1;
   }
   

   oppure con l implementazione equivalente ma molto più compatta

   double sign(double x){return (x==0.?0.:(x>0?1.:-1)); };
   

2. classe che erediti da FunzioneBase

   class Segno : public FunzioneBase {
     public :
       double Eval( double x) const { return (x==0.?0.:(x>0?1.:-1)); } ;
   }
   

   Usando questa funzione la condizione richiesta dall'algoritmo di bisezione sarà ( implementando segno come funzione libera )

if((sign(f(a))*sign(f(c)))<0)
  ......

else if((sign(f(b))*sign(f(c)))<0)
  ......

else
  ......

Classe astratta Solutore

La classe astratta Solutore potrebbe avere un metodo virtuale puro CercaZeri corrispondente alla chiamate dell'algoritmo che cercherà di determinare gli zeri di una generica FunzioneBase, passata come puntatore o come referenza: nell'esempio sono presentati entrambi i casi, se siete indecisi optate per l'uso di una referenza. Inoltre possiamo definire dei metodi per configurare la precisione richiesta: tale precisione può essere definita nel costruttore, tramite un metodo dedicato o direttamente nella chiamata al metodo CercaZeri. Stesso discorso vale per il numero massimo di iterazione dell'algoritmo. Come nel caso della FunzioneBase notate l'implementazione del distruttore come metodo virtuale e l'utilzzo della keyword override.

class Solutore {

 public:

// costruttori : potete aggiungere costruttori che ritenete opportuni e/o necessari

  Solutore() ;
  Solutore( double prec ) ;


  virtual ~Solutore() {;};

  // a scopo illustrativo trovate le due versioni una con puntatore FunzioneBase *
  // e una con reference FunzioneBase&. Chiaramente ne basta una !

  virtual double CercaZeriPointer(double xmin ,
                                  double xmax ,
                                  const FunzioneBase* f,
                                  double prec = 0.001,
                                  unsigned int nmax = 100) = 0;

  virtual double CercaZeriReference(double xmin ,
                                    double xmax ,
                                    const FunzioneBase& f,
                                    double prec = 0.001,
                                    unsigned int nmax = 100) = 0;

  void SetPrecisione(double epsilon) { m_prec = epsilon; }
  double GetPrecisione() { return m_prec;}

  void SetNMaxiterations(unsigned int n ) { m_nmax = n ; }
  unsigned int GetNMaxiterations () { return m_nmax ; } ;

  unsigned int GetNiterations () { return m_niterations ; } ;

 protected:

  double m_a, m_b;              // estremi intervallo di ricerca
  double m_prec;                // precisione richiesta
  unsigned int m_nmax;          // massimo numero di iterazioni permesse 
  unsigned int m_niterations;   // numero di iterazioni effettuate

};

La classe concreta Bisezione

L'implementazione dell'algoritmo di bisezione dovrà necessariamente avvenire costruendo una classe dedicata Bisezione che erediti da Solutore e implementi una versione concreta del metodo CercaZeri.

class Bisezione : public Solutore {

 public:

  Bisezione() {} ;
  Bisezione(double prec );
  virtual ~Bisezione();

  // a scopo illustrativo trovate le due versioni una con puntatore FunzioneBase * 
  // e una con reference FunzioneBase&. Ne basta una sola, meglio con la reference                                              
  // e comunque coerente con quanto dichiarato nella classe madre                                                             

  virtual double CercaZeriPointer(double xmin,
                                  double xmax,
                                  const FunzioneBase* f,
                                  double prec = 0.001 ,
                                  unsigned int nmax = 100) override ;

  virtual double CercaZeriReference(double xmin, double xmax,
                                const FunzioneBase& f,
                                double prec = 0.001,
                                unsigned int nmax = 100) override ;

};

Il metodo della bisezione

Data una funzione $f(x)$, si dice zero, o radice, di $f$ un elemento $x_{0}$ del suo dominio tale che $f(x_{0}) = 0$.

Per trovare numericamente gli zeri della funzione $f(x)$ possiamo utilizzare il metodo della bisezione: è l'algoritmo più semplice e consiste in una procedura iterativa che, ad ogni ciclo, dimezza l'intervallo in cui si trova lo zero. Dal teorema di Bolzano (o degli zeri) sappiamo che data una funzione $f(x)$ continua sull'intervallo chiuso $[a,b]$ a valori in $\mathbb{R}$ tale che $f(a)*f(b) < 0$, allora esiste un punto $x_{0}$ all'interno di $[a,b]$ tale che $f(x_{0}) = 0$.

Definiamo intervallo di incertezza di $f$ un intervallo $[a,b]$ che soddisfa il teorema di Bolzano. L'idea di base dell'algoritmo è che se esiste un intervalo di incertezza $[a,b]$ per una funzione $f(x)$, allora ne esiste uno più piccolo (esattamente la metà). L'algoritmo deve avere in input l'intervallo di incertezza di partenza $[a,b]$ ed una precisione (o tolleranza) $\epsilon$ con cui si vuole trovare lo zero di $f(x)$ tale che $|b - a| < \epsilon$.

Si parte quindi dividendo in due l'intervallo trovando il punto medio $c = a + 0.5*(b-a)$ per cui avremo due intervalli $[a,c)$ e $(c,b]$. Ora se $f(c) = 0$ siamo fortunati e abbiamo trovato lo zero. Altrimenti si deve valutare $f(a)*f(c)$ e $f(c)*f(b)$ e ripetere la procedura sull'intervallo in cui $f(x)$ cambia di segno. La procedura va ripetuta finchè la larghezza dell'intervallo finale non è minore di $\epsilon$.

Ci sono alcuni caveat da tenere in considerazione : * nell'implementazione delle condizioni di ricerca dall'intervallo di incertezza occorre prestare attenzione alle operazioni tra floating point soprattutto in prossimità della radice. Ad esempio le espressioni $f(a)*f(c)$ e $f(c)*f(b)$ hanno una buona probabilità di essere approssimata a zero dal momento che entrambi gli argomenti convergono a una radice di $f$. Per evitare questa eventualità è meglio valutare il prodotto dei segni $sign(f(a))*sign(f(c))$ e $sign(f(c))*sign(f(b))$. * un controllo utile da implementare è quello di contare il numero di iterazioni dell'algoritmo e stampare un avviso nel caso queste superino una soglia $nmax$ che può essere impostata dall'utente. In tal modo ci si accorge se ci sono possibili problemi nel ciclo (errori o richiesta di precisione troppo alta che innesca loop infiniti). * Se l'intervallo contiene più di una radice il metodo della bisezione ne troverà solo una !

Precisione sulle cifre significative

Poiché la precisione richiesta all'algoritmo è passata al programma runtime, abbiamo bisogno di determinare runtime quante cifre significative stampare nel nostro risultato. È facile rendersi conto che il numero di cifre significative è dato da

int cifre_significative = -log10(precision);

Per cui per impostare il numero di cifre significative nella scrittura a video il codice sarà

cout <<fixed;
cout <<setprecision(cifre_significative) <<zero <<endl;

Valori di default nelle funzioni

Avrete probabilmente notato che nella dichiarazione di CercaZeri abbiamo assegnato dei valori di default a prec e nmax.

virtual double CercaZeriReference(double xmin, double xmax,
                              const FunzioneBase& f,
                              double prec = 0.001,
                              unsigned int nmax = 100) override ;

Questo significa che le ulttime due variabili possono essere omesse nella chiamata nel main

bisettore.CercaZeriReference( xmin, xmax, f );

e il metodo utilizzerà i valori si 0.001 e 100 per prec e nmax.

ESERCIZIO 6.3 - Equazioni non risolubili analiticamente (da consegnare):

In problemi di meccanica quantistica che verranno studiati nel prossimo anno, ci si può imbattere in equazioni del tipo: $$x = tan(x)$$ È facile rendersi conto che tale equazione ha una soluzione in ciascuno degli intervalli $(n\pi,n\pi+\pi/2), n=1,2,3...$ Calcolare con una precisione di almeno $10^{-6}$ i valori delle soluzioni per $n$ da 1 a 20.

Suggerimento

Riscrivere l'equazione come $sin(x)-xcos(x)=0$.

ESERCIZIO 6.4 - Ricerca di zeri di una funzione senza uso del polimorfismo (facoltativo):

Si provi ad implementare un algoritmo di ricerca degli zeri di una funzione senza utilizzare il polimorfismo. Prendere come spunto le soluzioni indicate nelle trasparenze finali della lezione teorica. Si potrebbe codificare il metodo della bisezione in una funzione che accetti in input una std::function e modellizzare la funzione di cui si vuole cercare lo zero con una funzione lambda ( oppure utilizzando un funtore o ancora utelizzando una semplice funzione).