Laboratorio di Calcolo. 2019/20.

Esercitazione n. 5

1. Il programma sotto riportato integra le equazioni del moto per un oscillatore armonico in una dimensione.

   Verificare che implementi l' algoritmo di Verlet introdotto nell'  ultima  lezione.

   Provare a determinare un valore ottimale per l' intervallo  temporale di integrazione (dt). Verificare inoltre l' accuratezza dell' integrazione numerica confrontando i risultati per la posizione e velocità con i valori analitici della  soluzione (da far calcolare nel programma ed eventualmente stampare su un altro file  oppure da far calcolare da gnuplot direttamente).Si suggerisce di usare la forma A cos(omega t) + B sin(omega t) per la soluzione analitica in quanto così risulta immediato attribuire i valori alle costanti A e B in funzione dei dati iniziali  (A=x0, B=v0/omega).
   

    program harm
    implicit none
   
    real    :: massa=1.0,kappa=1.0 ! valori di default per massa e cost. elastica
    real    :: dt,ekin,epot
    real    :: pos,vel,vel_parziale,f
    integer :: nstep,it
    write(unit=*,fmt="(a)",advance="no")"delta t : "   ! il formato (a) chiede che il dato sia trattato come 
    read*,dt                                           ! caratteri e advance="no"  sopprime l'  inserimento del 
    write(unit=*,fmt="(a)",advance="no")"n.step: "     ! caratter "a capo"  alla fine della linea per cui 
    read*,nstep                                        ! la prossima operazione di lettura/scrittura inziera'
    write(unit=*,fmt="(a)",advance="no")"massa: "      ! sula stessa riga di schermo di quella corrente
    read*,massa
    write(unit=*,fmt="(a)",advance="no")"kappa: "
    read*,kappa
    write(unit=*,fmt="(a)",advance="no")"pos(0): "
    read*,pos
    write(unit=*,fmt="(a)",advance="no")"vel(0): "
    read*,vel
   
    it=0        ! step 0 : valori iniziali
    write(unit=1,fmt=*)it,it*dt,pos,vel
    epot =  0.5 * kappa * pos**2
    f    = - kappa * pos
    ekin =  0.5 * massa * vel**2
    write(unit=2,fmt=*)it,dt*it,ekin,epot,ekin+epot
   
    do it = 1,nstep
       pos = pos + vel * dt + 0.5* f/massa * dt**2
       vel_parziale = vel + 0.5 * dt * f/massa                !  prima parte della formula per le velocita'
       f    = - kappa * pos
       epot =  0.5 * kappa * pos**2
       vel = vel_parziale + 0.5 * dt * f/massa                !  la formula per le velocita' viene completata qui
       write(unit=1,fmt=*)it,it*dt,pos,vel
       ekin = 0.5 * massa * vel**2
       write(unit=2,fmt=*)it,it*dt,ekin,epot,ekin+epot
    end do
    end program harm
   
    

   Da notare che la variabile

   vel_parziale

   non è strettamente necessaria. Il valore parziale della velocità potrebbe essere immagazzinato nella stessa variabile vel che poi, dopo il ricalcolo della forza, verrebbe ridefinita:

   vel = vel + 0.5 * dt * f/massa
   f    = - kappa * pos
   epot =  0.5 * kappa * pos**2
   vel = vel + 0.5 * dt * f/massa 




