Lezione 8 : Equazioni differenziali

In questa lezione studieremo alcuni metodi per la ricerca di soluzioni di equazioni differenziali ordinarie relative al moto di sistemi fisici interessanti. Implementeremo un codice per la risoluzione numerica di queste equazioni con i metodi di Eulero e di Runge-Kutta. Per risolvere l'esercizio vedremo come è possibile definire le principali operazioni algebriche per i vector della STL. Questo ci permetterà di realizzare i metodi di integrazione di equazioni differenziali usando una notazione vettoriale, molto simile al formalismo matematico.

ESERCIZIO 8.0 - Algebra vettoriale:

Come prima cosa proviamo a dotare i vector della STL di tutte le operazioni algebriche che ci possono essere utili definendo opportunamente gli operatori +,*,/,+= . Dal momento che non possiamo modificare gli header files della classe vector, implementiamo questi operatori come funzioni libere in un header file apposito da includere quando necessario. Provate a scrivere un piccolo codice di test per verificare il corretto funzionamento delle operazioni tra vettori : potreste provate a realizzare un semplice programma che scriva a video le componenti della risultante di due forze (vettori tridimensionali) e le componenti del versore della risultante.

Algebra per i vector :

Potete utilizzare le funzioni definite nel file VectorOperations.h. Notate che tutte le funzioni che abbiamo implementato sono dichiarate inline: in questo modo le funzioni possono essere codificate in un header file che può essere incluso in diversi file .cpp senza creare problemi in fase di linking. In questo modo i vari .cpp possono anche essere compilati separatamente e in fase di linking il compilatore utilizzerà una sola versione della funzione compilata insieme a ciascun .cpp. Potreste poi provare l'algebra dei vettori con un semplice main del tipo seguente :

#include <vector>
#include <iostream>

#include "VectorOperations.h"

using namespace std;

int main() {

   vector<double> v1 { 3.,4.,6.} ;
   vector<double> v2 { 4.,5.,6.} ;

   double para = 2.;

   // provare la somma e prodotto per scalare, calcolare il versore della risultante v1+v2

}

ESERCIZIO 8.1 - Risoluzione tramite metodo di Eulero:

Nel primo esercizio implementeremo un codice per la risoluzione numerica dell'equazione differenziale che descrive il moto di un oscillatore armonico tramite il metodo di Eulero. Gli obiettivi principali di questo esercizio sono due :

studiare l'andamento della posizione in funzione del tempo t ( per t che va da 0 a 70 secondi ) per un fissato passo di integrazione 'h (si potrebbe costruire un TGraph per esempio) e confrontare l'errore commesso con la soluzione esatta.
studiare l'andamento dell'errore che si commentte utilizzando il metodo di Eulero quando confrontiamo la soluzione approssimata con la soluzione esatta nell'istante t=70s in funzione del passo di integrazione h in un intervallo compreso tra 0.1 e 0.001.

Per testare il metodo, risolviamo l'equazione differenziale: $\frac{d}{dt} \begin{pmatrix} x \\ v \end{pmatrix} = \begin{pmatrix} v \\ -\omega_{0}^{2} x \end{pmatrix} \qquad \qquad \begin{matrix} x(0)=0 \\ v(0)=1 \end{matrix}$

mettendo in grafico il valore della x in funzione del tempo ed eventualmente anche il suo errore rispetto alla soluzione esatta del problema $x=sin(t)$ . Utilizziamo $\omega_{0}$ =1 s $^{-1}$ . Si consiglia di svolgere l'integrazione per un certo numero di periodi, in modo da vedere se l'ampiezza di oscillazione rimane costante. Integrare fino a t=70 s permette di vedere circa 10 periodi.

Il metodo di Eulero

Consideriamo la seconda legge della dimanica di Newton

$a=\frac{d^{2}x}{dt^{2}} = \frac{F}{m}$

che è un'equazione differenziale del secondo ordine che può essere ridotta ad un'equazione differenziale del prim'ordine introducendo la variabile velocità

$\begin{cases} \frac{dx}{dt} = v \\ \frac{dv}{dt} = \frac{F}{m} \\ \end{cases}$

Il metodo di Eulero consiste nel calcolare lo stato della soluzione al tempo t+h dato quello ad un tempo t tramite le espressioni:

$\begin{cases} x(t+h) = x(t) + h \frac{dx}{dt} = x(t) + hv \\ v(t+h) = v(t) + h \frac{dv}{dt} = v(t) + h\frac{F}{m} \\ \end{cases}$

Struttura del programma :

Struttureremo la soluzione del problema in modo simile a quanto fatto nelle precedenti lezioni su ricerca degli zeri e integrazione numerica:

Per primna cosa costruiamo la classe astratta FunzioneVettorialeBase che contenga il metodo virtuale puro Eval. Aggiungiamo poi la classe concreta OscillatoreArmonico nella quale implementeremo il metodo Eval concreto relativo all'oscillatore armonico.
In seguito aggiungeremo la classa astratta EquazioneDifferenzialeBase che contenga il metodo virtuale puro Passo. Aggiungiamo poi la classe concreta Eulero nella quale implementeremo il metodo Passo relativo al metodo di Eulero.

Per comodità possiamo mettere tutte queste classi nello stesso header file che potrebbe avere l'aspetto seguente :

#ifndef __EquazioniDifferenziali_h__
#define __EquazioniDifferenziali_h__

#include "VectorOperations.h" 
#include <cmath>

using namespace std;

// ===========================================================================
// classe astratta, restituisce la derivata valutata nel punto x
// ===========================================================================

class FunzioneVettorialeBase {

public:
  virtual vector<double> Eval(double t, const vector<double> & x) const = 0;
};

// caso fisico concreto, oscillatore armonico

class OscillatoreArmonico : public FunzioneVettorialeBase {
public:
  OscillatoreArmonico(double omega0) {  m_omega0 = omega0;  } ;

  virtual vector<double> Eval(double t, const vector<double> & x) const {

  // .... scrivere implementazione 

  };

 private:
  double m_omega0;  
};

// ===========================================================================
// classe astratta per un integratore di equazioni differenziali
// ===========================================================================

class EquazioneDifferenzialeBase {

public:
  virtual vector<double> Passo(double t, 
                               const vector<double>& x, 
                               double h, 
                               const FunzioneVettorialeBase &f) const =0;
};

// integratore concreto, metodo di Eulero

class Eulero : public EquazioneDifferenzialeBase {

 public:

  virtual vector<double> Passo(double t, 
                               const vector<double> & x, 
                               double h, 
                               const FunzioneVettorialeBase &f) const override {

    // ... scrivere implementazione

  };

};

#endif

