SOMMARIO PPUNTI DI SOLUZIONE NUMERICA DI CIRCUITI A ...unina.stidue.net/Teoria dei...

M. de Magistris Appunti di Modelli numerici per circuiti - 22/06/2009

APPUNTI DI

MODELLI NUMERICI PER CIRCUITI

M. DE MAGISTRIS

A.A. 2008-2009

PARTE II:

ALGORITMI PER LA SIMULAZIONE CIRCUITALE


SOMMARIO

1. SOLUZIONE NUMERICA DI CIRCUITI A-DINAMICI: SISTEMI DI EQUAZIONI ALGEBRICHE ........................................................................ 3

SOLUZIONE NUMERICA DI SISTEMI ALGEBRICI LINEARI ............................................. 4 Metodo di Gauss .................................................................................................4 Sistemi triangolari e fattorizzazione LU .............................................................7

SOLUZIONE NUMERICA DI SISTEMI ALGEBRICI NON LINEARI ..................................... 9 Metodo di Bisezione ..........................................................................................10 Algoritmo di NewtonRaphson e simili.............................................................11 Condizionamento e criteri di arresto degli algoritmi iterativi ..........................15 Valutazione dellerrore dellalgoritmo di Newton-Raphson.............................18 Linearizzazione di elementi circuitali ed algoritmo di Newton Raphson..........19 Estensione dellalgoritmo di Newton-Raphson al caso n-dimensonale ............21 Algoritmi di punto fisso.....................................................................................22

2. SOLUZIONE NUMERICA DI CIRCUITI DINAMICI: INTEGRAZIONE DI EQUAZIONI ODE ....................................................................................... 27

CONSIDERAZIONI SULLERRORE NUMERICO............................................................ 28 ALCUNI ALGORITMI DI INTEGRAZIONE TEMPORALE ................................................ 30

Algoritmi di Eulero ...........................................................................................31 Algoritmo dei trapezi ........................................................................................33 Algoritmo di Gear .............................................................................................34

PROPRIET DI CONSISTENZA ................................................................................... 35 ZERO-STABILIT E TEOREMA DI LAX-RICHTMYER.................................................. 37 ASSOLUTA STABILIT (CASO LINEARE MONO-DIMENSIONALE)).............................. 38 ASSOLUTA STABILIT (CASO NON LINEARE MONO-DIMENSIONALE) ....................... 41 CONVERGENZA E STIMA DELLERRORE NUMERICO ................................................. 43

Convergenza ed errore di troncamento per Eulero esplicito ............................44 Errore di round off per Eulero esplicito .......................................................48

ALGORITMI RUNGE KUTTA..................................................................................... 50 3. INTEGRAZIONE DI EQUAZIONI ODE: CASO N-DIMENSIONALE .... 53

DIAGONALIZZAZIONE DELLE EQUAZIONI DINAMICHE ............................................. 54 CRITERI DI ASSOLUTA STABILIT NEL CASO N-DIMENSIONALE ............................... 55

Eulero esplicito .................................................................................................56 Eulero implicito.................................................................................................57 Trapezi ..............................................................................................................58 Algoritmo di GEAR ...........................................................................................59 Caso non lineare n-dimensionale......................................................................60

CENNI A CIRCUITI STIFF ED ALGORITMI A PASSO VARIABILE................................ 61

M. de Magistris Appunti di Teoria dei circuiti - 22/06/2009

1. Soluzione numerica di circuiti a-dinamici: sistemi di equazioni

algebriche Questa parte del corso dedicata ad introdurre i principali metodi

numerici che si incontrano nella simulazione circuitale. Abbiamo gi avuto modo di osservare che per circuiti non lineari, o che comunque abbiano un elevato numero di elementi, non pensabile di effettuare lanalisi con tecniche analitiche. Abbiamo anche osservato che nella discretizzazione di un problema circuitale essenzialmente sono da affrontare tre questioni fondamentali: la soluzione numerica di sistemi di equazioni algebriche lineari, la soluzione numerica di sistemi di equazioni algebriche non-lineari, la soluzione numerica di sistemi di equazioni differenziali.

Introdurremo dapprima il problema della soluzione numerica di circuiti lineari e non lineari a-dinamici. I metodi numerici che possono essere adoperati ricadono in due generali categorie, che sono i metodi diretti e quelli iterativi. Nei primi, com intuibile, lalgoritmo consta di un numero prefissato di passi che conduce alla soluzione esatta, a meno degli errori di approssimazione dovuti alla precisione limitata dellaritmetica numerica. Per converso, un metodo si dice iterativo quando la soluzione viene ricercata basandosi su di una successione di iterazioni, che in un numero non determinato a priori (ed eventualmente infinito) di passi ci conduce alla soluzione del problema. Pi realisticamente, in un numero finito (e possibilmente piccolo) di passi ci conduce ad una buona approssimazione della soluzione. La base teorica per il funzionamento dei metodi iterativi il principio delle contrazioni, cui accenneremo nel seguito.

Nel risolvere numericamente un problema sempre necessario poter stimare lerrore che si commette, vuoi a causa dellalgoritmo, vuoi per le conseguenze della rappresentazione dei numeri con un numero limitato di cifre. E altres importante saper valutare il costo computazionale dellalgoritmo che si vuole utilizzare.

4 Soluzione numerica di circuiti a-dinamici: sistemi di equazioni algebriche


Soluzione numerica di sistemi algebrici lineari In questo paragrafo riprendiamo brevemente le problematiche relative alla

soluzione numerica di sistemi di equazioni algebriche lineari, con cenni ai problemi di onere computazionale e di precisione, anche in considerazione della rappresentazione dei numeri con un numero finito di cifre.

Consideriamo il sistema lineare:

A =x b , (1.1)

con A matrice quadrata di ordine N, x vettore delle incognite e b vettore dei termini noti. Sappiamo dallalgebra che esso ammette ununica soluzione se la matrice A risulta invertibile, ovvero a rango pieno (rank(A)=N).

La soluzione del sistema (1.1) pu essere espressa attraverso la regola di Kramer:

11 12 1, 1 1 1, 1 1

1 2 , 1 , 1

..... ..... ..............................

..... ...., 1,2... ,

A A

+

+= = =

k k N

N N N k N N k NNkk

a a a b a a

a a a b a ax k N (1.2)

La formula di Kramer per di scarso interesse pratico in quanto il calcolo dei determinanti richiede unonere computazionale dellordine di (N+1)! che risulta assolutamente proibitivo per sistemi anche di modesta dimensione. Ad esempio, per un sistema di 20 equazioni in 20 incognite, e considerando un calcolatore da 1Gflops/s (floating point operations) si ha: 9 10 321!/10 5.1 10 s 1.6 10 anni= t

Metodo di Gauss Un altro metodo, ben noto, quello di eliminazione (di Gauss-Jordan),

che usiamo spesso anche nella soluzione manuale dei sistemi. Esso si basa sullidea di ridurre il sistema (1.1) ad uno equivalente (cio con la stessa soluzione) ma in forma triangolare:

U =x b , (1.3)

Metodo di Gauss 5


dove U una matrice triangolare superiore (Upper) e b un nuovo termine noto. Ci si ottiene mediante operazioni di riga sul sistema (1.1) che hanno appunto la propriet di non alterare la soluzione. Il metodo pu facilmente essere descritto in relazione allesempio:

(1) (1) (1) (1)11 1 12 2 13 3 1(1) (1) (1) (1)21 1 22 2 23 3 2(1) (1) (1) (1)31 1 32 2 33 3 3

+ + =+ + =+ + =

a x a x a x ba x a x a x ba x a x a x b

(1.4)

Al primo passo la riga 1 venga moltiplicata per 21 11/a a e sommata alla seconda, ed ancora moltiplicata per 31 11/a a e sommata alla terza per eliminare 1x dalle ultime equazioni, ottenendo:

(1) (1) (1) (1)11 1 12 2 13 3 1

(2) (2) (2)22 2 23 3 2(2) (2) (2)32 2 33 3 3

+ + =+ =+ =

a x a x a x ba x a x ba x a x b

(1.5)

dove posto (1)1

1 (1)11

= iiama

, si ha (2) (1) (1)1 1ij ij i ia a m a= e (2) (1) (1)

1 1= i i ib b m b .

Ai passi successivi il procedimento viene iterato fino ad ottenere un sistema triangolare:

(1) (1) (1) (1)11 1 12 2 13 3 1

(2) (2) (2)22 2 23 3 2

(3) (3)33 3 3

+ + =+ =+ =

a x a x a x ba x a x b

a x b (1.6)

ovvero in generale:

( ) ( )A =N Nx b . (1.7) Il sistema nella forma (1.7), nella quale lultima equazione di immediata soluzione, viene risolto facendo il percorso a ritroso con successive sostituzioni (Backward Substitutions) e trovando cos il valore di tutte le variabili.



possibile mostrare che in questo caso lonere computazionale dato da 3 2 32 22

3 3+ N N N . Considerando ancora una volta un sistema di 20

equazioni in 20 incognite, ed un calcolatore da 1Gflops/s si ha: 3 9 620 /10 5.3 10 s.= t

utile osservare che egli elementi ( )kkka nella (1.6) sono quelli per i quali

sono divisi i coefficienti della matrice nella procedura di eliminazione delle incognite. Per tale motivo prendono il nome di elementi Pivot (perno) e devono risultare non nulli perch la procedura possa proseguire correttamente. Consideriamo il seguente esempio:

1 2 3 1 2 3

1 2 3 1 2 3

1 2 3 1 2 3

1 1 1 1 1 1 1 11 1 2 2 0 0 1 11 2 2 1 0 1 1 0

+ + = + + =+ + = + + =+ + = + + =

x x x x x xx x x x x xx x x x x x

la cui soluzione esatta : 1 2 3 1x x x= = = . Lapplicazione diretta

