Introduzione: la funzione Cfr

La funzione Cfr presentata in questo lavoro risolve sitemi di equazioni differenziali di primo ordine con passo variabile, utilizzando i metodi di Runge-Kutta 3 e Runge-Kutta 4. I valori che la funzione riceve in ingresso sono:

�� - Sistema differenziale da risolvere

�� - Istante iniziale

�� - Istante finale

�� - Vettore dei valori iniziali

�� - Tolleranza

I valori restituiti sono:

�� - Vettore t dei tempi utilizzato

�� - Matrice con i valori delle soluzioni delle equazioni differenziali

Scarica la funzione cfr e i file di esempio

 

Controlli interni contro gli errori

La funzione, oltre a risolvere il sistema assegnato, � dotata di una serie di controlli interni relativi a situazioni critiche di calcolo. Nell�affrontare tali problemi abbiamo cercato di seguire lo stesso stile di Matlab: quando Cfr incontra un errore non si ferma, ma cerca di risolverlo modificando il passo di integrazione. Nel caso in cui ci� non sia possibile, stampa un warning, ma restituisce comunque un risultato. I due principali errori presi in considerazione sonoi:

 

1) Riduzione eccessiva di dt: in alcuni casi il passo di integrazione potrebbe diventare troppo piccolo, in termini assoluti o rispetto al vettore dei tempi t(k). Ci� potrebbe generare due errori:

- Se dt � troppo piccolo in termini assoluti, l�esecuzione del programma � troppo lenta. Abbiamo scelto come limite inferiore , che tiene conto dell�intervallo di integrazione.

- Se dt � troppo piccolo rispetto a t(k), potrebbe non essere sommabile a t(k), a causa dei limiti del sistema floating point. Il limite minimo � , ma per assicurare una somma pi� corretta abbiamo scelto .

Come soluzione a entrambi i problemi abbiamo scelto di lasciare dt al passo minimo, e procedere con l�integrazione. Viene stampato un warning indicante il mancato rispetto della tolleranza.

 

2) Problemi di somma tra u4 e e4: per u4 intendiamo l�integrale complessivo calcolato fino al passo k. Per e4 intendiamo la porzione di integrale calcolata al passo k. Questi due valori potrebbero avere ordini di grandezza molto diversi.

- : potrebbe verificarsi nel caso di funzioni con elevata pendenza. In questa situazione dt deve essere dimezzato. L'errore pu� essere corretto perch� non c'� rischio di entrare in un ciclo senza uscita.

- : potrebbe verificarsi per lunghi intervalli di integrazione. In questo caso dt dovrebbe essere raddoppiato, ma non viene corretto l'errore perch� al passo successivo il controllo sulla tolleranza potrebbe dimezzare nuovamente dt. In questo modo si entrerebbe in un ciclo senza uscita. Viene stampato un warning, e il calcolo prosegue.

 

Banchi di prova della funzione

Per testare il funzionamento di Cfr abbiamo scritto dei file (testbench.m) che utilizzano tale funzione e la funzione Ode45 per risolvere delle equazioni differenziali. Abbiamo scelto di presentare i risultati in una forma grafica, che ci sembra pi� intuitiva. Ciascun file testbench.m fornisce i seguenti risultati:

1) Grafico esatto della soluzione e della sua derivata;

2) Grafici della soluzione calcolata con Cfr e con Ode45;

3) Grafici della lunghezza dei passi per Cfr e Ode45;

4) Grafico dell'errore tra la funzione corretta e quella calcolata con Cfr e Ode45;

5) Numero di istanti calcolati e tempi di esecuzione di Cfr e Ode (tramite tic e toc).

 

La tolleranza utilizzata � sempre di . I grafici usano i seguenti colori:

Verde: soluzione dell�equazione differenziale, e sua derivata;

Rosso: risultati relativi a Cfr;

Blu: risultati relativi a Ode45.

 

Testbench A: oscillatore armonico.

In questo testbench � stato risolto il sistema differenziale:

 

 

Intervallo: [0� 2]

Tempo esecuzione Cfr: 0.0729

Tempo esecuzione Ode45: 0.2532

Numero di istanti calcolati da Cfr: 107

Numero di istanti calcolati da Ode45: 60

 

Considerazioni: Il sistema scelto intende testare il funzionamento di Cfr su funzioni regolari. Cfr svolge perfettamente il proprio compito, e meglio di Ode45: � pi� rapida nonostante utilizzi passi pi� piccoli, e rispetta meglio le tolleranze. Da notare che l�errore di Ode45 cresce col tempo, nonostante la funzione sia regolare.

 

 

 

 

Testbench B: radice

In questo testbench � stata risolta l�equazione differenziale:

���� �����

 

Intervallo: [0.00001� 5]

Tempo esecuzione Cfr: 0.0754

Tempo esecuzione Ode45: 0.1636

Numero di istanti calcolati da Cfr: 40

Numero di istanti calcolati da Ode45: 92

 

Considerazioni: L�equazione intende testare il funzionamento di Cfr su funzioni che hanno derivata prima con discontinuit� di secondo tipo (almeno uno dei due estremi tende a infinito). Sia Cfr che Ode45 presentano problemi nel caso in cui l�istante iniziale sia 0, cio� la discontinuit�, e non riescono a calcolare la soluzione. Nei grafici presentati la soluzione parte invece da valori molto vicini allo 0. Cfr impiega meno tempo e utilizza passi pi� lunghi, ma con precisione migliore di Ode45.

 

 