Una volta implementate le classi (l'implementazione di Eulero è semplicissima se si usano le operazioni di algebra vettoriale), un possibile programma per risolvere l'esercizio è il seguente:

#include "EquazioniDifferenziali.h"

#include <iostream>
#include <string>
#include <sstream>
#include <iomanip>

#include "TApplication.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TGraph.h"

std::string convert ( double h ) ;

int main (int argc, char** argv ) {

  if ( argc!=2 ) {
    std::cerr << "Usage: " << argv[0] << " <stepsize>" << std::endl;
    return -1;
  }  

  TApplication myApp("myApp",0,0);

  Eulero myEuler;

  OscillatoreArmonico osc(1.);

  double tmax = 70.;
  double h = atof(argv[1]);
  vector<double> x {0.,1.} ;
  double t = 0.; 

  TGraph myGraph ;
  int nstep = int(tmax/h+0.5);

  // evoluzione del sistema fino a 70 s  

  for (int step=0; step<nstep; step++) {
    myGraph.SetPoint(step,t,x[0]);
    x = myEuler.Passo(t,x,h,osc);
    t = t+h;
  }

  // grafici  

  TCanvas c ;
  c.cd();
  string title = "Oscillatore armonico (Eulero h = " + convert(h) + ")" ;
  myGraph.SetTitle(title.c_str());
  myGraph.GetXaxis()->SetTitle("Tempo [s]");
  myGraph.GetYaxis()->SetTitle("Posizione x [m]");
  myGraph.Draw("AL");

  myApp.Run();
}

std::string convert ( double h ) {

  int cifre_significative = -log10(h);
  std::ostringstream streamObj3;
  streamObj3 << std::fixed;
  streamObj3 << std::setprecision(cifre_significative);
  streamObj3 << h;
  std::string strObj3 = streamObj3.str();
  return strObj3;

} ;

Notate l'uso della funzione convert che ci permette di convertire il valore di h in una string in modo che possa entrare nel titolo del grafico.

Risultati attesi

Il metodo di Eulero non è molto preciso, in effetti con un passo di integrazione modesto si vede come esso possa risultare instabile, mostrando oscillazioni la cui ampiezza varia con il tempo. La figura seguente mostra l'andamento di x in funzione di t calcolato con un passo di integrazione di 0.1 s:

Eulero01

Per avere qualcosa di anche solo visivamente accettabile dobbiamo scendere a passi di almeno 0.0001 s:

Eulero0001

La figura seguente riporta l'errore accumulato dopo 70 s di integrazione per diversi valori del passo:

ErroreEulero

Si noti come la pendenza della curva sia 1 in una scala log-log, mostrando come l'errore ottenuto sia proporzionale al passo $h$ .

ESERCIZIO 8.2 - Risoluzione tramite Runge-Kutta (da consegnare):

Ora vogliamo ripetere l'esercizio 8.1 utilizzando però il metodo di risoluzione di equazioni differenziali di Runge-Kutta (del quarto ordine) e confrontare quindi in condizioni analoghe ( $t_{max}$ = 70 s) la precisione dei due metodi.

Per svolgere l'esercizio, basterà realizzare una nuova classe concreta RungeKutta a partire da EquazioneDifferenzialeBase.

Il metodo di Runge-Kutta

Il noto metodo di Runge-Kutta è un metodo del quarto ordine ed utilizza la seguente determinazione dell'incremento:

$\begin{cases} k_{1} = \frac{dx}{dt}(t,x) \\ k_{2} = \frac{dx}{dt}(t+h/2,x+k_{1}*h/2) \\ k_{3} = \frac{dx}{dt}(t+h/2,x+k_{2}*h/2) \\ k_{4} = \frac{dx}{dt}(t+h,x+k_{3}*h) \\ \\ x(t+h) = x(t) + (k_{1}+ 2*k_{2}+2*k_{3} + k_{4})*(h/6) \\ \end{cases}$

Risultati attesi

Il metodo di Runge-Kutta del quarto ordine è molto più preciso del metodo di Eulero, infatti, anche con un passo di integrazione di 0.1 s produce oscillazioni molto stabili:

La figura seguente riporta l'errore accumulato dopo 70 s di integrazione per diversi valori del passo:

Si noti come la pendenza della curva nella sua parte iniziale sia 4 in una scala log-log, mostrando come l'errore ottenuto sia proporzionale a $h^{4}$ . Quando l'errore di troncamento del metodo diventa minore degli errori di arrotondamento della macchina si vede che non c'è più alcun miglioramento nel ridurre il passo, anzi, il maggior numero di calcoli richiesto risulta in un peggioramento globale dell'errore.

ESERCIZIO 8.3 - Moto del pendolo (da consegnare):

Incominciamo ora con una carrellata di interessanti applicazioni dei metodi che abbiamo appena studiato. In questo esercizio proveremo ad implementare la risoluzione dell'equazione del moto del pendolo: per prima cosa possiamo provare a fare un grafico dell'andamento dell'ampiezza dell'oscillazione in funzione del tempo. La cosa più interessante che possiamo studiare con questo codice è l'andamento del periodo di oscillazione in funzione dell'ampiezza di oscillazione: in questo modo possiamo verificare che per angoli grandi le oscillazioni non sono più isocrone. La struttura logica dell'esercizio dovrebbe essere la seguente :

portare il sistema in una condizione iniziale $(\theta_{0}, 0)$
far evolvere il sistema usando il metodo di Eulero o di Runge-Kutta con passo di integrazione $h$ opportuno
calcolare il periodo di oscillazione del pendolo
riportare il sistema ad una condizione iniziale con ampiezza $\theta_{0}$ variata e ripetere la sequenza di operazioni

Si suggerisce far variare $\theta_{0}$ tra 0.1 e 3 radiandi in passi di 0.1 radianti.

Il moto del pendolo

L'equazione del pendolo è data dalla relazione

$\frac{d^{2}\theta}{dt^{2}} = - \left( \frac{g}{l} \right) sin \theta$

dove $g=9.8$ m/s $^{2}$ è l'accellerazione di gravità sulla superficie terreste mentre $l$ è la lunghezza del pendolo. L'equazione differenziale si può approssimare con quella di un oscillatore armonico per piccole oscillazioni, $sin\theta \sim \theta$ . In tal caso le oscillazioni risultano isocrone cioè con periodo indipendente dall'ampiezza delle oscillazioni. Questa però è solo un'approssimazione e per grandi oscillazioni bisogna usare l'equazione esatta. Essendo il moto non più armonico, il periodo di oscillazione dipende dall'ampiezza della stessa.

Calcolo del periodo :

In questo caso l'integrazione numerica dell'equazione differenziale non si può effettuare per un tempo predefinito, ma deve essere portata avanti fino a quando non si raggiunge una condizione compatibile con l'aver terminato l'oscillazione.

Una possibile soluzione consiste nel portare avanti l'integrazione fino a quando non si registra un cambiamento di segno della velocità angolare;
Siccome possiamo calcolare la velocità solo con granularità pari al passo di integrazione, possiamo migliorare la stima del periodo di oscillazione interpolando linearmente tra i punti $(t,v(t))$ e $(t+h,v(t+h))$ calcolando quando la retta ottenuta passa per lo zero; il tempo così calcolato corrisponde al semiperiodo dell'oscillazione.

Un frammento di codice che implementa questo algoritmo è:

double A=0.1*(i+1);
double v=0.;
double t = 0.; 
vector<double> x {-A , v} ;
while ( x[1]>=0.) {
    v = x[1];    
    x = myRK4.Passo(t,x,h,osc);
    t = t+h;
    std::cout << A << " " << x[0] << " " << t << std::endl;
}
t = t-v*h/(x[1]-v);
double period = 2*t ;

Controllare la formula di interpolazione !

Risultati attesi

La figura seguente illustra il periodo al variare dell'ampiezza per un pendolo con $l=$ 1 m:

si noti come per piccole oscillazioni, il periodo sia effettivamente quello atteso dall'approssimazione dell'oscillatore armonico: $T = 2 \pi \sqrt{ \frac{l}{g} } = 2.007 s$ ma aumenti significativamente per grandi ampiezze.

ESERCIZIO 8.4 - Oscillazione forzate e risonanza (da consegnare):

Implementare la risoluzione dell'equazione di un oscillatore armonico smorzato con forzante. Fare quindi un grafico dell'ampiezza di oscillazione della soluzione in condizioni stazionarie in funzione della frequenza della forzante per ricostruire la curva di risonanza. La struttura logica dell'esercizio dovrebbe essere la seguente :

costruiamo un oscillatore armonico forzato con smorzante con una pulsazione propria $\omega_{0} = 10$ rad/se una costante $\alpha=1/30$ s
impostiamo un valore della pulsazione della forzante $\omega$ (si consiglia di esplorare un range di $\omega$ tra 9 rad/s e 11 rad/s in passi da 0.05 rad/s)
mettiamo il sistema nella sua condizione iniziale x=0 e $v_{x}$ =0
facciamo evolvere il sistema usando il metodo di Runge-Kutta con passo di integrazione h opportuno fino all'esaurirsi del transiente
calcoliamo l'ampiezza dell'oscillazione
modifichiamo il valore della pulsazione della forzante $\omega$ , riportiamo il sistema alle condizioni iniziali x=0 e $v_{x}$ =0 e ripetiamo la sequenza di operazioni

Oscillatore armonico con forzante

L'equazione dell'oscillatore armonico smorzato con forzante è data dalla relazione

$\frac{d^{2}x}{dt^{2}} = - \omega_{0}^{2} x - \alpha \frac{dx}{dt} + sin (\omega t)$

Nell'esercizio proposto utilizzare come valori iniziali $\omega_{0} = 10$ rad/s e $\alpha$ = 1/30 s.

Risultati attesi

La figura seguente riporta l'andamento dell'ampiezza in funzione del tempo nel caso in cui la frequenza propria del sistema ( $\omega_{0}$ ) e quella della forzante ( $\omega$ ) coincidano $\omega_{0} = \omega = 10$ rad/s.

È bene ricordare che per determinare l'ampiezza bisogna aspettare che il transiente delle oscillazioni si esaurisca. Questo avviene con una costante di tempo pari a $1/\alpha$ (ovvero dopo un tempo $1/\alpha$ l'ampiezza di oscillazione si riduce di un fattore 1/e).