dellalgoritmo porta a ricavare, nel primo passaggio immediatamente x3=1. Analogamente sottraendo la prima equazione alla terza otteniamo x2+x3=0. Se il nostro algoritmo fosse scritto in questo modo molto ingenuo la procedura, a questo punto, sarebbe gi andata in crisi, in quanto alliterazione successiva si dovrebbe dividere per il coefficiente (2)22 0=a . Il problema per pu essere risolto (banalmente) scambiando le ultime due equazioni prima di effettuare il passo successivo.

Vogliamo ora mettere in evidenza ulteriori problemi associati con la soluzione del sistema in esame. Consideriamo allora una piccola perturbazione sui coefficienti del sistema:

1 2 3

1 2 3

1 2 3

1 1 1 11 1.0001 2 21 2 2 1

+ + =+ + =+ + =

x x xx x xx x x

Risolvendo con il procedimento appena illustrato, otteniamo per eliminazioni successive, le equazioni:

Sistemi triangolari e fattorizzazione LU 7


1 2 3

1 2 3

1 2 3

1 +1 +1 10 0.0001 1

0 0 9999 10000

=+ + =

+ + =

x x xx x x

x x x

e la soluzione risulta 1 2 31.0001, 1.0001, 1,0001x x x= = = , che in effetti vicina alla soluzione del sistema senza perturbazione sui coefficienti. Si potrebbe anche dire che la perturbazione sui coefficienti produce una perturbazione della soluzione dello stesso ordine di grandezza!

Supponiamo ora di ripetere le stesse operazioni ma con unaritmetica con un limitato numero di cifre significative dopo la virgola, per esempio n=3. Dallultima equazione ricaviamo:

310000 1.0001 1.0009999

= = =x

E sostituendo a ritroso nelle precedenti equazioni si ha la soluzione

1 2 30, 0, 1,000x x x= = = , che ora per differisce da quella originaria sulla prima e no sullultima cifra significativa!. Dunque osserviamo in questo caso che, pur avendo unaritmetica a quattro cifre significative, lerrore che commettiamo addirittura sulla prima cifra. Tali problemi, che dimostrano come possa essere delicato il trattamento numerico delle equazioni, possono essere superati mediante la tecnica del Pivot. Essa consiste individuare ad ogni passo i-esimo allinterno della colonna i-esima lelemento di valore massimo, effettuando poi la riduzione successiva rispetto a questultimo. Nel caso appena considerato, ad esempio, invece di usare le equazioni nel loro ordine quando si tratta di ridurre la seconda equazione andiamo a vedere qual lelemento pi grande e scambiamo la seconda equazione con la terza. Il risultato finale, sempre in aritmetica a tre cifre dopo la virgola, sar in tal caso

1 2 31.000, 1.000, 1,000x x x= = = . Si tratta ancora di una soluzione approssimata, che differisce dalla soluzione vera sulla quarta cifra decimale e lapprossimazione molto migliorata. Cio, attraverso un utilizzo pi meditato ed opportuno del metodo siamo riusciti a riportare lerrore nella soluzione allinterno dellerrore con cui rappresentiamo i numeri.

Sistemi triangolari e fattorizzazione LU Abbiamo gi osservato come nel metodo di Gauss si ottiene un sistema

triangolare, la cui soluzione data per sostituzione diretta (allindietro) delle



variabili calcolate gi calcolate, a partire dallultima. In generale considerato un sistema nella forma (triangolare inferiore o Lower):

11 1 1

21 22 2 2

31 32 33 3 3

0 00 ,

=

l x bl l x bl l l x b

(1.8)

la sua soluzione espressa da:

( )( )

1 1 11

2 2 21 1 22

3 3 31 1 32 2 33

= = =

x b lx b l x lx b l x l x l

(1.9)

ed il corrispondente algoritmo prende il nome di Forward substitution.

Per lalgoritmo di Gauss abbiamo gi visto come al termine delle iterazioni id riduzione si ottiene una matrice triangolare alta ( ) =NA U . Considerando allora le matrici di trasformazione gaussiana, che esprimono in forma matriciale le operazioni di riga per ogni iterazione di riduzione:

( 1) ( ) ,+ =k kkA M A (1.10)

si ha:

( )

( )1 1

11 1

........ ,

........

NN

N

A U M M A

A M M U LU

= =

= =

(1.11)

E possibile mostrare che la matrice L triangolare inferiore (Lower), a

partire dal fatto che le matrici kM sono tali. In tal caso il sistema diviene:

=LUx b (1.12)

La fattorizzazione =A LU non che un modo di reinterpretare lalgoritmo di Gauss; come tale non altera il costo computazionale. per importante osservare che una volta note L ed U, la soluzione del problema di partenza si ottiene risolvendo in successione i due sistemi triangolari:

Sistemi triangolari e fattorizzazione LU 9


, ,= =L Uy b x y (1.13) Pertanto la maggior parte dellonere computazionale risiede nella

fattorizzazione. Ci comporta consistenti vantaggi quando ci si trovi a dover risolvere pi volte lo stesso sistema con una diversa colonna di termini noti, cosa che come vedremo piuttosto frequente nella simulazione circuitale

Soluzione numerica di sistemi algebrici non lineari La soluzione di sistemi algebrici non lineari costituisce un sotto-problema

fondamentale nellanalisi generale dei circuiti, ed al tempo stesso un problema ben assestato nella matematica numerica. Al contrario del caso precedente (sistemi lineari) possiamo dire che essa generalmente si basa su algoritmi di tipo iterativo, come in precedenza descritti. Considerata lequazione:

( ) 0=f x (1.14)

e la sua soluzione:

: ( ) 0;=x f x (1.15)

un qualsiasi metodo iterativo basato sulla definizione di una successione:

( ) ( ): lim .

=k k

kx x x (1.16)

Val la pena osservare che lutilizzo di un qualsiasi metodo iterativo per la

soluzione di un problema, una volta che sia stata definita la legge che a partire da un certo valore alla k-esima iterazione ci permette di calcolare il corrispondente alla iterazione successiva, pone alcune questioni che meritano di essere indagate, ed in particolare: - la successione ( )kx deve essere ben definita, ovvero da qualsiasi

iterazione k deve essere possibile calcolare in modo univoco il punto successivo ( 1)kx + ;

- la successione ( )kx deve essere convergente, ovvero deve tendere alla soluzione esatta in (al pi) un numero infinito di passi.



Diremo che un algoritmo globalmente convergente se la convergenza non dipende dal valore di inizializzazione della successione. Viceversa diremo che solo localmente convergente se ci non verificato.

Per un generico algoritmo di tipo iterativo, se definiamo lerrore alla k-esima iterazione:

( ) ,= kke x x (1.17)

si pu definire lordine di convergenza come:

1lim ,+

=kpkk

eC

e (1.18)

dove p lordine di convergenza e C il fattore (costante asintotica) di convergenza.

Daremo conto di tali questioni almeno nel caso degli algoritmi di nostro maggiore interesse, assieme naturalmente alle considerazioni sullonere computazionale e sulla stima dellerrore.

Metodo di Bisezione Considerato un intervallo [ ],a b in cui la ( )f x (supposta continua) ha

