Lezione 2

Introduzione

In questa sessione di laboratorio introdurremo la risoluzione delle equazioni differenziali ordinarie con il metodo di Eulero e ne approfitteremo per introdurre alcuni concetti legati alla precisione numerica delle soluzioni.

Nei prossimi esercizi poi svilupperemo altri metodi e ne confronteremo la precisione con quella del metodo di Eulero.

Dal punto di vista delle proprieta` del linguaggio C++, approfondiremo il concetto di gestione dinamica della memoria nel caso di chiamate a funzione ed introdurremo il concetto di puntatore a funzione.

Il metodo di Eulero

Se consideriamo la legge della dinamica di Newton
a=d²x/dt²=F/m
questa e' un equazione differenziale del secondo ordine che puo` venire tradotta in un'equazione differenziale del primo ordine, introducendo come variabile ausiliaria la velocita`:
dx/dt=v
dv/dt=F/m

Il metodo di Eulero consiste nel calcolare la lo stato della soluzione ad un tempo t+h, dato quello ad un tempo t, tramite la relazione approssimata:
x(t+h)=x(t)+h*dx/dt=x(t)+h*v
v(t+h)=v(t)+h*dv/dt=v(t)+h*F/m

Il programma

In questo esercizio indroduttivo affronteremo per semplicita` solo problemi monodimensionali e quindi scriveremo:

1. un programma principale che chieda come dati di ingresso:
   1. la posizione e la velocita` di un punto materiale al tempo t=0
   2. il tempo totale per il quale si vuole integrare l'equazione
   3. il valore di h da usare.
   ed e poi stampi una serie di terne di valori di t, x e v con t che va da 0 al tempo massimo ad intervalli h.
2. diverse funzioni che soddisfino il prototipo:
   double* nomefunzione(double *x, double t);
   dove x e' un vettore contenente in x[0] la posizione ed in x[1] la velocita`, e la funzione restituisce un vettore che contiente nella locazione 0 la derivata della posizione e nella locazione 1 la derivata della velocita`. Ad ogni equazione differenziale corrispondera` una funzione diversa. Abbiamo previsto una variabile aggiuntiva t per le forze dipendenti dal tempo.
3. infine la funzione che implementa il metodo di Eulero, cioe` che, data una funzione per calcolare le derivate, fornisce le condizioni al tempo t+h date quelle al tempo t:
   void eulero(double& t, double h, int N, double *x,double* derivata(double *x, double t))
   dove all'uscita della funzione il vettore x contiene i nuovi valori e la variabile t e' stata incrementata della quantita` h.

Abbiamo scelto questa separazione del codice perche' nei prossimi esercizi vorremmo poter riutilizzare questo codice usando algoritmi diversi da quello di Eulero.

L'aggiunta del parametro N nella dichiarazione della funzione ci permettera` poi di riutilizzarla anche per equazioni differenziali non puramente monodimensionali. Per equazioni al second'ordine monodimensionali, avremo che N=2.

Qui di seguito seguono delle implementazioni di esempio di tali funzioni con alcuni commenti sui trucchetti che trovate e che verranno discussi in dettaglio a lezione

La funzione principale

eqdifferenz.cxx

Si noti, come i prompt per l'input siano diretti verso lo standard error cerr, anziche' cout.

Questo e' voluto, in modo da non avere scritte nel file di output qualora si ricorra alla ridirezione dello standard output. Di solito i due stream di uscita cerr e cout vanno entrambi al terminale e si mischiano tra loro. La ridirezione separa i due stream, mantenento il primo sul terminale ed il secondo sul file prescelto.

Infine si noti che nella chiamata a eulero, abbiamo passato non solo delle variabili, ma anche una funzione. L'utilizzo di puntatori a funzione sara` descritto meglio a lezione.

La funzione che descrive l'equazione differenziale

oscillatore.cxx

