Metodi numerici per i sistemi dinamici

Metodi numerici per i sistemi dinamici

Come si è visto un PLL del secondo ordine può essere descritto efficacemente tramite un sistema di equazioni differenziali non lineari del secondo ordine; le soluzioni di un sistema di questo tipo solo raramente possono essere scritte in forma esplicita, quindi nasce l'esigenza di trovare metodi alternativi per studiarne il funzionamento. Un metodo molto potente è sicuramente la simulazione del PLL per via numerica. La simulazione numerica di sistemi dinamici continui attraverso equazioni alle differenze porta a espressioni numeriche in cui le derivate vengono sostituite da metodi numerici di integrazione. Il problema principale nella simulazione numerica è la stabilità e la precisione dei risultati della simulazione, comparati con i risultati esatti.

Esistono vari metodi per la simulazione numerica e vari gradi di complessità realizzativi. Come verrà mostrato in seguito, metodi con un alto tasso di complessità danno risultati migliori in termini di stabilità e precisione, a scapito però del tempo di simulazione.

Introduzione

Prima di illustrare i metodi che verranno confrontati in questo capitolo occorre dare qualche definizione che verrà usata in seguito.

Il problema è quello di trovare un soluzione numerica del problema di Cauchy :

con

dove è il vettore delle variabili di stato, è il vettore degli ingressi.

Un generico metodo numerico ad un passo è espresso nella forma:

Dove dipende dalla funzione relativa al problema trattato e dal metodo usato, mentre è la risoluzione temporale del metodo numerico e si mantiene costante per tutta la simulazione, in modo che per ogni . Posto questa formula serve a calcolare iterativamente i valori che assumono le variabili di stato al tempo . Ovviamente questa approssimazione numerica introduce un errore ad ogni passo, da a , che può essere calcolato come la differenza tra il valore esatto e quello teorico che si ottiene usando il metodo numerico con il valore esatto , cioè :

con      

;

la quantità è chiamata errore locale di troncamento, la figura 1 spiega graficamente il significato

Figura 1 Errore locale di troncamento

Un metodo si dice coerente se vale :

Inoltre, si può definire l'ordine del metodo come il più grande intero positivo per cui vale:

.

Un'altra caratteristica importante di un metodo numerico è la convergenza che si verifica se vale la relazione :

dove è l'errore globale di troncamento al passo n+1, cioè l'errore che viene commesso dopo n+1 passi.

Figura 2 Errore globale di troncamento

Metodo di Eulero

Il più semplice procedimento per approssimare un sistema dinamico continuo con uno discreto è quello di Eulero. Le soluzioni del sistema discreto così ottenuto approssimano le soluzioni del sistema continuo. La differenza tra sistema analogico e discreto, cioè l'errore di discretizzazione, è una funzione rapidamente crescente con il passare del tempo, per cui la potenza di calcolo richiesta per ottenere soluzioni accurate è notevole.

Per ottenere il metodo di Eulero si parte dallo sviluppo in serie di X(t) in calcolato in

calcolato in

E lo si interrompe alla derivata prima, quindi, ricordando che e sostituendo con si ottiene la formula di Eulero

nella quale per comodità abbiamo omesso il vettore U(t) degli ingressi. Si nota subito che è un infinitesimo di ordine superiore al secondo rispetto a T, infatti

e quindi Eulero è un metodo del primo ordine.

In poche parole il metodo di Eulero consiste nell'aggiornare il valore di   calcolando la devirata prima nel punto , la figura 3 illustra questo concetto.

Figura 3.1 Metodo numerico di Eulero

Metodo di Runge-Kutta

Si è visto che Eulero è un'applicazione dello sviluppo in serie di Taylor in , calcolato su un passo di lunghezza T, proviamo a calcolarlo ora su due passi di lunghezza T/2:

con , per il primo passo e

per il passo intero; anche in questo caso risulta

combinando però opportunamente , possiamo eliminare l'errore di secondo ordine, infatti

Se confrontiamo l'espressione precedente con quella relativa al metodo di Eulero, osserviamo che è necessario valutare la funzione due volte per intervallo di tempo anziché una, per guadagnare un'unità sull'ordine di convergenza. Sembra quindi possibile ricavare algoritmi di ordine superiore aumentando il numero di valutazioni della funzione .

Si può però dimostrare che servono più di k valutazioni a tempi intermedi di per ottenere un ordine di convergenza uguale a k per k>4 (cfr. tabella 1) , per questo motivo l'algoritmo di calcolo più vantaggioso e anche più usato è quello di ordine 4, che prende il nome di Runge-Kutta classico o a 4 passi intermedi.

Tabella 1 Corrispondenza tra ordine e numero minimo di stati

La struttura generale dei metodi di Runge-Kutta è:

dove s è il numero di stati e con

I parametri , , sono reali e definiscono il metodo insieme a s.

Per quanto detto risulta

inoltre, tenendo conto che una condizione necessaria affinché il metodo converga e

si può dedurre che la condizione di coerenza equivale a :

