Laboratorio di Calcolo 2

proff. A. Andreazza, D. Galli, E. Spoletini, G. Tiana, A. Vairo

Universita` degli Studi di Milano

Anno Accademico 2005/2006

Lezione 5

Introduzione

Scopo di questa sessione di laboratorio e' utilizzare le classi di ROOT TMatrixD e TVectorD per svolgere alcune operazioni con matrici. In particolare le applicheremo allo studio di alcuni sistemi dinamici.

TMatrixD e TVectorD

Nella sintassi standard del C++ si possono devinire in generale dei vettori di un qualunque oggetto, ad esempio un vettore di file:

ifstream f[17];

Allo stesso modo si possono definire delle matrici multidimensionali di qualunque tipo di oggetto:

TH1F istogramma[5][5];

Siccome un vettore od una matrice sono strutture molto generali, il C++ non puo' prevedere su esse delle funzioni matematiche specifiche.

Tuttavia nei problemi di fisica spesso abbiamo a che fare con vettori e matrici con elementi numerici che operano secondo le regole dell'algebra lineare. In questi casi abbiamo a che fare con oggetti appartenenti ad una struttura matematica particolarmente ricca, contenente diverse operazioni, ad esempio il prodotto scalare di vettori o la moltiplicazione riga per colonna tipica delle matrici.

Per questo ha senso definire degli oggetti piu` complicati dei semplici vettori forniti dal C++: dei vettori numerici dotati della possibilita` di effettuare le operazioni matematiche che ci interessano. Nel nostro caso ROOT ci fornisce le classi:

TVectorD: vettore di numeri reali in doppia precisione
TMatrixD: matrice bidimensionale di numeri reali in doppia precisione

In queste dispense verranno indicati i metodi principali di queste classi: una descrizione completa si trova nella referenza di ROOT.

Notate che la classe TMatrixD e' una classe derivata dalla classe TMatrixDBase, per cui puo` usare i metodi di tale classe. Ad esempio, il metodo UnitMatrix(), che inizializza la matrice ad essere una matrice diagonale unitaria, e' tra i metodi di TMatrixDBase e non di TMatrixD.

TVectorD

Per dichiarare un vettore si puo` usare un costruttore che ha come unico argomento un intero indicante la dimensione del vettore:

TVectorD x(3); // x e' un TVectorD di dimensione 3

Si noti che il costruttore e' una funzione, che viene chiamata al momento dell'esecuzione del programma, e quindi ìl suo argomento puo` avere una dimensione variabile.

Ad esempio, il seguente frammento di codice si comporta come ci si aspetta:

int N;
cout << "Inserisci la dimenzione del vettore" << endl;
cin  >> N;
TVectorD vettore(N);  // OK, il vettore e' dimensionato correttamente

mentre il seguente codice e' sbagliato, anche se un compilatore tollerante a volte lo permetta:

int N;
cout << "Inserisci la dimenzione del vettore" << endl;
cin  >> N;
double vettore[N];    // NO! Non si sa a cosa il vettore e' stato dimensionato

Per accedere alle componenti del vettore, si utilizza l'operatore () (parentesi tonde). Funziona anche l'operatore [] (parentesi quadre), ma per similitudine con le matrici che verranno dopo, e' meglio utilizzare le parentesi tonde. Un vettore puo` quindi venire inizializzato ad esempio con:

TVectorD x(3);
x(0)=1.;
x(1)=0.;
x(2)=-1.;

Al momento della creazione, un vettore puo` anche essere inizializzato direttamente al valore di un altro vettore, sia usando un costruttore, che un'inizializzazione esplicita. I seguenti costrutti che inizializzano il nuovo vettore z ad avere lo stesso contenuto del vettore x definito in precedenza, sono del tutto equivalenti:

TVectorD z(x);
TVectorD z=x;

Ora segue l'elenco di alcuni metodi utili:

void Print(): stampa il contenuto del vettore
int GetNRows(): restituisce il numero di componenti del vettore
double Norm2Sqr(): calcola la norma del vettore (somma dei quadrati degli elementi)
double Sum(): effettua la somma degli elementi del vettore
double Max(): resituisce il valore piu` alto nel vettore
double Min(): resituisce il valore piu` basso nel vettore

I seguenti metodi invece alterano il contenuto del vettore, trasformandolo secondo le seguenti operazioni:

TVectorD& Abs(): per ogni elemento v(i)=|v(i)|
TVectorD& Sqr(): per ogni elemento v(i)=v(i)*v(i)
TVectorD& Sqrt(): per ogni elemento v(i)=sqrt(v(i))
TVectorD& Invert(): per ogni elemento v(i)=1./v(i)

TMatrixD

Analogamente ad un vettore, una matrice bidimensionale puo` venire dichiarata con dimensioni arbitrarie nel costruttore (eventualmente calcolate al run-time). Gli elementi della matrice si possono poi accedere utilizzando sia gli operatori () che []. Il seguente programmino calcola il triangolo di Tartaglia:


int dim;
cout << "Inserisci la dimensione del triangolo da calcolare:" << endl;
cin  >> dim;
TMatrixD triang(dim,dim); // matrice quadrata con lato pari a dim
triang(0,0)=1.                 // costruiamo la prima riga
for (int i=1; i<dim; i++) { // esemplificando l'accesso 
  triang[0][i]=0.;             // agli elementi di matrice usando
}                              // sia l'operatore () che l'operatore []
for (int j=1; j<dim; j++) {
  triang(j,0)=1.;
  for (int i=1  ; i<=j ; i++) triang(j,i)=triang(j-1,i-1)+triang(j-1,i);
  for (int i=j+1; i<dim; i++) triang(j,i)=0;
}

Infine, per verificare la struttura di una matrice sono a disposizione i seguenti metodi:

void Print(): stampa il contenuto della matrice
int GetNrows(): restituisce il numero di righe della matrice
int GetNcols(): restituisce il numero di colonne della matrice

Anche per le matrici si puo` effettuare l'inizializzazione come copia di un'altra matrice nel costruttore:


TMatrixD M(triang);


TMatrixD M=triang;

sono inizializzazioni equivalenti. I seguenti metodi invece alterano il contenuto della matrice, sostituendola con quella trasformata:

TMatrixD& Invert(double *determinante);: calcola l'inversa della matrice e mette in determinante il valore del determinante della matrice originale. Se non si vuole vedere il determinante, si puo` omettere l'argomento
TMatrixD& T(): calcola la matrice trasposta.
TMatrixDBase& UnitMatrix(): imposta la matrice ad essere una matrce identita` (chiaramente funziona solo se la matrice originaria e` quadrata).
TMatrixDBase& Zero(): azzera la matrice.

N.B.: il valore di ritorno e` un riferimento ad oggetti TMatrixD, in modo da permettere delle assegnazioni del tipo:
TMatrix nuova=vecchia.T();
ma bisogna ricordarsi che in questa operazione, anche la matrice vecchia viene modificata. L'inizializzazione corretta richiede di fare prima una copia della matrice originale:
TMatrix nuova=vecchia;
nuova.T();

Come esempio di verifica del funzionamento delle operazioni di cui sopra, potete provare ad eseguire la seguente catena di operazioni e verificare che dia il risultato atteso:


TMatrixD A(3,3);
A.UnitMatrix()
A.Print()
A(1,0)=0.5;
A(2,2)=0.5;
A.Print();
double det;
A.Invert(&det);
A.Print()
cout << det << endl;
A.T();
A.Print();

Algebra lineare

Con le classi TVectorD e TMatrixD si possono effettuare alcune operazioni di algebra lineare in maniera molto mnemonica.

Supponiamo di aver definito le seguenti variabili:

TVector x(N), y(N);
TMatrix A(M,N), B(I,M);
double a;

Allora:

x*y

restituisce un double pari il prodotto scalare tra x ed y

a*x

restituisce un TVectorD pari alla moltiplicazione del TVectorD x per lo scalare a. L'operazione x*=a e` a identica a x=a*x.

x+y

restituisce n TVectorD pari alla somma dei TVectorD x ed y. L'operazione x+=y e` a identica a x=x+y.

x-y

restituisce n TVectorD pari alla differenza dei TVectorD x ed y. L'operazione x-=y e` a identica a x=x-y.

A*x, x*A

prodotti destro e sinistro di matrice per vettore

x*=A

modifica il vettore x effettuando l'operazione di moltiplicazione matrice per vettore Ax. Si noti che:

la matrice A sembra apparire al posto sbagliato (dopo la x anziche' prima),
x viene modificato dalla moltiplicazione. Oltre a cambiare il suo valore, viene cambiata anche la sua dimensione (dopo l'operazione e' M e non piu` N).

B*A

restituisce una TMatrixD di dimensione IxN, pari alla moltiplicazione righe per colonne di B per A. Per matrici quadrate e` anche definita l'operazione B*=A, identica a B=B*A.

Equazioni differenziali (lineari, a coefficienti costanti)

Esistono numerosi problemi fisici significativi descritti da equazione del moto del tipo:
dx/dt=Ax+b

dove A e' una matrice e x e b sono dei vettori contenenti le componenti della velocita` e della posizione di un punto.

Un esempio classico e' quello dell'oscillatore armonico smorzato, la cui equazione:
d2x/dt2+adx/dt+w2x=0
puo` essere resa in forma matriciale come:
d/dt(x,v)=(0,1//-a,-w2>)(x,v)

Se al sistema aggiungiamo una forzante esterna:
d2x/dt2+adx/dt+w2x=F/m
e lo mettiamo in forma matriciale:
d/dt(x,v)=(0,1//-a,-w2)(x,v)+(0,F/m)
possiamo identificarlo con l'equazione generale precedente
A=(0,1//-a,-w2) b=(0,F/m)

Un altro esempio e' il moto di una particella in un campo elettromagnetico costante:
d2x/dt2=(q/m)(E+vxB)
che puo` venire espressa nella forma:
A=(0,0//0,B), b=(0,E)

Punti fissi

In un'equazione del tipo
dx/dt=Ax+b

sono dei punti fissi quelli per cui la derivata rispetto al tempo e' nulla. Questi sono le soluzioni dell'equazione
$dx/dt=Ax+b=0 x=A^{-1}b$

Se la matrice A non fosse invertibile, questo significherebbe che il moto presenta una deriva (si considerino ad esempio le matrici corrispondenti al moto rettilineo uniforme o uniformement accelerato, ma anche quella della particella in campo elettromagnetico uniforme).

Risoluzione dell'equazione differenziale

Un'equazione differenziale di questo genere si puo` risolvere numericamente con il metodo di Eulero, ovvero stimando il valore del vettore x ad un certo tempo, dal suo valore ad un tempo precedente piu` un incremento dato dalla derivata della funzione moltiplicata per la differenza di tempo tra i due punti:
x(t+dt)=x(t)+(Ax+b)dt

Il seguente frammento di codice effettua l'integrazione dell'equazione dell'oscillatore armonico smorzato, per un tempo di 10 secondi e valori di dt di 0.01 s e stampa un grafico della traiettoria nello spazio delle fasi (x,v):


double a, w, F;
cout << "inserisci il valore di omega: ";
cin  >> w;
cout << "inserisci il valore dello smorzamento: ";
cin  >> a;
cout << "inserisci il valore della forzante: ";
cin  >> F;
TMatrixD A(2,2);
A.Zero();     // Azzera la matrice
A(0,1)=1.;    // inizializza la matrice A
A(1,0)=-w*w;
A(1,1)=-a;
TVectorD b(2);
b(0)=0.;     // inizializza il vettore b
b(1)=F;
TVectorD x(2); // condizioni iniziali:
x(0)=1.;       // il punto e` a distanza 1 dall'origine
x(1)=0.;       // ed a velocita` nulla
TGraph traiettoria; // grafico dove mettere la traiettoria
traiettoria.SetPoint(0,x(0),x(1)); // primo punto: condizioni iniziali
double dt=0.001, tmax=10.;
int nstep=tmax/dt;
for (int i=1; i<=nstep; i++) {
  TVectorD h(x);
  h*=A;  // h=Ax 
  h+=b;  // h=Ax+b
  h*=dt; // h=(Ax+b)*dt
  x+=h;  // x(t+dt)=x(t)+(Ax+b)*dt
  traiettoria.SetPoint(i,x(0),x(1));
}
traiettoria.Draw("ALX");

Il problema grosso del metodo di Eulero e' che l'approssimazione al primo ordine non e' molto buona e bisogna trovare un buon compromesso tra la precisione ottenibile ed il passo di integrazione. Negli esercizio che seguono, provate ad aumentare il tempo di integrazione su di una soluzione nota per provare a vedere quanto varia con il tempo l'errore di integrazione.

Esercizi

Oscillatore armonico smorzato

Considerando l'equazione dell'oscillatore armonico smorzato con:
α=2
F/m=0.8
calcolare i punti fissi e provare ad vedere la forma delle traiettorie nello spazio delle fasi per diverse intensita` della forza dell'oscillatore armonico:

ω=0.7
ω=1.0
ω=2.0

Come vi spiegate la forma qualitativamente diversa tra i vari casi?

Oscillatori armonici accoppiati

Generalizzare la trattazione precedente al caso di oscillatori armonici accoppiati:
oscillatori armonici accoppiati

In particolare, provate tracciare le traiettorie nel piano (x,y) in alcuni casi particolari:

β=0, &gamma=0: scegliere le pulsazioni tali che:
- ω_x/ω_y sia un intero;
- ω_x/ω_y sia un numero razionale intermedio tra 1 e 2, come 3/2, 4/3...;
- ω_x/ω_y sia un numero irrazionale.
In tutti i casi dovete integrare per un tempo abbastanza lungo da includere molti periodi per entrambe le frequenze.
ω_x=ω_y=1:
- β=γ=0.9, 1., 1.1;
- β=-γ=0.1, -0.1;
Per saggiare i vari fenomeni, provate a sperimentare piu` di un condizione iniziale. Ad esempio sono interessanti i punti (1,1) ed (1,-1) con velocita` iniziale nulla ed il punto (1,0) con velocita` iniziali (0,1) e (0,-1).

Deriva di una particella in un campo e.m. uniforme

Considerate il moto nel piano (x,y) di una particella soggetto ad un campo elettromagnetico tale che:

q/m E_x=1,
q/m B_z=1,
tutte le altre componenti del campo nulle.

Potete prendere la particella inizialmente ferma in un punto. Vedrete che tendera` ad acquistare una velocita` di deriva a causa del campo elettrico che tende ad accelerarla, ma... va nella giusta direzione?