segno opposto agli estremi, per il teorema degli zeri esiste almeno una soluzione dellequazione ( ) 0=f x allinterno dellintervallo. Sulla base di ci si pu costruire un semplice algoritmo che dimezza di volta in volta lintervallo, individuando il sottointervallo ai cui estremi la ( )f x ha ancora segno diverso (a meno del caso che in un estrmo si annulli, nel qual caso si trovata la soluzione. Sinteticamente si ha:

( 1) ( ) ( 1) ( ) ( ) ( )

( 1) ( ) ( 1) ( ) ( ) ( )

( 1) ( 1) ( 1)

, se ( ) ( ) 0, se ( ) ( ) 0

( ) / 2

k k k k k k

k k k k k k

k k k

a a b x f x f aa x b b f x f bx a b

+ +

+ +

+ + +

= = = . (1.25)

Condizionamento e criteri di arresto degli algoritmi iterativi 15


Infine, dovremmo considerare il caso in cui esistano pi zeri della funzione ( )f x . In tal caso, ancora una volta, le propriet di convergenza nonch la soluzione alla quale eventualmente si converga dipenderanno dalla scelta del punto iniziale.

Sintetizzando ora i casi sin qui considerati, a proposito della convergenza dellalgoritmo di Newton, possiamo dire che possono verificarsi le seguenti circostanze:

1. lequazione ( ) 0=f x ha ununica soluzione lalgoritmo converge verso la soluzione lalgoritmo non converge 2. lequazione ( ) 0=f x ha pi soluzioni lalgoritmo converge verso una delle soluzioni (ad esempio quella pi

vicina alla stima iniziale (0)x ) lalgoritmo non converge

In entrambi i casi il problema della scelta del punto iniziale cruciale sia per la convergenza, che per la determinazione della soluzione nel caso ve ne sia pi di una.

Dal punto di vista pratico spesso si associa allalgoritmo di Newton qualche algoritmo pi semplice (ad es. bisezione) o per la stima del punto iniziale, o per sfuggire a situazioni di non convergenza locale.

Condizionamento e criteri di arresto degli algoritmi iterativi Vogliamo ora discutere la questione dei criteri di arresto. Ricordiamo

infatti che, teoricamente, in un numero infinito di iterazioni possiamo convergere alla soluzione esatta. Ma daltro canto, non possibile immaginare di realizzare un numero illimitato di iterazioni nella realt. Dobbiamo allora stabilire dei criteri di arresto delliterazione, che potranno essere genericamente di due tipi:

- valutazione del residuo, ovvero: ( )( ) kf x , - valutazione dellincremento (assoluto) di x: ( ) ( 1)

Condizionamento e criteri di arresto degli algoritmi iterativi 17


Esempio 1: sistemi lineari. Consideriamo come problema il classico sistema lineare A =x b , dove il vettore b assume il ruolo di dato. Si ha 1 ( )A G d= =x b e 1( )G d A = , per cui il numero di condizionamento relativo risulta:

1 1 1 1( ) ( )A A

K d A A A A A K A = = =b x xx x x

(1.32)

Esempio 2: equazione non lineare. Consideriamo come problema lequazione non lineare: ( , ) ( ) ( ) 0.F x d f x x d= = = (1.33) In tal caso 1( ) ( )x d G d= = e considerato che lequazione scalare si ha

anche che ( )1 ( ) 1 ( )d x = ; infine osserviamo che per la (1.33) ( ) ( )x f x = . Pertanto per i numeri di condizionamento, osservando che in

questo caso le norme si riducono ai valori assoluti, avremo:

( ) 1 ( ) ; ( ) 1 ( ) .absd

K d f x K d f xx

(1.34)

Dunque per il problema considerato la valutazione della derivata nel punto in esame fornisce una stima del condizionamento. Possiamo ora applicare quanto appena visto alla soluzione mediante un qualsiasi algoritmo iterativo del problema ( ) 0f x = al fine di valutare laccuratezza dei criteri di arresto. Consideriamo con le usuali notazioni x la soluzione del problema, ( )kx la soluzione approssimata alla iterazione k, e lerrore ( ) ( )k ke x x= . Inoltre definiamo il residuo ( ) ( )( )k kr f x= . Se supponiamo di fissare un criterio di arresto delliterazione basato sul residuo: ( ) ( )( )k kr f x = (1.35) dalla (1.34) abbiamo:



( ) ( )

( )( )

.( ) ( )

k kk

e f xe

x f x x f x

(1.36)

immediato concludere che:

a. se ( )( ) 1 kf x e

b. se ( )( ) 1 kf x e >> (criterio inaffidabile) I risultati appena riassunti sono illustrati dal punto di vista grafico nella seguente figura.

Figura 1.5: interpretazione grafica dei criteri di arresto per lalgoritmo di Newton-

Raphson

Valutazione dellerrore dellalgoritmo di Newton-Raphson Vogliamo, ora, considerare una questione piuttosto rilevante, che la valutazione dellerrore delalgoritmo. Per far ci, consideriamo ( )f x mediante lo sviluppo di Taylor, attorno a un generico punto x. Per lo sviluppo considerato il termine di Lagrange (al secondo ordine) esiste un punto , con

( , )x x , per il quale si ha leguaglianza:

21( ) ( ) ( )( ) ( )( )2

= + + f x f x f x x x f x x . (1.37)

Linearizzazione di elementi circuitali ed algoritmo di Newton Raphson 19


Ricordando ora che per ipotesi ( ) 0f x = , posto ( )= kx x si ha immediatamente:

( ) ( ) ( ) ( ) 21( ) '( )( ) ''( )( ) 02

+ + =k k k kf x f x x x f x x , (1.38)

da cui, dividendo tutto per ( )'( )kf x , si ricava facilmente

( )( ) '' ( ) 2

( )' ( ) ' ( )

( 1)

( ) 1 ( )( )( ) 2

k kk

k k

k

f x f x xx xf x f x

x

+

+ =

. (1.39)

Osserviamo che il termine tra parentesi a primo membro rappresenta

proprio ( 1)kx + ; pertanto, se definiamo la quantit ( )kke x x= come residuo (errore) alla k-esima iterazione, avremo:

( 1) 21 ( )1 ''( )2 '( )

kk k k

fx x e ef x

++ = = (1.40)

Siccome, se il metodo convergente si ha che ( )kx x e x quando k ; dunque possiamo affermare che:

121 ''( )lim .2 '( )

k

kk

e f x Cf xe

+

= = (1.41)

Ci che esprime il fatto che il metodo di Newton ha, come si dice, convergenza quadratica.

Linearizzazione di elementi circuitali ed algoritmo di Newton Raphson Lalgoritmo appena descritto passibile di una istruttiva interpretazione circuitale, la quale permette anche di comprendere meglio la modalit operativa di utilizzo di tali algoritmi allinterno di simulatori circuitali come ad esempio SPICE. A tale scopo riprendiamo nuovamente lo sviluppo in serie



di Taylor (per una funzione f che descrive questa volta la caratteristica di un generico resistore non lineare), e come precedentemente lo interpretiamo in forma iterativa: ( 1) ( ) ( ) ( 1) ( )( ) ( ) ( )( )k k ' k k kf x f x f x x x+ += + ; (1.42) Supponiamo, per fissare le idee, che il bipolo in questione sia un resistore non lineare con una caratteristica controllata in tensione, come mostrato in Figura 1.2a. Lequazione (1.42), particolarizzata al caso in esame, diviene: ( 1) ( ) ( ) ( 1) ( )( )( )k k ' k k ki i g v v v+ += + (1.43) Essa pu essere facilmente interpretata come la caratteristica di un generatore reale di corrente, come indicato in Figura 1.2b.

Figura 1.6: un bipolo controllato in tensione ed il relativo circuito equivalente discreto

alla k-esima iterazione Considerato dunque un elemento resistivo non lineare, possiamo interpretare la sua linearizzazione in relazione a due punti successivi della iterazione di Newton Raphson come appena visto. Se il circuito presenta pi elementi non lineari, ad ogni passo di iterazione ciascuno di essi viene sostituito dal corrispondente equivalente di Thevenin o di Norton. Tali generatori sono caratterizzati da valori di tensioni e di correnti determinati ai passi precedenti e dunque noti. Va a questo punto osservato esplicitamente che se invece di Newton Raphson si applica il metodo delle corde il valore della g non varia da una iterazione allaltra, e ci significa dover risolvere, nelle iterazioni successive, sistemi lineari con la stessa matrice dei coefficienti e con nuovi vettori di termini noti. Ci pu essere sfruttato per ridurre lonere computazionale quando si utilizza la fattorizzazione LU per la soluzione del sistema lineare.

Estensione dellalgoritmo di Newton-Raphson al caso n-dimensonale 21


Estensione dellalgoritmo di Newton-Raphson al caso n-dimensonale E istruttivo a questo punto andare verificare le problematiche che sorgono nellestensione al caso n-dimensionale dellalgoritmo di Newton. In forma vettoriale lequazione (1.14) diviene: ( ) =f x 0 (1.44) Considerato un generico punto 0x anche in questo caso possiamo sviluppare in serie la ( )f x :

( ) ( )0

0 0( ) ( )=

= +

x x

f xf x f x x x

x, (1.45)

dove questa volta ( )f x

x la matrice Jacobiana della funzione ( )f x .

Lalgoritmo iterativo pu formalmente essere scritto con la stessa forma del caso scalare, rappresentata dallespressione:

( )

1( 1) ( ) ( )( ) ( )

k

k k k

+

=

= x x

f xx x f xx

(1.46)

Nellimplementazione dellalgoritmo, si costretti in questo caso ad invertire la matrice Jacobiana ad ogni passo di iterazione, operazione onerosa che rende lento il metodo. Un modo equivalente pi conveniente riscrivere la (1.46) nella forma:

( )

( )( )

( 1) ( ) ( )

( 1) ( ) ( )

( ) ( ),

,

kk

k k k

k k k

+

=

+

=

= +

x xx

f x x x f xx

x x x

(1.47)

dove ad ogni iterazione va risolto un sistema lineare con matrice dei

coefficienti ( )

( )k

A=

=

x x

f xx

.



In ambedue le formulazioni il calcolo della matrice Jacobiana e la sua inversione (o equivalentemente la soluzione del problema lineare nella (1.47)) hanno un onere computazionale tanto pi rilevante allaumentare dellordine del sistema (1.44). Spesso pertanto si preferiscono schemi come quello delle corde gi introdotto, nel quale la matrice Jacobiana quella relativa al primo passo di iterazione, e dunque il costo computazionale delle iterazioni si abbassa molto non dovendo essere essa invertita o equivalentemente risolto il problema (1.47) ad ogni passo. Naturalmente ci viene pagato con un ordine di convergenza pi basso.

Algoritmi di punto fisso La maggior parte delle tecniche iterative per la soluzione di equazioni algebriche non lineari (ivi compreso lalgoritmo di NewtonRaphson) possono essere considerati come casi particolari dei cosiddetti algoritmi di punto fisso. Il principio di base molto semplice: si consideri infatti lequazione: ( )x F x= (1.48) ed, al solito, sia x~ una soluzione per la (1.48) ( ( )x F x= ). Essa viene detta anche punto fisso dellequazione, in quanto, se essa viene considerata in senso iterativo, quando si sostituisca in una certa iterazione il valore x in

( )F x si riottiene lo stesso punto; in altri termini, se pongo nx x= ottengo ancora 1nx x+ = .

Leggendo al solito la (1.48) in senso iterativo, ovvero 1 ( )n nx F x+ = , supponiamo che vi sia convergenza; ci vuol dire che, considerato un punto di partenza arbitrario 0x , 1x rappresenta unapprossimazione di x migliore rispetto ad 0x , ovvero 0 1x x x x >

Gli algoritmi di punto fisso hanno uninterpretazione geometrica molto semplice ed utile. Infatti, considerata la curva ( )F x , il punto fisso altro non che lintersezione tra questultima e la retta y x= .

In tale grafico, literazione facilmente rappresentabile; difatti, considerato un generico punto iniziale 0x , si trova il corrispondente valore

( )F x (primo tratto verticale), e lo si proietta sulla retta y x= (primo tratto orizzontale) per trovare (secondo tratto verticale) il corrispondente 1x . A questo punto si itera il procedimento.

Algoritmi di punto fisso 23


Figura 1.7

Nel caso considerato facile verificare che si converge verso la soluzione (che rappresentata dallintersezione delle due curve) indipendentemente dalla scelta del tentativo iniziale. In realt possono presentarsi anche altri casi, come graficamente mostrato nelle Figura 1.8 e Figura 1.9.

Figura 1.8



Figura 1.9

Vogliamo ora chiederci quali sono le condizioni per le quali si ha certamente che, posto 1 ( )n nx F x+ = si abbia la convergenza, ovvero: lim nn x x = . (1.49)

Esse sono fornite dal seguente Teorema di Banach Caccioppoli, detto anche principio delle contrazioni:

Se ( )F x una contrazione di N in N , ovvero

( ) ( ) L F x F y x y con NL 1, ,< x y (1.50) allora ( )F x ha un unico punto fisso x cui converge la sequenza delle iterate 1 ( )n n+ =x F x Inoltre lerrore alla iterazione nma limitato da:

1 01

n

nL

L

x x x x (1.51)

Il teorema appena enunciato fornisce le condizioni, che dipenderanno dalla

funzione ( )F x , per le quali garantita la convergenza. Naturalmente non sempre, nella soluzione di un problema circuitali, possibile analizzare a priori tali propriet.



Osserviamo infine che lalgoritmo di Newton definito dalla (1.22). pu facilmente essere interpretato nella classe degli algoritmi di punto fisso a patto che si ponga:

( )( )'( )

f xF x xf x

= . (1.52)

Naturalmente in tal caso le propriet di convergenza, ed eventualmente

lunicit dipenderanno dalle propriet di F(x), dunque in ultima analisi da f(x).


2. Soluzione numerica di circuiti dinamici: integrazione di

equazioni ODE Vogliamo ora introdurre il problema della soluzione numerica dei circuiti

dinamici, estendendo quindi al caso pi generale lanalisi numerica di un circuito. Anche in questo caso, cos come abbiamo fatto in precedenza il caso dei circuiti a-dinamici, lobiettivo sar quello di dare una panoramica sulle problematiche e sui principali algoritmi in uso, anche in vista del loro utilizzo allinterno dei simulatori circuitali pi diffusi.

Ricordiamo che il problema che vogliamo risolvere numericamente esprimibile nella forma:

0 0,

( , ),( )

tt=

=

x f xx x

(2.1)

e supporremo che esso risulti ben posto e stabile; ricordiamo anche che condizione sufficiente per lunicit della soluzione del problema (2.1) che la

( , )tf x risulti globalmente Lipschitziana. Queste premesse sono fondamentali quando ci accingiamo ad applicare al problema (2.1) un qualsivoglia schema numerico. evidente infatti che lo schema numerico, che deve comunque essere costituito da una successione ben definita, produce ununica successione di numeri che ha possibilit di fornire unapprossimazione della soluzione del problema (cio convergere alla soluzione del problema (2.1)) solo nel caso in cui tali ipotesi siano verificate.

Nelle considerazioni che seguono al posto del problema (2.1) faremo riferimento dapprima al caso scalare:

0 0

( , ),( ) .

x f x tx t x=

= (2.2)

28 Soluzione numerica di circuiti dinamici: integrazione di equazioni ODE


Ricordiamo che il problema differenziale (2.2) equivalente a quello integrale:

0

0( ) ( ) ( , ) .t

tx t x t f x t dt= + (2.3)

Tale forma equivalente risulta spesso pi comoda per costruire schemi numerici di soluzione, che difatti prendono il nome di algoritmi di integrazione temporale.

Considerazioni sullerrore numerico

Prima di addentrarci nella descrizione di qualsivoglia algoritmo numerico che permetta di analizzare (in modo approssimato) la dinamica dei circuiti, vale la pena fare alcune considerazioni sui diversi tipi di errore che, in generale, commettiamo nella soluzione numerica di un problema dinamico. Supponiamo allora che )(~ tx sia la soluzione esatta del problema (2.2), dunque il nostro riferimento. Qualsiasi schema numerico di soluzione prevede preliminarmente la discretizzazione dellintervallo temporale di interesse in una successione di istanti kt , distanziati di intervalli (piccoli) 1k k kt t t = . In particolare considerato un generico intervallo ( )0 , kt t t= e per semplicit assumiamo costante il passo di discretizzazione dellintervallo ponendo

0kk

t tt tk

= = e 0kt t k t= + .

Figura 2.1: discretizzazione dellasse dei tempi



Con le definizioni date ( )k kx t x rappresenta la soluzione esatta calcolata

nel generico istante kt ; naturalmente nessun metodo numerico consente di trovare esattamente i valori )(~ ktx . Se chiamiamo kx i valori ottenuti al passo k-esimo mediante il generico algoritmo che vogliamo considerare, possiamo facilmente definire lerrore globale ke nel modo seguente:

.k k ke x x= (2.4) Nel caso in cui si verifica che

0lim 0,kt e = si dice che lalgoritmo

considerato convergente. Naturalmente tale condizione fondamentale per una qualsivoglia utilit pratica dello schema numerico: un algoritmo convergente pu approssimare bene quanto si desideri la soluzione vera del problema quando si scelga t sufficientemente piccolo!

Lerrore globale definito dalla (2.4) (complessivo) ke per sua natura consta di due parti concettualmente distinte: una parte associata allalgoritmo in se, che intrinsecamente basato su approssimazioni (se cos non fosse potremmo disporre della soluzione esatta al problema (2.2)) e viene in genere chiamata errore di troncamento (il motivo di tale denominazione sar pi chiaro nel seguito); unaltra parte, logicamente disgiunta dalla precedente, invece associata al fatto che la rappresentazione dei numeri nel calcolatore con un numero limitato di cifre, e viene chiamato errore di arrotondamento (round off). Osserviamo esplicitamente che questultimo, essendo associato alla lunghezza della parola macchina, dipender a parit di algoritmo dal calcolatore che utilizziamo. Ci non vale invece per lerrore di troncamento, che dipende solo dallalgoritmo scelto. Per tali motivi i due tipi di errore vengono anche detti, rispettivamente, errore algoritmico ed errore di macchina.

Un possibile andamento dellerrore a titolo di esempio rappresentato in Figura 2.2. Lerrore, al crescere del tempo (e quindi del numero di passi di iterazione), va in questo caso ad accumularsi. Ha senso dunque definire, accanto allerrore globale ke , un errore locale ke :

k k ke x x= , con 1 1k kx x = . (2.5)



Figura 2.2: errore globale scomposto nelle sue aliquote di troncamento e di

arrotondamento

La differenza tra errore globale e locale risiede nel fatto che questultimo rappresenta lerrore che si commette in un singolo passo, a partire dallo stesso valore a quello precedente, mentre quello globale in qualche maniere ingloba anche la porzione di errore che si accumulata nei passi precedenti. Naturalmente anche allerrore locale pu applicarsi la precedente scomposizione in errore di troncamento e di arrotondamento.

Se lerrore locale tende a zero al tendere a zero di t si dice che il metodo numerico considerato consistente al problema (2.2), nel senso che localmente lapprossimazione pu essere migliorata quanto si vuole quando si scelga t sufficientemente piccolo. evidente che la consistenza si necessaria alla convergenza per un metodo numerico. Sulla relazione generale tra tali propriet torneremo pi avanti.

Alcuni algoritmi di integrazione temporale Possiamo introdurre alcuni degli algoritmi di integrazione temporale

partendo dalla forma integrale del problema (2.2); difatti se applichiamo al generico intervallo ( )1,k kt t + la (2.3):

11( ) ( ) ( , ) ,k

k

t

k k tx t x t f x t dt++ = + (2.6)

approssimando lintegrale nella (2.6) in diversi modi otteniamo i relativi algoritmi.

Algoritmi di Eulero 31


Algoritmi di Eulero Considerando per la (2.6) lapprossimazione:

11( ) ( ) ( ) ( ) ,k

k

t

k k ktx t x t f d f t t ++ = (2.7)

definiamo algoritmo di Eulero esplicito (Forward Euler)la relazione:

1( ) ( ) ( ) ,k k kx t x t f t t+ = + (2.8)

(laggettivo esplicito si riferisce al fatto che per il calcolo del punto 1( )kx t + sufficiente conoscere i dati relativi al passo k-esimo). In Figura 2.3 viene data uninterpretazione geometrica dellapprossimazione fatta, che consiste nellapprossimare larea sottesa dalla funzione ( ) ( , )x t f x t= con il rettangolo di area ( )kf t t .

Figura 2.3: interpretazione grafica algoritmo di Eulero esplicito

Val la pena osservare che allalgoritmo definito dalla (2.8) si pu anche

pervenire a partire dallo sviluppo della soluzione ( )x t in serie di Taylor attorno al punto kt , per calcolare in modo approssimato il valore in 1kt + , otteniamo:

1 1( ) ( ) ( )( ) ....k k k k kx t x t x t t t+ + + + (2.9) che nel nostro caso, avendo assunto t costante, si semplifica in:



1( ) ( ) ( ) , con ( ) ( , ).k k kx t x t x t t x t f x t+ + = (2.10)

Considerando ora per la (2.6) lapprossimazione:

11 1( ) ( ) ( ) ( ) ,k

k

t

k k ktx t x t f d f t t ++ + = (2.11)

definiamo algoritmo di Eulero implicito (Backward Euler) la relazione:

1 1( ) ( ) ( ) ,k k kx t x t f t t+ += + (2.12)

(laggettivo implicito si riferisce al fatto che per il calcolo del punto 1( )kx t + ora necessario risolvere unequazione implicita per la presenza del termine

1( )kf t + ). In Figura 2.3 viene data uninterpretazione geometrica dellapprossimazione fatta, che consiste nellapprossimare larea sottesa dalla funzione ( ) ( , )x t f x t= con il rettangolo di area 1( )kf t t+ .

Figura 2.4: interpretazione grafica algoritmo di Eulero implicito

Anche in questo caso possibile interpretare la (2.12) attraverso un

opportuno sviluppo in serie di Taylor attorno al punto intorno a 1kt + , ottenendo:

1 1( ) ( ) ( ) , con ( ) ( , ),k k kx t x t x t t x t f x t+ + + = (2.13)

Sulle propriet e le differenze tra tali algoritmi ci soffermeremo pi avanti. Anticipiamo sin dora che entrambi risultano convergenti e hanno un errore dello stesso ordine di grandezza. Lalgoritmo esplicito risulta pi semplice da

Algoritmo dei trapezi 33


implementare e con minore complessit computazionale rispetto a quello implicito; vedremo per pi avanti che questultimo presenta altri vantaggi.

Algoritmo dei trapezi Abbiamo gi visto che la differenza 1( ) ( )k kx t x t+ pu essere interpretata

tramite un integrale come nelle equazioni (2.6). Se si approssima lintegrale con la cosiddetta regola dei trapezi:

[ ]11 1( ) ( ) ( ) ( ) ( )2k

k

t

k k k kt

tx t x t f d f t f t ++ +

= + . (2.14) si ottiene la formula iterativa:

[ ]1 1( ) ( ) ( ) ( )2k k k ktx t x t f t f t+ +

= + + . (2.15)

La formula dei trapezi approssima lintegrale con il trapezio definito dai

valori in kt e 1kt + della funzione integranda, come mostrato nella Figura 2.5.

Figura 2.5 interpretazione grafica algoritmo dei trapezi

Lalgoritmo definito dalla (2.15) detto dei trapezi (Trapezoidal) o di

Crank-Nicholson. Esso presenta un errore pi basso rispetto a quelli di Eulero. un algoritmo di tipo implicito a causa della presenza del termine

1( )kf t + a secondo membro della (2.15).



Algoritmo di Gear possibile definire unampia categoria di algoritmi attraverso formule

ricorsive alle differenze. In particolare considerata la relazione:

10 1

n n

k i k i i k ii i

x a x b x+ = =

= + . (2.16)

Essa definisce, al variare dellordine n e dei coefficienti ia e ib una classe di algoritmi che prende il nome di algoritmi di Gear. Al contrario degli algoritmi visti sinora per il calcolo di 1kx + non sono sufficienti i valori di

,k kx f ma vengono utilizzati un certo numero di passi precedenti. Sulla base di quanto appena osservato per lalgoritmo di Gear possiamo

classificare gli algoritmi di integrazione nelle seguenti classi: - algoritmi espliciti - algoritmi impliciti - algoritmi ad un passo - algoritmi a pi passi. possibile scrivere in forma generale le espressioni degli algoritmi

secondo la classificazione introdotta. Ad esempio, se consideriamo gli algoritmi espliciti ad un passo avremo:

1( ) ( ) ( , , ; )k k k k kx t x t t t x f t+ = + . (2.17)

Per gli algoritmi impliciti ad un passo avremo:

1 1( ) ( ) ( , , , ; )k k k k k kx t x t t t x f f t+ += + . (2.18)

Per gli algoritmi a pi passi avremo:

10

( , , , , ; ) + + =

= n

i k i k k k n k n ki

a x t t x f f f t . (2.19)

Algoritmo di Gear 35


La funzione ( ) che compare nelle precedenti espressioni detta funzione incremento. Le sue propriet si riflettono in modo importante sulle propriet degli algoritmi, come vedremo pi avanti.

Propriet di consistenza Vogliamo ora studiare in modo puntuale la propriet di consistenza di un

algoritmo per poi andarla a verificare in via generale. Riferendoci al solito al problema di Cauchy definito dalla (2.2) e considerata la sua soluzione ( )x t ricordiamo che un algoritmo detto convergente se accade che

00lim ( ) k kt x x t x = . (2.20)

Abbiamo gi accennato al fatto che condizione perch la (2.20) sia

verificata che lo schema numerico sia consistente, ovvero che considerata lespressione della formula iterativa essa risulti soddisfatta con sufficiente accuratezza quando al suo interno sostituiamo i valori della soluzione ( )x t . Consideriamo ad esempio lalgoritmo di Eulero esplicito:

1 ,k k kx x f t+ = + (2.21)

Se nella formula iterativa operiamo le sostituzioni 1 1( )k kx x t+ += ed ( )k kx x t= otteniamo:

1( ) ( ) ( ( )) ( ) ( )k k k k kx t x t t f x t x t t x t+ = + = + . (2.22)

immediato osservare che la (2.22) coincide con lo sviluppo in serie di Taylor della 1( )kx t + attorno a kt arrestato al primo ordine. Tenuto conto che per tale sviluppo vale la relazione:

21( ) ( ) ( ) ( )k k kx t x t t x t O t+ = + + , (2.23)

e definendo a partire dalla (2.21) un residuo 1kr + come: 1 1( )k k k kx t x f t r+ += + + , (2.24)



si ha:

10

lim 0kt

rt+

=

. (2.25)

Pertanto il metodo (in questo caso di Eulero esplicito) consistente al

problema di Cauchy di partenza. Si osservi che ci equivale a dire che lerrore locale tende a zero al tendere a zero di t . evidente che tale condizione sia necessaria per la convergenza del metodo.

Landamento del residuo pu in generale essere utilizzato per definire lordine del metodo. Difatti si dice che un metodo di ordine p se accade che:

1 1 10 0lim 0; lim , k k

p pt t

r r Ct t+ +

+ = =

, (2.26)

con C costante limitata.

possibile dimostrare che per un metodo convergente di ordine p (per il quale lerrore globale GTE di ordine pt ) lerrore locale LTE di ordine

1pt + . Il discorso sin qui fatto in relazione allalgoritmo di Eulero esplicito si pu

facilmente generalizzare alla classe degli algoritmi ad un passo espliciti. Difatti, considerata lespressione generale (2.17), possiamo ad essa applicare il ragionamento appena sviluppato nel caso di Eulero esplicito:

1( ) ( ) ( , ( ), ( ( ); )+ = + k k k k kx t x t t t x t f x t t , (2.27)

osservando che sicuramente dovr essere:

0lim ( , ( ), ( ( ); ) ( , ( ))

=k k k k kt t x t f x t t f t x t , (2.28)

tenuto conto anche qui della (2.23) e definito il residuo come nella (2.24) si ha, per la classe dei metodi in esame:

10

lim 0kt

rt+

=

. (2.29)



Zero-stabilit e teorema di Lax-Richtmyer

Vogliamo ora formulare una richesta di stabilit agli algoritmi che andiamo a considerare. Ci evidentemente una condizione affinch eventuali perturbazioni, anche eventualmente dovute al solo errore di round-off non vadano ad amplificarsi al crescere del numero di iterazioni pregiudicando la convergenza, e dunque il senso stesso dellutilizzo dello schema numerico.

Considereremo per illustrare il concetto un generico schema esplicito ad un passo:

1( ) ( ) ( , , ; )k k k k kx t x t t t x f t+ = + . (2.30)

e definiamo due generiche perturbazioni:

0 0

0 0

; ;

; .

= + = +

= + = + k

k

x k k

x k k

x x

x x (2.31)

Se accade che:

max 0 00, 0 , ,

:

>

k k

k k

t t x x

C x x C (2.32)

Allora il metodo si dice zero-stabile. Osserviamo che la condizione posta equivale a stabilire che:

00lim 0 (nell'ipotesi che per 0 si ha 0).

= k kt x x t x (2.33)

Questo il motivo della denominazione di zero-stabilit. Osserviamo inoltre che ci si riferisce ad una relazione relativa allalgoritmo, e non al problema di Cauchy sottostante.

Esistono diverse condizioni che garantiscono la zero-stabilit di un algoritmo. In particolare per i metodi ad un passo espliciti vale il seguente teorema:

se la funzione incremento Lipschitziana rispetto ad x allora il metodo zero-stabile.



abbastanza facile mostrare in modo diretto che il metodo di Eulero esplicito zero stabile. In modo piuttosto analogo si pu procedere anche per Eulero implicito. Anche i metodi dei Trapezi e di Gear risultano zero-stabili.

La propriet di zero-stabilit di un algoritmo assume importanza fondamentale in virt del seguente teorema di Lax-Richtmyer, che qui enunciamo in forma piuttosto sintetica e grossolana:

+

consistenzaconvergenza

zero stabilit

Limplicazione da sinistra a destra piuttosto scontata; la dimostrazione del teorema nelle condizioni pi generali invece complessa e viene omessa, rimandando alla letteratura specializzata.

Assoluta stabilit (caso lineare mono-dimensionale)) La propriet di zero-stabilit valuta landamento con 0 t e su un

intervallo limitato la stabilit dellalgoritmo. Poich gli intervalli di analisi possono essere illimitati e comunque arbitrariamente grandi necessario considerare anche la stabilit dellalgoritmo al crescere dellintervallo di analisi, ovvero con t fissato e k . Si parla in questo caso di stabilit o assoluta stabilit dellalgoritmo.

Analizziamo lassoluta stabilit degli algoritmi sin qui considerati anzitutto in riferimento al problema (2.2) nel caso lineare ed omogeneo (osserviamo che per i problemi lineari lanalisi di stabilit conduce naturalmente a problemi omogenei):

00

( >0) ( )(0)

tx x x t x ex x

=

= = (2.34)

In riferimento al problema (2.34) possibile esprimere in forma esplicita il

valore al passo n-esimo in funzione del valore iniziale. a) Eulero esplicito La formula iterativa, in questo caso diviene:



1 k k k kx x t x t x+ = = , (2.35)

ed andando ad esplicitarla a partire dal passo 0, si ottiene: 1 0(1 )x t x= , 2

2 1 0(1 ) (1 )x t x t x = = e dunque, in generale: 0(1 )

nnx t x= (2.36)

b) Eulero implicito La formula iterativa, in questo caso sar: 1 1 1 k k k kx x t x t x+ + + = = (2.37)

ed andando ad esplicitarla a partire dal passo 0, si ottiene: 11 0(1 )x t x= + ,

22 0(1 )x t x

= + , dunque: 0(1 )

nnx t x

= + (2.38) c) Trapezi In questo caso si ha, per la formula iterativa:

1 1 1( ) ( )2 2k k k k k kt tx x x x x x+ + +

= + = + (2.39)

che al solito, seguendo il procedimento gi utilizzato nei casi precedenti, conduce facilmente allespressione:

01 0,51 0,5

n

ntx xt

= +

(2.40)

Analogamente si pu procedere per lalgoritmo di Gear. La questione della stabilit assoluta viene studiata in relazione al problema

lineare ed omogeneo definito dalla (2.34), che pu esser visto come la differenza di due problemi lineari con lo stesso forzamento e condizioni iniziali distinte. In tal modo lo studio della stabilit ricondotto alla verifica della convergenza a zero degli algoritmi per n . Sappiamo, infatti che:



0( ) , 0tx t x e = > , (2.41)

da cui immediatamente discende lim ( ) 0

tx t

= . Vediamo cosa accade nei tre

casi in esame. Consideriamo dapprima il caso a) ovvero Eulero esplicito; tenuto conto dellespressione (2.36) si ha:

0lim lim(1 )

nnn n

x t x

= (2.42)

dunque la condizione per la quale lim 0nn x = sar data da 1 1t < , e

considerando che 0t > , si ottiene facilmente: 2 2t < = . (2.43) La (2.43) esprime la condizione cosiddetta di stabilit per lalgoritmo di