2. Col programma dell' esercizio precedente studiare il caso di costante elastica k = 10.0 e massa  m = 1.0, determinando il valore ottimale per il passo di integrazione (dt) verificando la bontà della conservazione dell' energia. Va tenuto presente che per un sistema conservativo l' energia totale (energia potenziale + energia cinetica) dovrebbe essere una costante del moto. Nell' integrazione numerica la costanza dell' energia totale è verificata solo approssimativamente ad una dato ordine nell' intervallino elementare di integrazione dt. Pertanto questo va variato fino a trovare dei valori che garantiscano almeno un paio di ordini di grandezza in meno tra variazioni dell' energia totale e variazioni delle energie cinetica e/o potenziale.  Si consiglia di analizzare i dati graficamente mediante gnuplot. In particolare, il comando  plot 'fort.1' u 2:3 w l   permette di avere un grafico della legge oraria (tempo in secondi in colonna 2 e posizione in colonna 3). Il comando  plot 'fort.2' u 2:3 w l, '' u 2:4 w l, '' u 2:5 w l    grafica dati contenuti nelle colonne 3,4,5 del file fort.2 (rispettivamente energia cinetica, potenziale e totale) in funzione del tempo in secondi (colonna 2) sullo stesso grafico. L 'indicazione che i dati successivi al primo insieme provengono dallo stesso file fort.2 è contenuta nell' indicazione di file  vuota (  attenzione:  il simbolo ''  nell' esempio di comando gnuplot precedente  è due volte il carattere '  e non una volta il carattere ").
   
   Riduzioni ulteriori  del valore di dt  migliorano l' accuratezza ma a scapito dell' efficienza (troppi passi per coprire un intervallo temporale fissato).
   Cercare di verificare quantitativamente  l' andamento dell'ampiezza delle  fluttuazioni dell' energia totale come funzione del passo di integrazione dt.
   
   Graficare mediante gnuplot anche  la traiettoria nello spazio delle fasi (cioè la curva ottenuta congiungendo le coppie (pos(t),vel(t))). A che curva assomiglia ? (Suggerimento: pensare alla condizione di conservazione dell' energia) Cosa succede a questa curva se prendiamo un time step troppo grande ?


3. Modificare la massa nel programma dell' esercizio precedente in modo che valga 0.01 kg. Ripetere l' integrazione delle equazioni del moto.  Occorre variare qualche parametro per garantire la stessa precisione?




4. Gli array permettono di avere variabili collettive. Corrispondono alle variabili con indici (per ciascuno dei tipi primitivi o per ciascun derivato ). Queste possono essere pensate come insiemi di variabili etichettate con un nome comune e i cui elementi sono selezionabili mediante gli indici. Un array Fortran con un indice (array di rango 1) è definito da una istruzione

         tipo, dimension(N)    :: A

   dove tipo rappresenta il tipo dati  dell' array (ovvero real, integer, character, logical, complex) ed  N deve essere una costante intera  (eventualmente una costante simbolica intera) esprimente il numero di componenti del vettore.    Un vettore può   esser letto componente per componente con un'istruzione iterata    del tipo:

             do I=1,N
                read*,A(I)
             end do

   oppure mediante un' unica istruzione:

           read*,A

   o mediante un do implicito, utile quando non si vogliono leggere tutte le componenti dell' array, ma con la possibilità di immettere i dati sulla stessa linea:

           read*,(A(i),i=1,N-3)

   Lo stesso risultato può essere ottenuto mediante l' utilizzo di una sezione di array:

   read*,A(1:N-3)
   
       

   
   Scrivere un programma   che dichiari  un array x a 10  componenti e  ne legga le componenti, con uno qualsiasi dei metodi sopra esposti. Dopo la lettura il programma  dovrebbe scrivere su schermo l' array x due volte: prima un elemento per riga, poi tutto sulla stessa riga.
Per svolgere l'esercizio occorre tener presente che ad ogni istruzione print*,... corrisponde una nuova riga di output.

Il prodotto scalare tra vettori N-dimensionali  a e b (di componenti a_i e b_i) è calcolabile come
        ∑^N_i=1 a_i * b_i

Il modulo di un vettore è la radice quadrata del prodotto scalare del vettore con se stesso. Il coseno dell' angolo tra 2 vettori può  essere ottenuto dividendo il prodotto scalare dei due vettori per il prodotto dei moduli. Il Fortran dispone di una funzione acos(x)  che permette di ottenere l' ampiezza (in radianti) dell' angolo il cui coseno è  x (la funzione trigonometrica  arccos (x))
Inoltre le funzioni sum(x) e dot_product(x,y) calcolano la somma delle componenti dell' array x e il prodotto scalare tra i due array x e y.
Implementare un programma che

            

• legga le N componenti di a e b (definendo N come costante simbolica)
              
• calcoli i prodotti scalari  a.a, b.b e a.b sia senza, sia utilizzando le funzioni intrinseche sum (primo modo) e dot_product (secondo modo). Tenere presente che a*a è un array le cui componenti sono i prodotti delle componenti di a; la funzione sum(a) agisce su un array calcolando la somma delle sue componenti);
            
