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:
_file/image014.gif)
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.
_file/image017.gif)
_file/image018.gif)
Testbench
B: radice
In questo testbench � stata risolta l�equazione
differenziale:
����
�����_file/image024.gif)
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.
_file/image025.gif)
_file/image026.gif)
Testbench
C: valore assoluto
In questo testbench � stata risolta l�equazione
differenziale:
����
�����_file/image031.gif)
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�.
_file/image032.gif)
_file/image033.gif)
Testbench
D: esponenziale
In questo testbench � stata risolta l�equazione
differenziale:
����
�����_file/image038.gif)
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.
_file/image039.gif)
_file/image040.gif)
Testbench
E: iperbole
In questo testbench � stata risolta l�equazione
differenziale:
����
�����_file/image045.gif)
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.
_file/image046.gif)
_file/image047.gif)
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.