Per trovare i coefficienti , , per l'algoritmo di Runge-Kutta di 4° ordine occorre imporre che risolva esattamente equazioni differenziali che abbiano come soluzioni polinomi di 4° grado, il che ci assicura che il metodo sia di ordine massimo. I coefficienti possono essere efficacemente rappresentati attraverso le tavole di Butcher. Omettendo la dimostrazione RIFERIMENTO BIBLIOGRAFICO MANCANTE si trovano i seguenti valori:

Quindi possiamo riassumere il metodo di Runge-Kutta a 4 stadi come segue

Con questo algoritmo la funzione viene valutata 4 volte per intervallo di tempo, una nel punto iniziale, due volte nel punto intermedio e ancora una nel punto finale, e i quattro valori così ottenuti vengono combinati linearmente per calcolare il nuovo valore , tutto questo è efficacemente illustrato in figura 3.2

Figura Errore. Nel documento non esiste testo dello stile specificato. Metodo di Runge-Kutta del 4° ordine

E' facile rendersi conto che il metodo di Runge-Kutta del 1° ordine coincide con il metodo di Eulero sei vengono scelti i coefficienti in modo tale che c=0, a=0 e b=1 . Come abbiamo detto per applicare R-K del 4° ordine servono più valutazioni di e, quindi, sembrerebbe che una maggior precisione viene pagata in termini di maggior onere di calcolo. Sebbene in assoluto questo sia giusto, non lo è se si ragiona in termini di precisione dell'algoritmo. Prendiamo per esempio un semplicissimo caso scalare:

che ha come soluzione e supponiamo di essere interessati al valore che assume la variabile al tempo =10s, che ha come soluzione esatta . La tabella 2 riassume i risultati ottenuti.

Tabella Confronto tra R-K del 4° ordine e Eulero

Dalla tabella si desume che nonostante il numero di calcoli a parità di passo di discretizzazione T sia maggiore nel metodo di R-K, per ottenere buoni risultati con il metodo di Eulero serve una discretizzazione decisamente più fine, il che porta ad un aumento dell'onere di calcolo da parte del simulatore.

Stabilità numerica dei metodi

La convergenza è una condizione necessaria affinché un metodo numerico produca "buoni" risultati, infatti si riferisce alla bontà dell'approssimazione su piccoli passi (), nella pratica invece è importante usare un valore del passo di discretizzazione non troppo piccolo per evitare che la simulazione diventi troppo lunga. Nasce quindi la necessità di capire per quali valori di T gli errori di discretizzazione si mantengono limitati con il procedere dei calcoli.

Questa proprietà, che riguarda la sensibilità del metodo agli errori, va sotto il nome di stabilità numerica.

Per semplicità ci possiamo riferire ad un semplice problema nel piano complesso:

con e , la cui soluzione, con c una costante che dipende dalle condizioni iniziali. Poiché , è naturale richiedere che la soluzione numerica abbia un comportamento analogo a quello della soluzione continua, cioè deve valere:

con qualsiasi, e .Prendiamo in esame il metodo di Eulero :

Affinché sia rispettata la condizione per la stabilità, deve essere verificata la disuguaglianza

quindi deve essere interno al cerchio di raggio unitario con centro [-1, 0] nel piano complesso.

Figura 3.2 Regione di assoluta stabilità del metodo di Eulero

Tutto questo può essere generalizzato per un metodo di Runge-Kutta qualunque. Cioè se applichiamo un algoritmo R-K problema , si ottiene un'equazione della forma

dove è detta funzione di stabilità. Se allora il metodo è assolutamente stabile.

Introducendo il vettore , si può dimostrare RIFERIMENTO BIBLIOGRAFICO MANCANTE che la funzione di stabilità di un metodo di Runge-Kutta è della forma:

,

che nel caso di R-K del 4° ordine si specializza nella formula:

Figura 3.3 Regione di asintotica stabilità per il metodo di Runge-Kutta del 4°ordine

Fino ad adesso ci siamo riferiti ad un caso scalare e lineare, se invece siamo interessati allo studio di un sistema dinamico non lineare a più variabili possiamo ottenere ugualmente risultati importanti sull'assoluta stabilità studiando il sistema linearizzata intorno ad un punto di equilibrio.

Prendiamo il caso

Supponiamo che sia un punto di equilibrio del sistema dinamico, nel senso che

allora in un intorno di vale la relazione

dove è la matrice Jacobiana di calcolata nel punto di equilibrio.

Prendiamo, ora, in considerazione gli autovalori di J, , possiamo dire che un metodo numerico, applicato ad un sistema dinamico non lineare, è assolutamente stabile localmente se

dove la funzione di stabilità è quella ricavata per un sistema monodimensionale lineare.

E' molto interessante notare che se siamo in presenza di un sistema dinamico lineare, il polinomio caratteristico di J , cioè ,, è nella stessa forma del denominatore della funzione di trasferimento del sistema e che quindi i valori che assumono gli autovalori di J, , corrispondono ai poli stessi del sistema.