• calcoli il modulo di a e di b (usando la radice quadrata dei prodotti scalari a.a e b.b)
            
• calcoli l' angolo tra a e b in radianti.

        Verificare cosa si ottiene nel caso di tre coppie di vettori 2D di componenti (1,0) e (0,1) ,   (1,1) e (1,1)  e   (1,1) e (-1,1).



5. A partire dal programma del punto 1,  creare in un nuovo file  un programma che dia la soluzione numerica per il problema del pendolo semplice ( eq. del moto: L x'' = - g sin(x) , con L= lunghezza del pendolo e g= accelerazione di gravita' e x angolo del pendolo rispetto alla verticale, x''  denota la derivata seconda rispetto al tempo). Verificare l'  aumento del periodo di oscillazione con l'  aumento dell'  elongazione massima.  Verificare come cambiano i grafici nel piano di fase (vel  in ordinata e pos in ascissa, con i nomi delle variabili originali dell' oscillatore armonico) all'  aumentare dell'  angolo  di elongazione massima. Verificare come evolve la traiettoria all'  aumentare della velocita'  iniziale per fissato angolo iniziale. 
   Studiare anche il comportamento del pendolo nel caso di oscillazioni che partano in prossimità del punto di massimo dell' energia potenziale (punto di equilibrio instabile) corrispondente ad un angolo di π con approssimazioni per difetto crescenti: 3.14, 3.1415, 3.14159.
   Suggerimento: a parte la forma diversa del termine di forza, la struttura dell' equazione non è troppo diversa da quella per l' oscillatore armonico a patto di assegnare alla lunghezza il ruolo che aveva la massa e alla accelerazione di gravità  quello che aveva la costante di forza. Pertanto una modifica minimale del programma corrisponde a reinterpretare le variabili massa e kappa come nuovi nomi per lunghezza e accelerazione di gravità g rispettivamente, mentre pos e vel divengono l' angolo x e la corrispondente velocità angolare. Conviene anche verificare dimensionalmente le formule per la nuova posizione e la nuova velocità per controllare di aver trattato correttamente il termine con l'accelerazione. Ovviamente, per il calcolo dell' energia potenziale  e cinetica occorrono le formule corrispondenti e, se le variabili di nome "massa"  e "kappa"  sono state riutilizzate per rappresentare lunghezza e accelerazione di gravità, occorre anche introdurre una nuova variabile per la massa fisica del punto materiale vincolato a muoversi a distanza L  dal centro. Può  essere un utile ripasso di meccanica verificare che l' energia totale è data da: 
   

   L*m*( 0.5*L*vel**2 + g*(1-cos(pos))  )

1. Modificare il programma per l' integrazione numerica dell' equazione del moto dell' oscillatore armonico in modo che invece di utilizzare l' algoritmo di Verlet, usi quello di Groot e Warren (vedi note on-line) e come prima cosa verificare che il programma dia gli stessi risultati del caso Verlet per l'oscillatore armonico. Successivamente applicarlo al caso dell'equazione del moto di un oscillatore armonico smorzato e confrontare i risultati numerici con la soluzione analitica.
   
   
2. Modificare il programma per l'  integrazione delle equazioni del moto dell'  oscillatore armonico 1D per integrare quelle di un corpo che si muova   sotto l'  azione di un forza -k*x**3   (derivante da un'  energia  potenziale  k*x**4/4; k= costante positiva ).  Determinare come il periodo del moto dipende dall'ampiezza delle oscillazioni eseguendo la determinazione del periodo per elongazioni iniziali (ampiezze) diverse. In particolare cercare di individuare come tenere sotto controllo l' incertezza su queste due quantità. Come dipende dall' intervallino di integrazione dt la precisione rispettivamente del periodo e dell' ampiezza ?
   
   Anche in questo caso dare un' occhiata alla forma delle curve nel piano di fase. Come differiscono da quelle dell' oscillatore armonico ? Suggerimenti: per controllare la massima ampiezza può essere utile far partire il corpo dalla quiete. Con che precisione ritornera' al punto di partenza ? Per determinare il periodo si può partire con l' individuare gli intervalli di tempo ai cui estremi la posizione o la velocità cambiano di segno. Qual' è l' incertezza risultante sul periodo ? come si può ridurla ?