Testbench C: valore assoluto

In questo testbench � stata risolta l�equazione differenziale:

���� �����

 

Intervallo: [-6� 4]

Tempo esecuzione Cfr: 0.0434

Tempo esecuzione Ode45: 0.1466

Numero di istanti calcolati da Cfr: 24

Numero di istanti calcolati da Ode45: 44

 

Considerazioni: In questo tipo di funzione � presente una discontinuit� di I specie nella derivata della soluzione, e un punto angoloso nella soluzione. Sia� Cfr che Ode45 ignorano tali elementi, e proseguono il calcolo della soluzione come se non fossero presenti. L�errore, come previsto, cresce enormemente a partire dall�ultimo punto corretto calcolato prima della discontinuti�.

 

 

 

 

Testbench D: esponenziale

In questo testbench � stata risolta l�equazione differenziale:

���� �����

 

Intervallo: [0� 5]

Tempo esecuzione Cfr: 2.7301

Tempo esecuzione Ode45: 1.6990

Numero di istanti calcolati da Cfr: 1190

Numero di istanti calcolati da Ode45: 3788

 

Considerazioni: L�equazione intende testare il funzionamento di Cfr su funzioni con derivata prima stiff. Sia Cfr che Ode45 diventano poco efficienti in questo caso. Cfr risulta peggiore di Ode45 in termini di tempi di esecuzione e di istanti calcolati. Tuttavia Ode45 presenta un errore maggiore e oscillante. Ci� era prevedibile perch� questa funzione non � volta alla soluzione di problemi di questo tipo. Cfr si mantiene invece pi� stabile e rispetta la tolleranza, nonostante risenta anch�essa del problema.

 

 

 

Testbench E: iperbole

In questo testbench � stata risolta l�equazione differenziale:

 

���� �����

 

 

Intervallo: [-0.1� 0.1]

Tempo esecuzione Cfr: 2.1499

Tempo esecuzione Ode45: 0.2504

Numero di istanti calcolati da Cfr: 939

Numero di istanti calcolati da Ode45: 196

 

E� stato stampato il seguente avviso: � Warning: La tolleranza non � stata rispettata perch� il passo ha raggiunto il valore minimo imposto di 10^-6*(tf-t0) �.

 

Considerazioni: L�equazione intende testare il funzionamento di Cfr su funzioni con derivata prima con discontinuit� di seconda specie, e l�intervento dei sistemi di sicurezza e correzione dell�errore. Entrambe le funzioni riescono a superare la discontinuit� con successo. Tuttavia Cfr � molto pi� lenta di Ode45, e calcola pi� passi, ma risulta anche pi� precisa, come � possibile vedere sia dal grafico della soluzione che da quello dell�errore.

Come previsto, in Cfr interviene un meccanismo di sicurezza: per poter calcolare la soluzione in prossimit� dello 0, dove la funzione diventa molto ripida, Cfr cerca di ridurre eccessivamente il passo. Poich� interviene il sistema di limitazione su dt, il calcolo della soluzione diventa inaccurato dopo la discontinuit�.

E� stato necessario utilizzare un intervallo di integrazione molto piccolo per evitare tempi troppo lunghi di calcolo.

 

 

Conclusioni

La funzione Cfr riceve in ingresso gli stessi valori della Ode45 (a eccezione della tolleranza), e fornisce gli stessi risultati. Nel caso di funzioni regolari e senza problemi particolari (oscillatore armonico nel testbenchA) ha dimostrato di essere pi� efficiente, in termini di tempi di esecuzione, passi e tolleranza. Nel caso di funzioni stiff o con discontinuit� Cfr e Ode45 danno problemi simili. Tuttavia Ode45 non si blocca mai, e spesso da risultati migliori(esponenziale nel testbenchD), mentre in altri prevale Cfr (logaritmo nel TestBenchE). Ancora una volta emerge il legame tra efficienza e semplicit�: Ode45 permette di affrontare meglio un insieme pi� ampio di funzioni diverse, anche se sul gruppo ristretto di funzioni regolari Cfr � pi� efficiente.

Spesso l�errore commesso da Cfr, soprattutto per le funzioni regolari, � di molto inferiore alla tolleranza. Possiamo spiegare ci� con la scelta, in fase di programmazione, di utilizzare sempre i valori di Runge Kutta 4 per la soluzione: tali valori devono infatti essere sempre calcolati per controllare il passo di integrazione. Essendo pi� precisi di Runge-Kutta 3, sono stati utilizzati come soluzione.

Durante la scrittura della funzione sono emerse vari possibili miglioramenti all�algoritmo, che non abbiamo per� implementato. Tra questi citiamo:

- Modifica variabile del passo: in Cfr il passo � sempre dimezzato o raddoppiato. Si potrebbe pensare a un algoritmo che aumenti o diminuisca il passo di una quantit� ottimale rispetto alla situazione.

- Gestione delle discontinuit�: una sezione di Cfr che sia in grado di riconoscere le funzioni irregolari e permetta di dividere l�intervallo di integrazioni in sottointervalli privi di discontinuit�.

 

In conclusione, possiamo affermare che Cfr integra perfettamente funzioni regolari, con scarsa efficienza problemi stiff, e genera errori nel caso di funzioni discontinue.