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.