Eulero esplicito: perch la soluzione approssimata converga al risultato esatto in questo caso necessario che il passo di discretizzazione t scelto sia in una prefissata relazione con la frequenza naturale (o con la costante di tempo ).

Vediamo ora cosa accade per il caso b), ovvero Eulero implicito. Tenuto conto dellespressione (2.38) per nx , avremo:

0lim lim(1 ) 0 ( 0)

nnn n

x t x t

= + = > . (2.44)

In questo caso si ha che lalgoritmo come usa dirsi intrinsecamente (o

incondizionatamente) stabile. Infine, analizzando il caso c) (ovvero lalgoritmo dei trapezi), ci si rende

conto facilmente che:

01 0,5lim lim 0 ( 0).1 0,5

n

nn n

tx x tt

= = > +

(2.45)



Assoluta stabilit (caso non lineare mono-dimensionale) Proviamo ora ad analizzare la questione dellassoluta stabilit pi in generale e dunque in riferimento al caso non lineare: Pertanto consideriamo il sistema:

0 0

( )( )

x f xx t x=

= (2.46)

Iniziamo con lo studiare la stabilit dellalgoritmo di Eulero esplicito.

Consideriamo due sequenze che differiscono per le sole condizioni iniziali:

1 11 1

' ( ' ) ''' ( '' ) ''n n n

n n n

x f x t xx f x t x

= += +

(2.47)

Il metodo risulter stabile se sar verificata la successione di disuguaglianze: 1 1 0 0 2 2 1 1' '' ' '' ; ' '' ' '' ......x x x x x x x x (2.48) Considerata la disuguaglianza n-esima, sostituendo le (2.47) abbiamo: [ ]1 1 1 1 1 1( ' ) ( '' ) ' '' ' '' + n n n n n nt f x f x x x x x (2.49)

Utilizzando il teorema di Lagrange possiamo scrivere:

1 1 1 1( ) ( ' ) ( '' ) ( ' '' )n n n n

x

df xf x f x x xdx =

=

, (2.50)

dove 1 1'n nx x . Sostituendo la (2.50) nella (2.49) abbiamo dunque:

( )1 1 1 1( ) 1 ' '' ' ''n n n ndt f x x x xdx

+

(2.51)

ovvero:

( ) 1 1 + dt fdx

(2.52)



che, tenuto conto dei valori assoluti conduce alle condizioni:

2( ) 2 ( )

( ) 0 ( ) 0

d dt f fdx dx td dt f fdx dx

(2.53)

ed infine:

2( )

x

tdf x

dx =

(2.54)

Siccome durante la dinamica non sappiamo a priori come varia , per t fisso la condizione di stabilit in generale data se si considera la massima

derivata, max

( )

df xdx

.

Vediamo ora cosa accade nel caso di Eulero implicito. La stabilit si analizza analogamente al caso precedente, a partire due sequenze:

11

' ( ' ) ''' ( '' ) ''