Si consiglia quindi di integrare l'equazione differenziale per un tempo pari ad almeno dieci volte $1/\alpha$ , in modo da raggiungere una situazione in cui l'oscillazione è stabile, e poi valutare l'ampiezza. Anche in questo caso si può assumere di aver raggiunto il massimo dell'oscillazione nel momento in cui si trova un punto in cui la velocità cambia di segno.

Una curva di risonanza è illustrata in figura:

ESERCIZIO 8.5 - Moto in campo gravitazionale (facoltativo):

Implementare la risoluzione dell'equazione del moto di un corpo in un campo gravitazionale.

Verificare, nel caso del sistema Terra-Sole, che il periodo di rivoluzione della Terra intorno al Sole sia effettivamente di un anno e che l'orbita sia periodica (utilizzare un passo h di integrazione opportuno). Calcolare quindi il rapporto tra perielio ed afelio.
Provare ad aggiungere una piccola perturbazione al potenziale gravitazionale (ad esempio un termine proporzionale ad $1/r^{3}$ nella forza) e verificare che le orbite non sono più chiuse ma formano una rosetta.

Moto in campo gravitazionale

Nell'implementare il moto di un corpo in un campo gravitazionale utilizzare le seguenti condizioni:

costante gravitazionale G = 6.6742 10 $^{-11}$ m $^{3}$ kg $^{-1}$ s $^{-2}$
massa del Sole $M_{S}$ = 1.98844 10 $^{30}$ kg
distanza Terra-Sole al perielio = 147098074 km
velocità terra al perielio = 30.287 km/s