Come esempio e` implementata una funzione che descrive un oscillatore armonico.

Si noti che, se vogliamo restituire un vettore contenente dei dati, dobbiamo allocarlo usando l'operatore new. Se avessimo definito derivata semplicemente con:

double derivata[2];

Allocando della memoria invece, la variabile locale derivata di tipo double* viene distrutta all'uscita della funzione, ma non la memoria allocata, che rimane accessibile anche successivamente.

Per poter utilizzare questa funzione nel programma principale, ricordarsi di costruire anche il file oscillatore.h.

La funzione eulero

eulero.cxx

eulero vuole come argomento un puntatore ad una funzione con il prototipo definito nella sezione precedente. Si noti la sintassi particolare da usare quando si vuole passare come argomento non una variabile, ma una funzione.

Si noti l'utilizzo dell'operatore delete per liberare la zona di memoria allocata all'interno della funzione oscillatore.

Se non lo facessimo, ad ogni iterazione il programma richiederebbe sempre maggiore memoria fino a quando non finisce con l'esaurire la memoria disponibile sul PC, con gravi conseguenze. Questo tipo di problema di chiama memory leak ed e' uno degli aspetti cui prestare molta attenzione al momento della scrittura di programmi in C++.

Per toccare con mano questo fenomeno, una volta impostato l'esercizio, provate a cancellare il delete ed a far partire l'integrazione d ell'equazione differenziale per un tempo lungo. Mentre questa gira, aprite un'altra finestra terminale e date il comando top. Dovreste trovare il vostro processo che gira prendendosi una grossa parte del tempo di CPU, e con una memoria occupata che aumenta gradualmente con il tempo di esecuzione.

Se reinserite il delete e ripetete la procedura, dovrete osservare che invece la frazione di memoria occupata rimane costante.

Anche in questo caso, per utilizzare questa funzione, avete bisogno di creare il file eulero.h.

L'esercizio

Preparare gli header file e scrivere un makefile per compilare le implementazioni di esempio riportate sopra o quelle che avete scritto voi.

Risolvere numericamente l'equazione differenziale per un oscillatore armonico di frequenza unitaria, integrando su 10 periodi partendo dalla condizioni iniziali:
x=0
v=1 var
iando il valore di h da 0.1 a 10^-5.

Fare un grafico della funzione x(t) risultante. Dovreste ottenere figure simili alle seguenti:

h=10^-2

h=10^-4

Dovreste osservare che la precisione aumenta al diminuire di h, come ci si aspetta dallo sviluppo al primo ordine sfruttato per il metodo di Eulero.

Per mettere in evidenza la dipendenza dell'errore da h e dal tempo, possiamo sfruttare il fatto che conosciamo la soluzione di questa equazione differenziale:
x(t)=A*cos(2πx+Φ)
per valori di A e Φ opportuni a seconda delle condizioni iniziali.

Fare un grafico della differenza tra i valori veri della soluzione e quelli calcolati numericamente, per i casi di cui sopra. Dovreste ottenere dei grafici del tipo:

h=10^-2

h=10^-4

Anche in questo caso, ecco delle macro di esempio:

• grafico.C
• errore.C

Il metodo di Eulero migliorato

Una versione migliorata del metodo di Eulero consiste nel determinare l'incremento usando la derivata prima della funzione non nel punto di partenza, ma nel punto di mezzo tra quello di partenza e quello stimato dal metodo di Eulero.

Dobbiamo quindi creara una nuova funzione eulerom, con la stessa interfaccia di eulero, ma che realizza l'incremento del punto secondo la seguente relazione:
k₁=dx/dt(t,x)
k₂=dx/dt(t+h,x+k₁*h)
x(t+h)=x(t)+(h/2)*(k₁+k₂)
dove si vuole sottolineare che x e k sono dei vettori che contengono gli incrementi di tutte le variabile dell'equazione differenziale del primo ordine corrispondente (quindi posizioni e velocita`).

Ha senso comparare il metodo di Eulero e quello migliorato, a parita` di peso computazionale. Siccome il metodo di Eulero migliorato richiede di valutare le derivate prime ad ogni passo dovremmo comparare i due metodi nel caso in cui il passo del metodo di Eulero migliorato e` doppio di quello del metodo di Eulero.

Provare a ripetere gli ultimi passi dell'esercizio precedente con h=0.02 e h=0.002 e confrontare gli errori con quelli del metodo di Eulero per passi h=0.01 e h=0.001. Verificate anche l'ordine di consistenza del metodo.

Metodo Runge-Kutta del 4^o ordine.

Il piu` noto metodo di Runge-Kutta del 4^o ordine utilizza la seguente determinazione dell'incremento:
k₁=dx/dt(t,x)
k₂=dx/dt(t+h/2,x+k₁*h/2)
k₃=dx/dt(t+h/2,x+k₂*h/2)
k₄=dx/dt(t+h,x+k₃*h)
x(t+h)=x(t)+(h/6)*(k₁+2k₂+2k₃+k₄)

Dobbiamo quindi creara una nuova funzione rungekutta, con la stessa interfaccia di eulero, ma che realizza l'incremento del punto secondo la precedente relazione. Siccome ora il numero di valutazioni della funzione e` 4, valutatene la precisione per h=0.04 e h=0.004 e confrontatela con le condizioni analoghe per i metodi precedenti, valutando anche l'ordine di consistenza del metodo.

Applicazioni

Quelli che seguono sono alcuni dei problemi che di fisica del primo e secondo anno che potete provare a risolvere numericamente, usando i metodi di integrazione di cui sopra. Affrontatene completamente almeno uno.