= += +

n n n

n n n

x f x t xx f x t x

(2.55)

Con passaggi analoghi a quanto gi visto si ottiene in questo caso:

1 ( )1; 1 1( )1 x

x

df xtdxdf xt

dx

=

=

(2.56)

dunque:

( ) 2

( ) 0

x

x

df xtdx

df xtdx

=

=

(2.57)

Infine, nel caso del metodo dei trapezi, a partire dalla formula iterativa



1 1( ) ( )2 2n n n nt tx f x x f x

= + , (2.58)

con analoghi passaggi si ottiene la condizione di stabilit:

( )12 1

( )12

t dfdx

t dfdx

+

(2.59)

Essa risulta sempre verificata se ( ) 0df xdx

< , mentre se ( ) 0df xdx

> non mai

verificata.

Convergenza e stima dellerrore numerico Vogliamo ora valutare, sempre in riferimento agli algoritmi precedenti,

lerrore introdotto dalla soluzione numerica. Faremo ci anzitutto con un esempio molto semplice. Consideriamo allora, in riferimento sempre al problema espresso dalle eq. (2.34), il caso 1 = ; consideriamo alcuni istanti, ad esempio t=1, t=5, t=10, ed andiamo a calcolare nei vari casi le soluzioni.

1) 0.2, 50t n = =

t 1 5 10 Analitica 0.36788 0.0067379 4.54e-005 EulForw. 0.32768 0.0037779 1,43E-01 EulBack. 0.40188 0.010483 0.00010988 Trapezi 0.36665 0.0066259 4,39E-01

2) 0.02, 500t n = =

t 1 5 10 Analitica 0.36788 0.0067379 4.54e-005 EulForw. 0.36417 0.006405 4,10E-01 EulBack. 0.37153 0.0070788 5,01E-01 Trapezi 0.36787 0.0067368 4,54E-01



E immediato verificare che, nel caso degli algoritmi di Eulero, al passare da 0.2t = a 0.02t = (diminuzione di un ordine di grandezza, con un corrispondente aumento di iterazioni) si ottiene una diminuzione dello stesso ordine dellerrore. Ci corrisponde al fatto essi sono algoritmi del primo ordine. Per lalgoritmo dei trapezi, invece, la diminuzione di un ordine di grandezza del t comporta una diminuzione di due ordini di grandezza dellerrore! I risultati sono riassunti nella seguente tabella:

Eulero Trapezi

t=0.1. t =0.01 t =0.1. t =0.01 2a cifra 3a cifra 4a cifra 6a cifra 3a cifra 4a cifra 5a cifra 7a cifra 5a cifra 6a cifra 7a cifra 9a cifra

Figura 2.6

Convergenza ed errore di troncamento per Eulero esplicito Dopo aver mostrato con un esempio numerico quel che accade, vogliamo

ora valutare o almeno stimare in modo analitico lerrore che si commette rispetto alla soluzione esatta, ci che ci consentir di verificare al tempo stesso la convergenza del metodo. Riferendoci ancora una volta ad un caso lineare, questa volta non omogeneo (con forzamento 0b costante), avremo:

00 0( )

x x bx t x

= + =

(2.60)

Convergenza ed errore di troncamento per Eulero esplicito 45


cui corrisponde la soluzione:

0 00( )tb bx t x e

= +

. (2.61)

La formule iterativa per Eulero esplicito in questo caso diviene:

1 0(1 )k kx t x b t+ = + (2.62)

Tenuto conto che 0x x b= + , la successione dei valori fornita dallalgoritmo :

1 0 0

22 1 0 0 0 0

1

0 00

(1 )

(1 ) (1 ) (1 )

(1 ) (1 )n

n in

i

x t x b t

x t x b t t x t b t b t

x t x b t t

=

= +

= + = + +

= +

(2.63)

Consideriamo anche, parallelamente alla successione definita

dallalgoritmo, quella definita negli istanti n t mediante lespressione dello sviluppo in serie con il resto di Lagrange:



( )

( )

21

20 1

22

20 0

2 21 2

1( ) (0) (0)2

1 (1 ) (0) ( )2

1(2 ) ( ) ( )2

(1 ) (0) (1 )1 1 (1 ) ( ) ( )2 2

( ) (1

x t x x t x t

t x b t x t

x t x t x t t x t

t x t b t b t

t x t x t

x n t

= + + =

= + +

= + + =

= + + +

+ +

=1

0 00

2

1

) (1 )

( )(1 )2

nn i

i

nn i

ii

t x b t t

t x t

=

=

+ +

+

(2.64)

Confrontando ora le espressioni di x2 e (2 )x t osserviamo che esse

differiscono per i termini che contengono le derivate seconde; in particolare 2( )x lerrore introdotto dalla seconda iterazione, mentre il termine

contenente 1( )x conserva la memoria dellerrore dovuto alliterazione precedente.

Possiamo in generale definire lerrore algoritmico (o di troncamento) ( )nte al passo n come la differenza tra i termini n-esimi delle due successioni:

2

( )

1( ) ( )(1 ) .

2

=

= =

nn n i

t n ii

te x n t x x t (2.65)

utile riscrivere la (2.65) liberandosi dalle derivate; ci pu essere fatto

ricordando che 02 bxxx == . In tal caso lespressione di ( )nte

diviene:

( )( )( )

( )( ) ( )

2( ) 2

01

22 20

1 1

12

1 12 2

=

= =

= =

nn in

t ii

n nn i n i

ii i

te x b t

b tt x t t (2.66)

Convergenza ed errore di troncamento per Eulero esplicito 47


Infine, ponendo 1 t = possiamo semplificare ulteriormente lespressione (2.66). Difatti si riconosce facilmente che:

( ) 2 11

11 1 .... ,1

nnn i n

it

=

= + + + =

(2.67)

dunque sostituendo nellespressione dellerrore (2.66) si ha:

( )22 2

( ) 0

1

1 .2 2 1

=

=

nn

n n it i

i

b tte x (2.68)

Linterpretazione della (2.68) interessante: essa pu essere letta come la

somma di un errore transitorio e di un errore di regime. Il secondo termine semplice da studiare, in funzione di e t . Il primo termine, a causa della presenza di )( ix (che non noto) non pu essere valutato esattamente. Considerato il modulo dellerrore di troncamento ( )nte possibile effettuare

alcune maggiorazioni per studiare comunque landamento asintotico dellerrore. In particolare si ha:

22 2

0( )

1

2 20

max

1( )( )2 2 1

1 ,2 1

nnn n i

t ii

n

b tte x

bt x

=

+

+

(2.69)

dove la prima disuguaglianza ottenuta maggiorando la (2.68) con la somma dei moduli, e la seconda considerando maxx in luogo di ( )ix nella sommatoria e poi mettendo in evidenza tale valore.

Se a questo punto ricordiamo che per la condizione di stabilit per Eulero esplicito deve essere 2 /t < , e ci implica 0 1 , possiamo verificare la convergenza dellalgoritmo effettuando il limite per n della espressione trovata nella (2.69) e che maggiora il modulo dellerrore:



2 2max 0 max 0

2max 0 max 0

1 1 1 1lim ( ) ( )2 1 2 1

1 1( ) ( ).2 2

n

nt x b t x b

tt x b x bt

+ = + =

= + = +

(2.70)

Lalgoritmo dunque convergente e lerrore di troncamento va a zero con legge lineare in funzione di t (metodo del primo ordine).

Errore di round off per Eulero esplicito Abbiamo sin qui valutato lerrore che intrinsecamente si commette

nellutilizzare lalgoritmo, errore che abbiamo chiamato errore di troncamento o algoritmico. Vogliamo ora in qualche modo valutare lerrore di round off. Per semplicit facciamo riferimento al caso lineare omogeneo:

0 0

,( ) .

x xx t x

= =

(2.71)

Considerando sempre lalgoritmo di Eulero esplicito e denominiamo (al

solito) 1,.... nx x la successione dei valori con esso calcolati, immaginando di utilizzare unaritmetica a precisione illimitata. Come visto in precedenza la successione in questo caso:

1 0

22 1 0

0

(1 )

(1 ) (1 )

(1 )nn

x t x

x t x t x

x t x

=

= =

=

(2.72)

Il risultato fornito da una qualsiasi implementazione reale, con precisione

questa volta limitata, sar affetto dallerrore di round off. Possiamo indicare con 1,.... nx x , per distinguerla dalla prima, la sequenza ottenuta in questo secondo caso. Potremo allora esprimere:

1 0 1(1 )x t x x= + (2.73)

Errore di round off per Eulero esplicito 49


dove 1x proprio larrotondamento dovuto allaritmetica a precisione limitata. Iterando, anche in questo caso, avremo poi:

22 1 2 0 1 2(1 ) (1 ) (1 )x t x x t x t x x = + = + + (2.74)

ed, in generale

01

(1 ) (1 )n

n n in i

ix t x t x

=

= + (2.75) Lerrore di arrotondamento (round-off) ( )nre al passo n ha dunque

lespressione:

( )1

(1 ) =

= = n

n n ir n n i

ie x x t x (2.76)

Anche in questo caso possiamo effettuare una maggiorazione:

max1 1

1(1 ) (1 )1

nn nn i n i

i ii i

t x t x x

= =

(2.77)

Facendo il solito limite n otteniamo

max max( )lim1

= =

nrn

x xe

t (2.78)

Lespressione (2.78) molto significativa: difatti essa mostra come lerrore di round off diverga per 0t . Tale circostanza fa capire chiaramente come, quando si operi con aritmetica a precisione limitata, anche senza considerare lonere computazionale, non si possa pensare di far diminuire arbitrariamente il passo di discretizzazione t .

Mettendo ora insieme i risultati ottenuti a proposito dei due tipi di errore considerati, possiamo dunque tracciare un grafico, se pur qualitativo, dellandamento dellerrore complessivo in funzione del passo di integrazione t .



prevale lerrore di round off

t

min

t ottimo

prevale lerrore di troncamento

Figura 2.7: andamento dellerrore frobale in funzione di t per Eulero esplicito

E immediato osservare che esiste una limitazione intrinseca, se pur dipendente dalla lunghezza della parola macchina con cui i numeri vengono rappresentati, alla precisione che possiamo sperare di ottenere diminuendo lampezza dello step t .

Algoritmi Runge Kutta Al termine di questa sezione vogliamo accennare ad una importante

categoria di algoritmi che sono prendono il nome di Runge-Kutta. possibile introdurre lidea di base a partire dalle seguenti considerazioni. Consideriamo la formula di Eulero esplicito (2.8) che riportiamo per comodit:

1 ( , ).+ = + k k k kx x t f x t (2.79)

Sappiamo che la stima (locale) fornita dalla (2.79) accurata nellordine 2( )O t , e ci pu essere visto come effetto del fatto che la ( )f valutata in

un solo estremo dellintervallo 1+ = k kt t t (in questo caso kt ). Supponiamo ora di utilizzare la (2.79) come prima stima dellincremento 1 1+= k kk x x ; a questo punto applichiamo nuovamente la (2.79) dove questa volta per la

( )f viene valutata in corrispondenza del punto 1,2 2

+ +k kktt x :



11 ( , ).2 2+

= + + +k k k kk tx x t f x t (2.80)

possibile mostrare che in tal modo si ottiene una stima pi accurata rispetto a quella data da Eulero esplicito. In particolare si ha:

1

12

31 2

( , ),

, ,2 2

( ).+

=

= + +

= + +

k k

k k

k k

k t f x tk tk t f x t

x x k O t

(2.81)

Questo un primo esempio di un metodo Runge-Kutta (a due stadi). Un altro esempio rappresentato dal metodo di Heun, che si pu ottenere a partire da quello dei trapezi utilizzando unapprossimazione di 1 1( , )+ +k kf x t ottenuta utilizzando Eulero esplicito:

[ ]1 ( , ) ( ( , )) .2k k k k k k ktx x f x t f x t f x t+

= + + + (2.82)

Esso prende anche il nome di metodo dei trapezi modificato. Osserviamo che la (2.82) produce un metodo esplicito al contrario della formula dei trapezi. Entrambi gli schemi (2.81) e (2.82) possono essere scritti come:

11

.+=

=Ns

k k s ss

x x c k (2.83)

dove Ns indica il numero di stadi del metodo, che poi corrisponde al numero di valutazione della funzione ( )f .

qui opportuno sottolineare alcune caratteristiche e propriet generali dei metodi RK:

- i metodi RK sono ancora ad un passo, ma richiedono pi valutazioni della ( )f allinterno del singolo passo. Ci rende pi facile, rispetto ad algoritmi multistep, implementare il passo adattativo

- per quanto riguarda le propriet di convergenza valgono quelle generali per i metodi ad un passo

- possibile dimostrare che un metodo RK (esplicito) ad Ns stadi non pu avere ordine maggiore di 4Ns .



A causa di questultima propriet molto diffuso lo schema RK del 4 ordine:

( )

1

12

22

1 3 1

51 1 2 3 4

( , ), , ,2 2 , ,2 2

( , ),

2 2 ( ).6

+

+

=

= + +

= + +

= +

= + + + + +

k k

k k

k k

k k

k k

k f x tt k tk f x t

t k tk f x t

k f x t k ttx x k k k k O t

(2.84)

possibile dimostrare che tale procedura generalmente conveniente rispetto alle (2.81) o (2.82) se a parit di accuratezza si riesce a mantenere un passo t doppio. Tale condizione risulta generalmente verificata dal punto di vista pratico.


3. Integrazione di equazioni ODE: caso n-dimensionale

Abbiamo sin qui trattato la soluzione numerica di circuiti dinamici del primo ordine. Vogliamo ora considerare la generalizzazione di quanto visto al caso dei sistemi di equazioni differenziali ordinarie, cui perveniamo quando consideriamo circuiti dinamici di ordine superiore al primo. Per la valutazione della complessit e dellonere computazionale ha senso riferirsi al caso lineare; ricordiamo che un circuito dinamico lineare di ordine N descritto da equazioni nella forma: C D (t)= +x x f ovvero -1 -1C D C (t) A (t)= + = +x x f x b , (3.1) dove la matrice C, in generale non diagonale se vi sono accoppiamenti mutui. In tal caso lespressione dellalgoritmo di Eulero esplicito diviene: 1 ( )n n n n n ntA t I tA t+ = + + = + + x x x b x b (3.2) Limplementazione della (3.2), anche nel caso n-dimensionale, molto agevole in quanto non si deve effettuare alcuna inversione di matrice o soluzione di sistema lineare. Ci ha un ovvio riflesso anche sullonere computazionale, che facilmente valutabile per la singola iterazione. Per lalgoritmo di Eulero implicito, la formula iterativa assume la forma: ( )-11 1( )+ += + n n nI tA tx x b (3.3) Per lalgoritmo dei trapezi, invece si avr: -11 1( 0.5 ) ( 0.5 )( )+ += + + +n n n nI tA I tAx x b b (3.4)

54 Integrazione di equazioni ODE: caso n-dimensionale


In entrambi i casi necessario realizzare una inversione di matrice (o equivalentemente la soluzione di un sistema lineare) ad ogni passo di iterazione.

Diagonalizzazione delle equazioni dinamiche Vogliamo ora analizzare la questione della stabilit dei metodi considerati nel caso di circuito lineare di ordine qualsiasi, ovvero dellequazione: ( ) ( )= +A t tx x b (3.5) Al solito, ci vuol dire di considerare due soluzioni vicine che differiscono per la condizione iniziale, e verificare se la loro distanza tende a diminuire. Per far ci utile considerare allora, lequazione per le differenze: date le due soluzioni '' e ' xx consideriamo ''' xxx = ; in tal caso il problema diviene:

( )0 0 =

=

At

x xx x

(3.6)

Consideriamo, per la matrice A che supponiamo di ordine n, linsieme dei suoi autovalori, ovvero delle soluzioni dellequazione 0A I = , cui perveniamo studiando il problema: A i i i=u u . (3.7) con lipotesi di soluzione diversa da quella banale. Fissato un autovalore i il corrispondente vettore iu si dice un autovetture per il sistema. Ad autovalori distinti corrispondono autovettori indipendenti. In tale ipotesi la matrice: 1 2( , ,..... )nT = u u u , (3.8) ha rango pieno e dunque invertibile. Dalla (3.7) si ha anche che: 1 2 1 1 2 2( , ,.... ) ( , ,..... )n n nAT A = =u u u u u u , (3.9) Introducendo la matrice diagonale 1diag( .... )n = la (3.9) pu scriversi:



1 1 2 2 n( , ,..... )nAT T = = u u u (3.10) da cui si ricava facilmente che: -1A T T= (3.11) A questo punto il sistema (3.6) pu scriversi come 1T T = x x (3.12) Se poniamo ora -1T= y x e dunque anche -1y T= x , avremo = y y . (3.13) Con tale cambiamento di variabili (base) il sistema risulta in forma diagonale, pur mantenendo gli stessi autovalori. Dunque in questa forma il sistema costituito da n equazioni disaccoppiate, per cui possibile considerare ciascuna equazione separatamente, ovvero: i i iy y= ( ni .....1= ) (3.14)

Criteri di assoluta stabilit nel caso n-dimensionale Il cambiamento di base appena considerato a questo punto di fondamentale importanza al fine di studiare la stabilit. Difatti ora il problema potr essere ricondotto ad n problemi scalari come quelli affrontati in precedenza. La condizione di stabilit pu essere espressa in generale come ( 1) ( )+ k ky y

mentre per lasintotica stabilit dovr verificarsi ( 1) ( )+

Eulero implicito 57


Pu pertanto accadere che soluzioni stabili nella realt risultino instabili nellalgoritmo se t non sufficientemente piccolo. Invece non potr mai accadere che soluzioni intrinsecamente instabili ( 0) > vengono stabilizzate numericamente, essendo la circonferenza sempre nel semipiano

0 , qualunque sia il t considerato.

Eulero implicito Ricordiamo che la condizione di stabilit per il caso di Eulero implicito espressa dalla 1 1t , e con lo stesso procedimento considerato prima, facilmente si perviene alla condizione: ( )2 2 21 1t t + : (3.17) Anche in questo caso possiamo dare uninterpretazione geometrica della situazione. Dovremo in tal caso considerare la circonferenza di equazione:

2

22

1 1t t

+ = , (3.18)

ed i punti di stabilit per lalgoritmo saranno dati, questa volta, da quelli esterni alla circonferenza, come mostrato in Figura 3.2.

t

regione di stabilit

Figura 3.2: regione di stabilit assoluta per Eulero implicito



Osserviamo che, in questo caso, per i modi stabili ( )0 < si ottengono sequenze stabili qualsiasi sia la scelta di t . Pu per capitare che anche nel caso di modi instabili ( )0 > si ottengono sequenze stabili. In tal caso solo riducendo t si pu sperare di includere i corrispondenti autovalori nella regione di instabilit.

Trapezi La condizione di stabilit espressa, in tal caso da

2 2 2

2 2 2

1 0.5 (1 0.5 ) 0.251 11 0.5 (1 0.5 ) 0.25

t t tt t t

+ + +

+ (3.19)

ovvero:

( )( )

22

22

21

2t

t

+ + +

(3.20)

Essa suscettibile dellinterpretazione grafica di Figura 3.3.

Figura 3.3: regione di stabilit assoluta per il metodo dei Trapezi

Algoritmo di GEAR 59


Algoritmo di GEAR Abbiamo gi introdotto lalgoritmo di Gear in un paragrafo precedente. Riportiamo qui per comodit lespressione della formula iterativa (2.16):

10 1

n n

k i k i i k ii i

x a x b x+ = =

= + . (3.21) E possibile mostrare in tal caso che, in relazione allordine n dellalgoritmo, si hanno le figure di stabilit riportate in Figura 3.4.

Figura 3.4: regioni di stabilit per lalgoritmo di GEAR



Caso non lineare n-dimensionale Consideriamo infine il problema:

0 0

( )( )t=

=

x f xx x

(3.22)

Andiamo anzitutto a verificare la stabilit degli algoritmi. Partiamo dal caso di Eulero esplicito: 1 ( )n n nt = +x f x x (3.23) Consideriamo al solito due sequenze che differiscono per la condizione iniziale: [ ]1 1' '' ( ' ) ( '' ) ' ''n n n nt+ + = + n nx x f x f x x x . (3.24) Anche in questo caso, per studiare la stabilit, si applica lespressione:

( ' ) ( '' ) ( ' '' )n n n ndd =

= x

ff x f x x xx

(3.25)

Dunque lequazione alle differenze pu scriversi come:

n

x1n xx

xfx

+=

=+ I

)(

ddt

. (3.26) Si avr in tal caso che 1n n+ x x se gli autovalori della matrice

Jacobiana appartengono al cerchio di raggio 1t

e centro 1 ,0t

(come nel

caso lineare). Nel caso invece dellalgoritmo di Eulero implicito: 1 1( )n n nt x+ + =x f x , (3.27)

Caso non lineare n-dimensionale 61


la sequenza alle differenze verifica lequazione:

1( )- n n

dI td +=

=

x

f x x xx

. (3.28)

Essa risulta stabile se gli autovalori di della matrice Jacobiana non cadono nel cerchio di raggio 1 t e centro ( )1 ,0t . Infine, per quanto riguarda il metodo dei trapezi:

1( ) ( )2 2n nt tf x f x+ +

= +n 1 nx x (3.29)

si ottiene

1I- 2 2n nt d t dI x

d d +

= +

f fxx x

(3.30)

Essa stabile se gli autovalori della matrice Jacobiana (considerando la linearizzazione ad ogni passo) hanno sempre parte reale negativa.

Cenni a circuiti stiff ed algoritmi a passo variabile Abbiamo visto che la stabilit di un algoritmo pone in genere delle condizioni sullo step temporale t . A parte i problemi di onere computazionale legati alla riduzione dello step, sappiamo che lerrore di arrotondamento diverge per

0t . Dunque, per diversi motivi da un punto di vista realizzativo vi linteresse allutilizzo di un t quanto pi grande possibile. In realt, il problema se possibile ancor pi complesso. Difatti, tenuto conto di dinamiche che possono avere tassi di variabilit assai differenti in successive regioni temporali, il problema che si pone quello di poter avere addirittura un time step t variabile che possa seguire le esigenze di stabilit e precisione lungo la dinamica. Il caso classico nel quale si pone in modo stringente tale esigenza quello di circuiti stiff, ovvero circuiti che presentano costanti di tempo assai differenti le une dalle altre.



Consideriamo allora il circuito rappresentato in Figura 3.5. Per esso immediato ricavare lequazione di stato (si tratta ovviamente di un circuito del primo ordine):

Figura 3.5

[ ]1RC

sc c s

dvv v vdt

= + (3.31)

ovvero

[ ]1( )( ) ds tx x s t

dt= + (3.32)

La soluzione del problema data da 1( ) ( )tx t ke s t= + , con ( ) (0)k x o s= (3.33) Consideriamo allora 2( ) 1 ts t e = 1 2( ) ( ) 1t tx t x o e e = + , (3.34) e supponiamo che 1 e 2 differiscano per svariati ordini di grandezza. La soluzione del problema in tal caso quella rappresentata in Figura 3.6. Ebbene, quando nella soluzione esiste una componente veloce ed una lenta, il problema si dice STIFF. Supponiamo, ad esempio,

1 e 10 29

1 == ; In tal caso la prima componente decade (ed dunque trascurabile) per 95 10t s mentre la seconda raggiunge lo stato di regime per 5t s .

Caso non lineare n-dimensionale 63


Figura 3.6

Immaginiamo, a questo punto, di voler utilizzare lalgoritmo di Eulero forward: sappiamo che la condizione di stabilit data da:

9102t (3.35)

e con questo time step dovremmo risolvere un intervallo di 5s! E immediatamente evidente come ci richieda un numero di passi elevatissimo, in pratica inaccettabile. Daltro canto, se proviamo a violare la condizione (3.34) considerata, facile mostrare che si ottiene una sequenza instabile, come mostrato in Figura 3.7

Figura 3.7

Sulla base di quanto visto a proposito della diagonalizzazione delle equazioni circuitali, si intuisce facilmente come considerazioni analoghe possano valere nel caso di circuiti lineari di ordine n con autovalori molto diversi fra loro. Infine, a proposito di circuiti dinamici non lineari, diremo che ci troviamo in



un caso stiff se la matrice Jacobiana, nei punti di lavoro, ha gli autovalori che differiscono molto fra loro. La soluzione ai problemi proposti lutilizzo di algoritmi con passo adattativo, che siano in grado di valutare dinamicamente la necessaria ampiezza del time step tale da garantire la stabilit e laccuratezza desiderata con il minimo onere computazionale possibile.

SOMMARIO PPUNTI DI SOLUZIONE NUMERICA DI CIRCUITI A ...unina.stidue.net/Teoria dei...

Documents

Transcript of SOMMARIO PPUNTI DI SOLUZIONE NUMERICA DI CIRCUITI A ...unina.stidue.net/Teoria dei...