Risultati attesi

Nel caso di potenziale gravitazionale standard dovremmo ottenere una traiettoria di questo tipo :

Aggiungendo un termine di perturbazione $\alpha/r^{3}$ al potenziale gravitazionale la traiettoria diventa il seguente:

Visualizzare evoluzione temporale con ROOT :

In questo caso potrebbe essere interessante visualizzare l'evoluzione temporale della traiettoria della terra in modo dinamico. Per fare questo possiamo semplicemente visualizzare il grafico della traiettoria ogni volta che viene aggiunto un nuovo punto. Facciamo in modo che un nuovo punto venga aggiunto dopo una breve attesa. Vediamo alcuni frammenti di codice :

Per prima cosa ggiungiamo alcune librerie utili.

#include <chrono>
#include <thread>
#include "TSystem.h"

La visualizzazione dinamica si potrebbe implementare in questo modo

  // creazione TGraph e TCanvas

  TGraph tRosetta ;
  TCanvas CRosetta("CRosetta","CRosetta",600,600);
  tRosetta.GetXaxis()->SetTitle("x (m)");
  tRosetta.GetYaxis()->SetTitle("y (m)");

  for ( unsigned int npoint = 0; npoint< nstep ; npoint++ ) {

    // aggiungi il punto, disegna e aggiorna

    tRosetta.SetPoint( npoint , x[0]  , x[1] );
    tRosetta.Draw("ALP"); 
    CRosetta.Update();
    gSystem->ProcessEvents();

    // attendi 

    std::this_thread::sleep_for(std::chrono::milliseconds(1));

    // calcola il nuovo punto

    x = integratore->Passo(t,x,h,f);
    t+=h;

 }

}

La stessa struttura può essere utilizzata per tutti i casi che vogliamo studiare.

ESERCIZIO 8.6 - Moto di una particella carica in un campo elettrico e magnetico uniforme (facoltativo):

Implementare la risoluzione dell'equazione del moto di una particella carica in un magnetico uniforme. Disegnare la traiettoria della particella e determinarne il diametro dell'orbita. Cosa succede se si aggiunge un campo elettrico con componente lungo l'asse z pari a $E_{z}$ = -1000 V/m ?

Moto in campo elettrico e magnetico uniformi

Il moto di una particella carica in un campo elettrico e magnetico uniformi risente della forza di Lorentz e pertanto è descritto dall'equazione

$m \frac{d^{2} \overrightarrow{x} }{dt^{2}} = q\overrightarrow{E} + q\overrightarrow{v} \times \overrightarrow{B}$

che si può riscrivere in forma matriciale come

$\frac{d}{dt} \begin{pmatrix} x \\ y \\ z \\ v_{x} \\ v_{y} \\ v_{z} \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & \frac{q}{m} B_{z} & -\frac{q}{m} B_{y} \\ 0 & 0 & 0 & -\frac{q}{m} B_{z} & 0 & \frac{q}{m} B_{x} \\ 0 & 0 & 0 & \frac{q}{m} B_{y} & -\frac{q}{m} B_{x} & 0 \end{pmatrix} \begin{pmatrix} x \\ y \\ z \\ v_{x} \\ v_{y} \\ v_{z} \end{pmatrix} + \begin{pmatrix} 0\\ 0 \\ 0 \\ \frac{q}{m} E_{x} \\ \frac{q}{m}E_{y} \\ \frac{q}{m}E_{z} \end{pmatrix}$

Consideriamo il moto nel piano (x,y) di un elettrone in un campo magnetico costante con i seguenti valori:

q = -1.60217653 $10^{-19}$ C
m = 9.1093826 $10^{-31}$ kg
v $_{x}$ (0) = 8 $10^{6}$ m/s
B $_{z}$ = 5 $10^{-3}$ T
$E_{z}$ = -1000 V/m
tutte le altre componenti di campi e velocità iniziali sono nulle.

Questi parametri corrispondono grosso modo all'apparato sperimentale per la misura di $e/m$ del laboratorio del II anno.

Risultati attesi

Per prima cosa potremmo cercare di disegnare la traiettoria nel piano (x,y) :

e poi ( a titolo di esempio ) l'andamento della coordinata z in funzione del tempo :

è infine possibile visualizzare la traiettoria della particella in tre dimensioni :

Per rappresentare la traiettoria dell'elettrone in 3D si può utilizzare il TGRaph2D di ROOT.