La difficolta` sta nel saper trovare un passo di integrazione sufficientemente piccolo per verificare la stabilita` del risultato (ad esempio si potrebbe integrare con un certo passo e poi con uno dimezzato, e vedere che il risultato non cambia)

Altra piccola difficolta` e` l'aumento di dimensione per alcuni problemi.

Infine ci sono situazioni in cui la vostra equazione differenziale dipende da un parametro e volete poterlo cambiare dall'interno di un programma. Un metodo e' quello di utilizzare delle variabili globali all'interno di un file e definire delle funzioni che le cambiano. Ad esempio, il file
osc_forzata.cxx
definisce la funzione osc_forzata che implementa l'equazione differenziale di un oscillatore armonico con una forzante di ampiezza unitaria e frequenza variabile attraverso la funzione set_frequenza. Chiaramente adesso avremo un header file
osc_forzata.h

Per inciso, ecco un'altra ragione per dividere i nostri programmi in diversi file: la variabile globale f e` nota solo alle funzioni osc_forzata e set_frequenza, quindi tutti le altre funzioni che costruiremo potranno usare delle variabili con lo stesso nome senza creare conflitti.

Pendolo

L'equazione differenziale esatta per un pendolo e` data dalla relazione
d<sup>2</sup>θ/dt<sup>2</sup>=-(g/l)sinθ
dove g=9.8 m/s<sup>2</sup> e` l'accelerazione di gravita` sulla superficie terrestre e l e` la lunghezza del pendolo (supponiamo che sia 30 cm).

Questa equazione differenziale diventa quella di un oscillatore armonico per piccole oscillazioni per cui sinθ~θ, risultando in oscillazioni isocrone, ovvero con periodo indipendente dall'ampiezza.

Fare un grafico del periodo in funzione dell´ampiezza di oscillazione e verificare che per grandi angoli le oscillazioni non sono piu` isocrone.

Oscillazioni forzate e risonanza

Impostate l'equazione di un oscillatore armonico smorzato con una forzante:
d²x/dt²=-ω₀²x-αdx/dt+sinωt
e fate un grafico dell'ampiezza dell'oscillazione stazionaria in funzione della frequenza dell'oscillatore, ricostruendo la curva di risonanza.

Come valore dei parametri potete usare quelli del sistema sperimentale massa/molla del laboratorio del primo anno:

• ω₀= 10 rad/s
• α= (1 /30 s)

Notate che, a seconda delle condizioni iniziali che mettete, dovrete aspettare che si esaurisca il transiente.

Moto in campo gravitazionale

Impostare l'equazione differenziale per il moto di un corpo in un campo gravitazionale.

Usando le seguenti condizioni:

• costante di gravitazione=6.6742 10^-11 m³kg^-1s^-2
• massa del sole=1.98844 10³⁰ kg
• distanza terra-sole al perielio=147 098 074 km
• velocita` al perielio = 30.287 km/s

Aggiungere una piccola perturbazione al potenziale gravitazionale (ad esempio un termine nella forza proporzionale a 1/r³) ed osservare che ora le orbite non sono piu` chiuse, ma formano una rosetta.

Moto di una particella carica in un campo elettrico e magnetico uniforme

Una particella carica in un campo elettromagnetico risente di una forza di Lorentz data dall'equazione:

che si puo` esprimere in forma matriciale come

Consideriamo il moto nel piano (x,y) di un elettrone in un campo magnetico costante, tale per cui

• q=-1.60217653 10^-19 C
• m=9.1093826 10^-31 kg
• v_x= 8 10⁶ m/s
• B_z=5 10^-3T
• tutte le altre componenti dei campi e della velocita` iniziale sono nulle.

Disegnare la traiettoria della particella e determinarne il diametro.

Aggiungere alla situazione precedente un campo elettrico:

• E_x=10000 V/m

Soluzione di equazioni algebriche

• Metodo di Dirichlet
  per equazioni della forma x=f(x)
  partendo da una soluzione di prova x₀ si utlizza la successione
  x_n+1=f(x_n)
  fino a quando ||x_n+1-x_n|| diventa minore di una tolleranza fissata.
• Metodo della secante
  per equazioni della forma f(x)=b
  partendo dagli estremi x₀ e x₁ dell'intervallo in cui si va a cercare la soluzione, si utilizza la successione
  x_n+1=x_n-1+(b-f(x_n-1))(x_n-x_n-1)/(f(x_n)-f(x_n-1))
  nuovamente iterando fino a quando la distanza tra due elementi successivi diventa minore della tolleranza.

Per applicare il metodo di Dirichlet, può essere utile riscrivere l'equazione per la i-esima radice:
x=atan(x)+(i-1)π

Relazione