5. Equazioni algebriche non lineari - Elettrotecnica · Solo in qualche caso è possibile risolvere...

G. Miano, Appunti del Corso di Modelli Numerici per i Campi, 2009

5. Equazioni algebriche non lineari

La soluzione dei sistemi di equazioni algebriche non lineari riveste un ruolo fondamentale nella modellistica numerica. Ad esempio, i sistemi a parametri concentrati non lineari in condizioni stazionarie sono governati da equazioni algebriche non lineari. Più in generale, la soluzione numerica di equazioni differenziali ed integrali non lineari si riduce alla soluzione di sistemi di equazioni algebriche non lineari.

Il problema che considereremo in questo capitolo è il seguente. Data una funzione f = f x( ) definita in D !"

n e avente codominio C !"n ,

f x( ) ! f

1x1, x

2,..., xn( ), f2 x

1, x

2,..., xn( ),..., fn x

1, x

2,..., xn( )

T

"C # $n , (5.1)

bisogna cercare un vettore

u = u1,u

2,...,u

n

T

!D "#n (5.2)

tale che

f u( ) = 0!

f1u

1,u

2,...,un( ) = 0,

f2u

1,u

2,...,un( ) = 0,

...

fn u1,u

2,...,un( ) = 0.

"

#

$$

%

$$

(5.3)

Il problema (5.3) è un sistema di

n equazioni algebriche non lineari in

n incognite 1. Le soluzioni u sono dette anche “zeri” della funzione f . In queste Lezioni non affronteremo 1 In queste Lezioni non considereremo il problema più generale in cui il numero di equazioni (vincoli) è diverso dal numero di incognite (gradi di libertà).

2

G. Miano, Appunti di Modelli Numerici per i Campi, 2009

il problema più generale della soluzione di equazioni algebriche non lineari con soluzioni complesse. Un’equazione del tipo (5.3) può avere una sola soluzione, più di una soluzione o addirittura nessuna soluzione. Si intuisce che la soluzione dell’equazione (5.3) esiste ed è unica se:

a) il codominio di f contiene lo zero; b) la funzione f definisce un omeomorfismo2 del dominio D nel codominio C .

Lasciamo al lettore la verifica di questo risultato. In questa Nota illustreremo alcune tecniche per calcolare numericamente le soluzioni di un sistema di equazioni algebriche non lineari. In generale, i metodi di soluzione di un’equazione sono classificati in metodi diretti e metodi iterativi. I metodi diretti, in assenza degli errori di arrotondamento3, danno la soluzione esatta in un numero finito di operazioni. Nei metodi iterativi, invece, la soluzione è ottenuta come limite di una successione, in generale, infinita di problemi più semplici risolvibili attraverso metodi diretti. Solo in qualche caso è possibile risolvere un sistema d’equazioni algebriche non lineari con un metodo diretto4. I metodi numerici per la soluzione del problema (5.3) sono, in generale, di tipo iterativo. Prima affronteremo il problema della soluzione di un’equazione algebrica non lineare in una sola incognita (equazioni scalari), poi estenderemo i risultati ad un sistema di n equazioni in n incognite (equazioni vettoriali). 5.1 Equazioni scalari

Consideriamo l’equazione scalare f u( ) = 0 , (5.4)

dove

f = f x( ) è una funzione ad un solo valore definita in D ! " , continua e differenziabile. Se il codominio C di f contiene lo zero, l’equazione (5.4) ha almeno una soluzione; inoltre, se f è monotona in D allora la soluzione è unica. In questo caso la funzione f definisce un omeomorfismo del dominio D nel codominio C . La funzione f è monotona se !f

5 non cambia segno in D , ovvero se !f x( ) " 0 in D . Ciò

2 Una funzione f :D ! R

n" C ! R

n definisce un omeomorfismo di D in C se la corrispondenza definita da f tra D e C è biunivoca, f e l’inversa f !1 sono continue in D e C , rispettivamente. 3 Gli errori di arrotondamento sono quelli introdotti dal fatto che le macchine per il calcolo numerico possono rapprersentare solo numeri reali con un numero finito di cifre. 4 La soluzione analitica esiste per le equazioni algebriche scalari di secondo, terzo e quarto grado. 5 Con il simbolo !f indichiamo la derivata prima di f rispetto alla variabile x .

3


equivale a richiedere che !f x( )"1

# $ < % in D . Questo risultato è un caso particolare

del teorema di Hadamard, [5.1]. Lasciamo la lettore la verifica di questa proprietà. Se f non è monotona la (5.4) può avere più di una soluzione. Esempio 5.1

Si consideri un circuito composto da un generatore ideale di tensione E , un resistore

lineare di resistenza

R e un resistore non lineare

N controllato in tensione 6, Figura 5.1. L’equazione caratteristica di

N può essere così rappresentata i = g v( ) , (5.5)

dove

i e

v sono rispettivamente l’intensità di corrente e la tensione del bipolo (con i versi di riferimento concordi con la convenzione dell’utilizzatore).

Il resistore

N potrebbe essere, ad esempio, un diodo a giunzione pn o un diodo tunnel, Figura 5.2. Per il diodo a giunzione pn una buona approssimazione della funzione

g è data da

g v( ) = IS ev /VT !1( ) (5.6)

per tensioni molto più grandi della tensione di breakdown;

IS è l’intensità di corrente di

saturazione inversa e

VT è la tensione termica. Valori tipici di

IS e

VT sono,

rispettivamente,

10!15A e

26 mV . Per il diodo tunnel una buona approssimazione della funzione

g è data da

g v( ) = a1v + a

2v2+ a

3v3 (5.7)

per v ! 0 , dove

a1

= 6.0 !"1,

a2

= !4.5 "V( )!1 e

a3

= 1.0 !V2( )

"1

Figura 5.1 Un circuito non lineare

6 Un resistore si dice che è controllato in tensione (coorente) se l’intensità di corrente (tensione) è univocamente determinata per ogni valore di tensione (intensità di corrente) .

4


Figura 5.2 Curve caratteristiche del diodo a giunzione pn (a) e diodo tunnel (b).

L’equazione caratteristica del bipolo composto dalla serie generatore di tensione –

resistore è

i =E ! v

R. (5.8)

Combinando le equazioni (5.5) e (5.8) si ottiene l’equazione per la tensione del resistore non lineare

g v( ) +v

R!E

R= 0 . (5.9)

Ecco un esempio d’equazione algebrica non lineare scalare. La soluzione V dipende dai parametri R ed E , V =V R,E( ) .

Figura 5.3 La retta di carico definita dalla (5.8) interseca sempre in un solo punto la curva caratteristica del diodo a giunzione pn (a); la retta di carica può intersecare in tre punti distinti la curva caratteristica del diodo tunnel (b).

♦

5


Se il resistore non lineare è un diodo a giunzione pn l’equazione (5.9) ha una ed una sola soluzione per qualsiasi valore dei parametri; invece, nel caso di un diodo tunnel l’equazione (5.9) può avere una sola soluzione o tre soluzioni distinte, a seconda dei valori dei parametri. Ciò può essere facilmente verificato attraverso il metodo grafico illustrato in Figura 5.3.

♦ Assumiamo, per ora, che l’equazione (5.4) abbia una e una sola soluzione

nell’intervallo a,b[ ] di ! . A partire da un valore iniziale x0 opportunamente scelto, un

metodo iterativo per risolvere la (5.4) genera, secondo una legge opportuna x

k= G x

k!1( ) , (5.10) una successione di valori

x1, x

2,..., x

k,... (5.11)

che, in generale, converge alla soluzione

u per

k!" . Non essendo possibile iterare all’infinito la (5.10) attraverso una macchina di calcolo, occorre definire un criterio d’arresto per stabilire a quale iterazione N si ha una stima accettabile della soluzione,

u ! xN

. (5.12) Il problema del criterio d’arresto sarà trattato ampiamente più avanti per ognuno dei metodi iterativi presentati.

Esistono diversi metodi iterativi per risolvere un’equazione algebrica non lineare. In queste Lezioni illustreremo le caratteristiche principali di tre metodi:

• il metodo di bisezione; • il metodo di Newton-Raphson; • il metodo di Picard.

Il metodo di Newton-Raphson e il metodo di Picard si estendono senza alcuna difficoltà alle equazioni vettoriali, mentre l’estensione del metodo della bisezione presenta non poche difficoltà. 5.1.1 Metodo della bisezione

Se la funzione f x( ) verifica la condizione

f a( ) f b( ) < 0 , (5.13)

l’equazione (5.4) ha almeno una soluzione nell’intervallo a,b[ ] . Si assuma che essa sia unica. L’idea alla base del metodo della bisezione è molto semplice: si dimezza l’intervallo

I0

= a,b[ ] e si cerca il semi-intervallo

I1 ai cui estremi i valori di

f x( ) hanno

6


segni opposti. Una volta trovato l’intervallo

I1, lo si dimezza e si sceglie il semi-intervallo

I2 ai cui estremi i valori di

f x( ) hanno segni opposti, e cosi via. Questo il metodo della bisezione.

Figura 5.4 Descrizione grafica del metodo di bisezione. Questa procedura è descritta graficamente in Figura 5.4, dove con

xk abbiamo

indicato il punto medio dell’intervallo

Ik. La soluzione esatta è rappresentata

dall’intersezione della curva della funzione f con l’asse delle ascisse. Si evince immediatamente che la successione

x0, x1, x

2,..., x

k,... (5.14)

converge ineluttabilmente verso la soluzione

u per

k!" ,

limk!"

xk= u , (5.15)

in quanto l’intervallo

Ik per ogni

k contiene la soluzione

u e la sua lunghezza ad ogni passo viene dimezzata.

Indichiamo con Ik

la lunghezza dell’intervallo Ik(con abuso di notazione) e

introduciamo l’errore ek così definito

e

k! x

k" u . (5.16)

Essendo

7


ek!1

2Ik

, (5.17)

l’iterazione può essere arrestata al primo intero N tale che

IN! " , (5.18)

dove ! è un parametro assegnato. Per l’errore al passo N si ha

eN!"

2. (5.19)

Questo può essere un criterio di arresto per l’iterazione di bisezione.

Nel Programma 5.1 proponiamo una implementazione al calcolatore del metodo della bisezione utilizzando il linguaggio di programmazione MATLAB , [5.2].

Le variabili di ingresso sono: “f” è la stringa che rappresenta, in forma simbolica, la funzione f x( ) , “a” e “b” sono gli estremi dell’intervallo I

0, “toll” è l’errore ! e “nmax”

è il numero massimo di iterazioni disposti ad effettuare per determinare la soluzione. Le variabili di uscita sono: “zero” contiene la soluzione calcolata e “niter” il numero di

iterazioni effettuate per determinarla secondo il criterio di arresto (5.19).

Programma 5.1: Bisezione function [zero,niter]=bisezione(f,a,b,toll,nmax) %Metodo della bisezione x=[a,(a+b)*0.5,b];fx=eval(f);niter=0;I=(b-a)*0.5; while I>=toll&niter<=nmax niter=niter+1; if sign(fx(1))*sign(fx(2))<0 %Il comando “sign(fx)” serve per trovare il segno di fx; x(3)=x(2); x(2)=x(1) +(x(3)-x(1))*0.5; fx=eval(f); I=(x(3)-x(1))*0.5; elseif sign(fx(2))*sign(fx(3))<0 x(1)=x(2); x(2)=x(1) +(x(3)-x(1))*0.5; fx=eval(f); I=(x(3)-x(1))*0.5; else, x(2)=x(find(fx==0));I=0;end %Il comando “find(fx==0)” serve per trovare gli indici del vettore fx corrispondenti a %elementi nulli; end; zero=x(2); return

8


Esercizio 5.1

Risolvere con il metodo della bisezione il circuito di Figura 5.1 assumendo che il resistore non lineare sia un diodo a giunzione pn, R = 1! e E = 1V . Ripetere il calcolo variando i parametri del circuito.

♦

Con il metodo della bisezione è possibile determinare tutte le soluzioni reali di un’equazione algebrica scalare scegliendo in modo opportuno l’intervallo iniziale I

0.

Esercizio 5.2 Risolvere con il metodo della bisezione il circuito di Figura 5.1 assumendo che il

resistore non lineare sia un diodo tunnel e considerando due insiemi di valori per i parametri: E = 12V , R = 6! ; E = 15V , R = 6! . Nel primo caso il circuito ha una sola soluzione, mentre nell’altro ne ha tre. Ciò può essere facilmente verificato attraverso il metodo grafico.

♦ Problema 5.1 Il circuito di Figura 5.1 può essere considerato come un sistema ingresso-uscita,

assumendo, ad esempio, come ingresso la tensione del generatore e come uscita l’intensità di corrente. Determinare la caratteristica ingresso-uscita al variare della resistenza R sia nel caso in cui N è un diodo a giunzione pn, sia nel caso in cui è un diodo tunnel.

♦

Problema 5.2 Studiare con il metodo della bisezione il circuito raddrizzatore rappresentato in Figura

5.5 al variare dei parametri Em,! ,R . In particolare, determinare l’andamento del valore

medio dell’intensità di corrente e delle ampiezze delle prime tre armoniche al variare dei parametri E

m,! ,R .

Figura 5.5 Circuito raddrizzatore

♦

9


5.1.2 Metodo di Newton Raphson

Indichiamo con G la curva che rappresenta funzione f x( ) sul piano x, y( ) , Figura 5.6. La soluzione u dell’equazione (5.4) è rappresentata dal punto P di intersezione della curva G con l’asse delle ascisse. Se si dispone di una “stima” approssimata x

k della

soluzione, un modo semplice ed efficiente di ottenere una stima successiva xk+1

è quello di considerare il punto di intersezione della retta

R

k, tangente alla curva G nel punto

Pk ! xk , f xk( )( ) , con l’asse delle ascisse. L’equazione di R

k è

y x( ) = f xk( ) + !f xk( ) x " xk( ) . (5.20)

(Con il simbolo !f indichiamo la derivata prima di f rispetto ad x .) Si sceglie come stima successiva x

k+1 la soluzione dell’equazione

y xk+1( ) = 0 . (5.21)

La soluzione della (5.21) è

xk+1 = xk !f xk( )"f xk( )

(5.22)

purché !f xk( ) " 0 . E’ come se avessimo “approssimato” la funzione f x( ) con la sua espansione in serie di Taylor nell’intorno di x

k ed arrestato lo sviluppo al secondo

termine. La legge (5.22) consente di calcolare una successione di valori x1, x

2,... a partire

da una stima iniziale x0

. Questo è il metodo di Newton-Raphson. Si osservi che la soluzione u dell’equazione (5.4) è un punto fisso della funzione di iterazione

! xk( ) = xk "f xk( )#f xk( )

, (5.23)

perché ! u( ) = u . La Figura 5.6 mette bene in evidenza, seppure solo qualitativamente, il meccanismo di convergenza alla base del metodo di Newton-Raphson. In seguito mostreremo che l’iterazione generata dalla (5.22) converge sempre alla soluzione u per k!" se x

0 è

sufficientemente vicino alla soluzione. Il metodo di Newton-Raphson, confrontato con gli altri due metodi, converge più rapidamente alla soluzione.

10


Figura 5.6 Interpretazione grafica del metodo di Newton-Raphson.

Osservazione

Il metodo di Newton-Raphson converge in un solo passo alla soluzione esatta se la

funzione f è lineare affine in x , cioè f x( ) = ax + b . (5.24)

Al primo passo si ha x

1= !b / a che è proprio la soluzione dell’equazione. In generale,

invece, è necessario un numero infinito di iterazioni. ♦

Esistono due criteri di arresto semplici ed efficienti per il metodo di Newton-Raphson.

Essi si basano sul controllo della variazione della soluzione approssimata tra due iterate successive

!xk" x

k# x

k#1, (5.25)

e sul controllo del residuo dell’equazione (5.4) così definito

rk ! f xk( ) . (5.26) Per x

k= u si avrebbe sia !x

k+1= 0 che r

k= 0 . Se si usa il controllo della variazione

della soluzione il processo si arresta quando

!xN" # (5.27)

dove ! è un valore assegnato. Invece, se si usa il controllo del residuo il processo si arresta quando

11


rN! " . (5.28)

Figura 5.7

In prossimità della soluzione u dalla (5.22) abbiamo che rk = ! "f u( )#xk+1 . (5.29)

E’ evidente, allora, che i risultati ottenuti con i due criteri di arresto possono essere abbastanza diversi. Nei casi in cui !f u( ) >> 1 (Figura 5.7a) un valore piccolo di !x

N

non esclude che rN

possa essere molto grande, mentre nei casi in cui !f u( ) << 1

(Figura 5.7b) un valore molto basso del residuo non esclude che !xN

possa essere molto elevato. Invece, nei casi in cui !f u( ) " 1 , !x

N e r

N sono dello stesso ordine di

grandezza. A differenza del criterio di arresto che abbiamo introdotto per il metodo della bisezione, i criteri di arresto che abbiamo appena descritto non consentono, almeno in modo diretto, di avere una stima dell’errore e

N. Su questa questione ritorneremo in

seguito (vedi § 5.2) quando studieremo le proprietà di convergenza.

Nel Programma 5.2 è riportata una implementazione in ambiente MATLAB del metodo di Newton-Raphson [5.2], utilizzando come criterio di arresto il controllo della variazione della soluzione (5.27).

Le variabili di ingresso sono: “f” è la stringa che rappresenta, in forma simbolica, la funzione f x( ) , “df“ è la stringa che rappresenta la derivata prima sempre in forma simbolica, “x0” è il dato iniziale, “toll” è l’errore ! e “nmax” è il numero massimo di iterazioni disposti a fare.

Le variabili di uscita sono: “zero” contiene la soluzione calcolata e “niter” il numero di iterazioni effettuate per determinarla secondo il criterio di arresto (5.27).

12


Programma 5.2: Newton-Raphson function [zero,niter]=newton1D(f,df,x0,toll,nmax) % Metodo di Newton-Raphson scalre x=x0; fx=eval(f);dfx=eval(df); niter=0; diff=toll+1; while diff>=toll&niter<=nmax niter=niter+1;diff=-fx/dfx; x=x+diff; diff=abs(diff); fx=eval(f); dfx=eval(df); end; zero=x; return

Esercizio 5.3

Risolvere il problema descritto nell’Esercizio 5.1 con il metodo di Newton-Raphson. Ripetere il calcolo variando x

0.

♦ Esercizio 5.4

Risolvere con il metodo di Newton-Raphson l’equazione

arctan x( ) = 0 . (5.30)

Questa equazione nell’intervallo !" / 2,+" / 2[ ] ha la soluzione x = 0 . Verificare che l’iterazione di Newton-Raphson non sempre converge alla soluzione.

♦

Con il metodo di Newton-Raphson è possibile determinare tutte le soluzioni di un’equazione algebrica non lineare scegliendo opportunamente il valore iniziale x

0.

Esercizio 5.5

Risolvere il problema descritto nell’Esercizio 5.2 con il metodo di Newton-Raphson.

♦

13


Figura 5.8 L’iterazione di Newton-Raphson non converge se x

0= x

b.

Come vedremo nella prossima sezione, il metodo di Newton-Raphson converge alla

soluzione indipendentemente dal valore di x0

se la concavità di f x( ) non cambia. Invece, se cambia la concavità di f x( ) l’iterazione di Newton-Raphson converge solo se x0

è sufficientemente vicino alla soluzione. In Figura 5.8 sono illustrati due casi. In Figura 5.8a , l’iterazione può dare origine ad una sequenza oscillante limitata che non converge o converge lentamente alla soluzione, mentre in Figura 5.8b l’iterazione diverge se x

0 non è sufficientemente vicino alla soluzione. Questi tipi di problemi possono essere

risolti con il metodo della continuazione o con opportune modifiche dello stesso metodo di Newton-Raphson.

5.1.2.1 Metodo di Newton-Raphson modificato con controllo del residuo

L’“impantanamento” dell’iterazione di Newton-Raphson prodotto da oscillazioni,

come, ad esempio, accade nel caso illustrato in Figura 5.8a) è dovuto al fatto che !f cambia segno. In queste situazioni il valore assoluto del residuo oscilla, quindi r

k

durante alcuni passi cresce e in altri decresce. Un modo semplice per ovviare a questo problema consiste nel bloccare il valore di !f nella (5.22) appena che r

k inizia a

crescere al valore che aveva prima dell’inizio della crescita, Figura 5.9. La successione x0, x1, x

2, x

3, x

4,... generata in questo modo converge alla soluzione u . Questa variante

del metodo di Newton-Raphson prende il nome di metodo di Newton-Raphson modificato. Esso converge più lentamente del metodo di Newton-Raphson. Una volta superato l’impantanamento si può riattivare il metodo di Newton-Raphson se si vuole una velocità di convergenza maggiore.

14


Figura 5.9 Metodo di Newton-Raphson modificato: le rette

M

k con k = 3, 4,... hanno la stessa

pendenza della retta R2. Si noti che r !x

3( ) > r x2( ) .

♦ Esercizio 5.6

Risolvere il problema descritto nell’Esercizio 5.1 con il metodo di Newton-Raphson

modificato. Confrontare la velocità di convergenza del Newton-Raphson modificato con la velocità di convergenza del metodo di Newton-Raphson riportando l’andamento del residuo in funzione del numero di iterazioni.

♦

5.1.2.2 Metodo della continuazione Si consideri l’equazione f u;!( ) = 0 (5.31)

e si assuma che la soluzione u vari con continuità al variare del parametro ! . Si assuma, inoltre, che per un dato valore di ! , che indichiamo con !

0, sia nota la soluzione

dell’equazione (5.31) e indichiamola con u0

. Attraverso il metodo di Newton-Raphson possiamo determinare la funzione

u = u !( ) (5.32) a partire dalla soluzione u

0= u !

0( ) .

15


Figura 5.10 Metodo della continuazione

Si consideri la soluzione u1 dell’equazione (5.31) per !

1= !

0+ "!

1. Se !"

1 è

sufficientemente piccolo, la soluzione u1 è sufficientemente vicina alla soluzione u

0 e,

quindi, può essere determinata con il metodo di Newton-Raphson assumendo come dato iniziale proprio u

0, Figura 5.10. In questo caso, indipendentemente dalle proprietà della

funzione f , l’iterazione non solo converge alla soluzione, ma converge anche rapidamente. Una volta determinata u

1 possiamo determinare la soluzione u

2

dell’equazione (5.31) per !2= !

1+ "!

2 applicando di nuovo il metodo di Newton-

Raphson e assumendo, questa volta, come dato iniziale proprio u1. Se !"

2 è

sufficientemente piccolo, la soluzione u2

è sufficientemente vicina a u1 e, quindi,

l’iterazione converge rapidamente alla soluzione. Iterando questa procedura possiamo determinare le soluzioni dell’equazione (5.31) per ! appartenente a un dato intervallo. Questo è il metodo della continuazione. La soluzione u = u !( ) è costruita per continuazione analitica a partire dalla soluzione u

0= u !

0( ) . Si osservi che, nella procedura della continuazione possiamo utilizzare un qualsiasi metodo iterativo per la soluzione dell’equazione (5.31) purché sia garantita la convergenza alla soluzione quando il dato iniziale è sufficientemente vicino ad essa. Problema 5.3

Risolvere i Problemi 5.1 e 5.2 con il metodo della continuazione e il metodo di Newton-Raphson.

♦

16


5.1.4 Metodo di Picard

Il metodo di Picard si basa sulla riformulazione della (5.4) in un’equazione del tipo u = ! u( ) , (5.33)

quindi la soluzione dell’equazione (5.4) è determinata cercando il punto fisso u della funzione di iterazione ! x( ) . Una possibile funzione di iterazione è

! x( ) = x + f x( ) . (5.34) E’ evidente che questa non è l’unica possibile, la scelta è molto più ampia. Questo è un punto molto importante e su di esso torneremo in seguito.

Il punto fisso della funzione di iterazione ! x( ) può essere determinato attraverso lo schema iterativo

xk+1

= ! xk( ) , (5.35)

che graficamente si interpreta come in Figura 5.11. Questo è il metodo di Picard (o metodo del punto fisso).

Figura 5.11 Iterazione di punto fisso: l’iterazione converge al punto per qualsiasi valore di x

0

(a) e per nessun valore di x0

(b).

Non tutte le funzioni di iterazione ! x( ) garantiscono la convergenza verso il punto fisso u . In Figura 5.11 sono riportati due casi: in Figura 5.11a l’iterazione converge al punto fisso u indipendentemente dal valore di x

0, mentre nel caso rappresentato in

17


Figura 5.11b l’iterazione diverge comunque si scelga il valore di x0

, anche quando x0

è arbitrariamente vicino al punto fisso.

Esercizio 5.7

Verificare che l’iterazione x

k= !x

k"1+ b (5.36)

converge al punto fisso

u =b

1! " (5.37)

se ! < 1 . ♦

Rispetto al metodo di Newton-Raphson il metodo di Picard è più semplice da

implementare (in particolare, come poi vedremo, per le equazioni vettoriali), ma come abbiamo appena visto può non convergere alla soluzione anche quando il valore iniziale x0

è arbitrariamente vicino alla soluzione. Inoltre, come in seguito vedremo, quando il metodo di Picard converge esso converge più lentamente del metodo di Newton-Raphson.

Esercizio 5.8

Verificare che con questa scelta della funzione di iterazione

! x( ) = x "f x( )

#f x( ) con !f x( ) " 0 , (5.38)

il metodo di Picard dà il metodo di Newton-Raphson.

♦

Per il metodo di Picard si possono utilizzare gli stessi criteri di arresto introdotti per il metodo di Newton-Raphson. Il residuo dell’iterazione è dato da

r

k= ! x

k( ) " xk . (5.39) Essendo ! x

k( ) = xk+1 , nel metodo di Picard il residuo definito in base alla (5.39) è uguale all’incremento

rk= !x

k. (5.40)

18


Nel Programma 5.3 è riportata una implementazione in ambiente MATLAB del metodo di Picard utilizzando come criterio di arresto il controllo della variazione della soluzione.

Le variabili di ingresso sono: “f” è la stringa che rappresenta, in forma simbolica, la funzione di iterazione, “x0” è il dato iniziale, “toll” è l’errore ! e “nmax” è il numero massimo di iterazioni accettate.

Le variabili di uscita sono: “zero” contiene il punto fisso e “niter” il numero di iterazioni effettuate per determinarla secondo il criterio di arresto (5.27). Programma 5.3: Picard function [zero,niter]=picard1D(f,x0,toll,nmax) % Metodo di Picard x=x0; fx=eval(f); niter=0; diff=toll+1; while diff>=toll&niter<=nmax niter=niter+1;diff=fx; x=x+fx; diff=abs(diff); fx=eval(f); end; zero=x; return Esercizio 5.9

L’equazione cos x ! x = 0 (5.41)

ha una sola soluzione. Determinarla con il metodo di Picard. Risolvere la stessa equazione con il metodo della bisezione e il metodo di Newton-Raphson, e studiare la velocità di convergenza dei tre metodi rappresentando l’andamento dei residui al variare del numero di iterazioni. Quale è il metodo più veloce?

♦

Esercizio 5.10

Verificare che non è possibile risolvere il problema dell’Esercizio 5.1 con il metodo di Picard perché l’iterazione non converge per nessun valore del dato iniziale.

♦ Nella prossima sezione ci occuperemo del problema della convergenza e, in

particolare, del problema della convergenza del metodo di Picard. Il metodo di Picard può essere riformulato in modo tale da evitare il problema incontrato nell’Esercizio 5.10.

19


5.2 Proprietà di convergenza dei metodi iterativi scalari Si assuma che l’equazione (5.4) abbia una sola soluzione u in a,b[ ] . Sia x

0, x1, x

2,..., x

k,... (5.42)

la successione generata risolvendo l’equazione (5.4) attraverso uno dei tre metodi iterativi (bisezione, Picard, Newton-Raphson) che abbiamo descritto nel precedente paragrafo. Ora affronteremo le seguenti questioni (fondamentali per ogni metodo iterativo):

i) la successione converge alla soluzione u ? ii) con quale velocità l’errore decresce? iii) che errore si commette arrestando l’iterazione dopo un numero assegnato di

passi? L’errore assoluto e

k al passo k è definito nella (5.16), l’incremento !x

k è definito nella

(5.25) e il residuo rk è definito nella (5.26). Se

lim

k!"

ek= 0 (5.43)

l’iterazione converge alla soluzione u . In queste condizioni si ha immediatamente che

limk!"

rk= 0 (5.44)

e lim

k!"

#xk= 0 . (5.45)

Se l’iterazione converge alla soluzione al crescere di k l’errore diminuisce. Si

consideri a tal proposito il rapporto tra l’errore in valore assoluto al passo k e l’errore, sempre in valore assoluto, al passo k !1 ,

!k"

ek

ek#1

. (5.46)

Se esiste un intero p e una costante positiva C tali che

limk!"

#k = C ek$1p$1( ) (5.47)

p prende il nome di ordine di convergenza del metodo iterativo e C prende il nome di

costante asintotica di errore. Dalla (5.47) si ha che per N sufficientemente grande (errore asintotico)

20


eN ! C eN "1

p . (5.48) La (5.48) è la legge che governa l’andamento asintotico dell’errore. L’errore di un metodo iterativo con ordine di convergenza p più elevato si riduce più velocemente dell’errore di un metodo con p più basso.

I criteri di arresto che abbiamo impiegato nel metodo di Picard e nel metodo di Newton-Raphson si basano sul controllo del valore assoluto dell’incremento !x

k o sul

controllo del residuo rk

. Pertanto, è importante cercare una relazione tra l’errore ek

e questi due parametri.

Analizzeremo insieme le iterazioni generate con il metodo di Picard e il metodo di Newton-Raphson perché, come abbiamo già visto, l’iterazione di Newton-Raphson può essere interpretata come un particolare metodo di Picard. L’iterazione generata attraverso la bisezione sarà trattata separatamente perché, a differenza degli altri due metodi, con questo metodo è possibile avere una stima diretta dell’errore.

5.2.1 Metodo della bisezione

Nel metodo della bisezione attraverso la relazione (5.17) è possibile stimare il valore

assoluto dell’errore ek. Essendo, per costruzione,

Ik=1

2Ik!1

(5.49)

si ha

Ik=1

2

!"#

$%&1

Ik'1 =

1

2

!"#

$%&2

Ik'2 = ... =

1

2

!"#

$%&k

I0

(5.50)

dove I

0= b ! a( ) . Combinando le (5.17) e (5.50) otteniamo

ek!1

2

"#$

%&'k+1

b ( a( ) . (5.51)

Essendo lim

k!"

0.5( )k

= 0 , dalla (5.51) si ha che limk!"

ek= 0 e, quindi, l’iterazione generata

dalla bisezione converge alla soluzione u . Esercizio 5.11

Verificare che il numero di iterazioni N che garantisce eN! " , dove ! è un numero

positivo assegnato, è dato da

21


N > log2b ! a"

#$%

&'(!1

)

*+

,

-. . (5.52)

♦

Se considerassimo la (5.17) con il segno di eguaglianza avremmo una sovrastima

dell’errore. Indichiamo con !ek questa sovrastima. Per essa abbiamo

!ek=1

2Ik

. (5.53)

Utilizzando la (5.49) dalla (5.53) otteniamo

!ek

!ek!1

=1

2. (5.54)

Di conseguenza l’ordine di convergenza del metodo della bisezione non è inferiore all’unità e la corrispondente costante asintotica di errore è uguale a 0.5 . In conclusione, il metodo di bisezione è un metodo del primo ordine. 5.2.2 Metodo di Picard

Analizziamo ora le proprietà della successione (5.42) generate dall’iterazione di punto

fisso x

k= ! x

k"1( ) . (5.55) Per ipotesi esiste un solo valore u ! a,b[ ] tale che

u = ! u( ) (5.56) Prima di iniziare a studiare le proprietà di convergenza, enunceremo un teorema che utilizzeremo spesso, il teorema della media per le funzioni scalari. Se la funzione ! = ! x( ) è differenziabile in a,b[ ] , allora per ogni x, y ! a,b[ ] esiste un 0 ! t ! 1 tale che

! x( ) "! y( ) = #! $( ) x " y( ) . (5.57) dove ! = x + t y " x( ) 7. Posto

7 Questo teorema è un caso particolare dell’espansione in serie di Taylor: data una funzione f = f x( )

continua e n volte derivabile nell’intervallo a,b[ ] esiste un numero reale 0 ! t ! 1 tale che

22


L = sup

x! a,b[ ]"# x( ) (5.58)

dalla (5.57) si ha che per ogni x, y ! a,b[ ]

! x( ) "! y( ) # L x " y . (5.59) Essendo, per ipotesi, ! differenziabile in a,b[ ] , la costante L è limitata se l’intervallo a,b[ ] è limitato. In questo caso si dice che ! è lipschitziana8. Se L < 1 la funzione ! si

dice che è L-contrattiva in a,b[ ] . Esercizio 5.12

Verificare che la funzione di iterazione nell’Esercizio 5.9 è L-contrattiva nell’intervallo !" ,+"( ) (aperto sia a destra che a sinistra).

♦

Esercizio 5.13

Verificare che la funzione di iterazione nel problema affrontato nell’Esercizio 5.10 non è L-contrattiva.

♦

5.2.2.1 Convergenza

Sottraendo membro a membro le (5.55) e (5.56) otteniamo e

k= ! x

k"1( ) "! u( ) . (5.60) Utilizzando la (5.58) dalla (5.60) abbiamo

ek! L e

k"1. (5.61)

Allora è immediato verificare che f x( ) = f a( ) + !f a( ) x " a( ) + 1

2!!f a( ) x " a( )

2

+ ...+ 1

n"1( )! fn"1( )

a( ) x " a( )2

+ Rn x( ) dove ! = a + t x " a( ) e

Rn x( ) = 1

n! fn( ) !( ) x " a( )

n . L’espressione di Rnx( ) è la formula del resto di Lagrange.

8 Si dice che la funzione f = f x( ) è lipschitziana nel dominio a,b[ ] se esiste una costante K limitata tale

che f x( ) ! f y( ) " K x ! y per ogni x, y ! a,b[ ] .

23


e

k! L

ke0

(5.62) per k = 1,2,... . Se la funzione di iterazione è L-contrattiva in a,b[ ] si ha L < 1 e dalla (5.62) segue che

limk!"

ek= 0 . (5.63)

In conclusione, se la funzione di iterazione è L-contrattiva nell’intervallo a,b[ ] che contiene il punto fisso u e x

0! a,b[ ] , l’iterazione di Picard converge al punto fisso u .

La funzione di iterazione è L-contrattiva in a,b[ ] se il modulo della sua derivata prima è inferiore a 1 in a,b[ ] . 5.2.2.2 Velocità di convergenza

Dal teorema del valore medio sappiamo che esiste un tk!1

" 0,1( ) tale che ! x

k"1( ) "! u( ) = #! $k"1( ) xk"1 " u( ) , (5.64)

dove !

k"1 = xk"1 + tk"1 u " xk"1( ) . Allora, combinando le (5.60) e (5.64) otteniamo

ek

ek!1

= "# $k!1( ) . (5.65)

Siccome lim

k!"

xk= u si ha che

lim

k!"xk#1 + tk#1 u # xk#1( )$% &' = u (5.66)

e, quindi,

limk!"

#$ %k&1( ) = #$ u( ) . (5.67)

Combinando le (5.65) e (5.67) otteniamo

limk!"

ek

ek#1

= $% u( ) . (5.68)

24


Dunque, se !" u( ) # 0 l’iterazione di Picard ha un ordine di convergenza unitario, quindi è un metodo del primo ordine come quello della bisezione. La costante asintotica di errore è uguale a !" u( ) . 5.2.2.3 Stima dell’errore

Ora faremo vedere che ek può essere espresso in funzione di !x

k+1. Potendo scrivere

x

k! u = x

k+1! u( ) ! x

k+1! x

k( ) (5.69) abbiamo che

ek= ! x

k( ) "! u( )#$ %& " 'xk+1 . (5.70)

Dal teorema del valore medio abbiamo che esiste un tk! 0,1( ) tale che

! xk( ) "! u( ) = #! $

k( ) xk " u( ) (5.71) dove !

k= x

k+ t

ku " x

k( ) . Sostituendo la (5.71) nella (5.70) abbiamo

ek=

!xk+1

"# $k( ) %1

. (5.72)

E’ evidente, allora, che se !" # +1 in a,b[ ] e limk!"

#xk= 0 , allora lim

k!"

ek= 0 e l’iterazione

converge al punto fisso u . Esercizio 5.14

Dimostrare che limk!"

#xk= 0 se la funzione di iterazione è L-contrattiva nel dominio

a,b[ ] e x0! a,b[ ] .

♦ Siccome lim

k!"#$ %

k( ) = #$ u( ) , esiste un N tale che per k ! N si ha !" #k( ) $ !" u( ) ,

quindi

eN!

"xN +1

1# $% u( ). (5.73)

Questa relazione ci consente di stimare l’errore e

N per un dato !x

N +1.

25


Osservazione

Abbiamo mostrato che il metodo di Picard può essere impiegato solo se la funzione di iterazione! è L-contrattiva nel dominio che contiene il punto fisso. Affinché la funzione ! data dalla (5.34) sia L-contrattiva è necessario che !f x( ) sia minore di zero. Questo vincolo rende praticamente inapplicabile il metodo di Picard per risolvere un’equazione del tipo (5.9). Questa limitazione può essere superata osservando che esiste una riformulazione più generale dell’equazione (5.4) come problema di punto fisso. L’equazione (5.4) può essere riscritta come

x = x +! f x( ) , (5.74)

dove ! è un parametro, per il momento arbitrario, diverso da zero. In questo caso la funzione di iterazione è data da

! x( ) = x +" f x( ) . (5.75) Possiamo pensare di scegliere ! in modo tale che questa funzione di iterazione sia L-contrattiva. Imponendo che la (5.75) verifichi la condizione !" x( ) < 1 si ottiene

1+! "f x( ) < 1 . (5.76) La (5.76) è verificata se

!2 < "f x( )# < 0 . (5.77) Esercizio 5.15

Risolvere il problema dell’Esercizio 5.10 con la funzione di iterazione (5.75), scegliendo il valore di ! in modo tale che ! sia L-contrattiva almeno in un intorno della soluzione.

♦ 5.2.3 Metodo di Newton-Raphson

La scelta più generale di funzione di iterazione per risolvere l’equazione (5.4) è data da

! x( ) = x +" x( ) f x( ) (5.78)

dove, ora, ! = ! x( ) è una funzione diversa da zero e limitata nell’intervallo a,b[ ] , da scegliere in modo tale che !" x( ) < 1 almeno in un intorno della soluzione u . Essendo

26


!" = 1+ !# f +# !f (5.79)

è conveniente scegliere ! in modo tale che

! = "1

#f (5.80)

purché !f " 0 e limitato. Con questa scelta l’iterazione di punto fisso (5.55) si riduce proprio all’iterazione di Newton-Raphson e

!" = #f x( )

!f x( ). (5.81)

In questo caso si ha !" u( ) = 0 , quindi esiste un intorno J del punto fisso u in cui !" x( ) < 1 . Scegliendo

x0!J l’iterazione converge alla soluzione u .

Riportiamo un teorema attraverso cui è possibile stimare la regione a cui deve appartenere il dato iniziale per avere un’iterazione convergente.

Teorema di convergenza globale [5.3] Sia f x( ) una funzione continua con la sua derivata prima e seconda definita nell’intervallo chiuso e limitato a,b[ ] . Se:

1. f a( ) f b( ) < 0 2. !f x( ) " 0 #x $ a,b[ ] 3. !!f x( ) " 0, oppure !!f x( ) # 0 $x % a,b[ ]

4. f a( )

!f a( )< b " a e

f b( )

!f b( )< b " a

allora il metodo di Newton-Raphson converge all’unica soluzione u in a,b[ ] , per ogni scelta di x

0! a,b[ ] .

La condizione 4. significa che la tangente negli estremi dell’intervallo a,b[ ] interseca l’asse x all’interno dell’intervallo. Problema 5.4

Dimostrare il Teorema di convergenza globale per il metodo di Newton-Raphson. ♦

27


Ora mostreremo che l’ordine di convergenza per il metodo di Newton-Raphson è

uguale a 2. Il metodo di Newton-Raphson è un metodo del secondo ordine. Sviluppiamo la funzione ! nell’intorno di x = u utilizzando l’espansione di Taylor al secondo ordine. Esiste un 0 ! t ! 1 tale che

! x( ) = ! u( ) + "! u( ) x # u( ) +1

2""! $( ) x # u( )

2 (5.82)

dove ! = x + t u " x( ) . Essendo ! u( ) = u e !" u( ) = 0 dalla (5.82) otteniamo

! xk( ) " u =

1

2##! $

k( ) xk " u( )2 , (5.83)

ovvero

ek+1

=1

2!!" #

k( )ek2 . (5.84)

Essendo lim

k!"#k= u , dalla (5.84) otteniamo

limk!"

ek+1

ek

=1

2##$ u( ) e

k. (5.85)

Pertanto, l’ordine di convergenza del metodo di Newton-Raphson è uguale a 2. La costante asintotica di errore è uguale a 0.5 !!" u( ) . E’ immediato verificare che

!!" u( ) =!!f u( )

!f u( ). (5.86)

Dalla (5.85) segue che

ek=

1

KKe

M

2k per k ! M (5.87)

dove K = 0.5 !!" u( ) e M è un numero intero sufficientemente elevato. E’ evidente, allora, che l’errore tende a zero solo se e

M< 1 / K . Ciò è in accordo con il fatto che l’iterazione

di Newton-Raphson converge alla soluzione solo se il punto di partenza è ad essa sufficientemente vicino.

Infine, resta da determinare la relazione tra l’errore ek

e l’incremento !xk+1

nell’intorno della soluzione. Dalla (5.73) si ha

28


eN! "x

N +1 (5.88)

perché in questo caso !" u( ) = 0 . 5.3 Sistemi di equazioni algebriche non lineari

Il metodo di Picard e il metodo di Newton-Raphson si estendono senza alcuna difficoltà ai sistemi di equazioni algebriche non lineari. L’estensione del metodo della bisezione, invece, presenta difficoltà insormontabili. In questo paragrafo estenderemo il metodo di Picard e il metodo di Newton-Raphson alla soluzione di un sistema del tipo (5.3).

Esempio 5.2

Si consideri il circuito resistivo non lineare illustrato in Figura 5.12. Risolviamolo applicando il metodo dei potenziali ai nodi.

Il funzionamento del transistore è descritto attraverso il modello di Ebers-Moll. Si ha ib = gb wb ,wc( ) = 1!"F( ) IES e

wb /VT !1( ) + 1!"R( ) ICS ewb !wc( )/VT

!1( ) (5.89)

ic = gc wb ,wc( ) = !F IES ewb /VT "1( ) " ICS e

wb "wc( )/VT"1( ) (5.90)

dove V

T= 26 mV , !

R= 0.5 ÷ 0.8 , !

F= 0.99 , I

ES! I

CS! 10

"12÷10

"10A a 25°C . Per

ispezione diretta del circuito otteniamo le equazioni

ib=E1! w

b

R1

, (5.91)

ic=E2! w

c

R2

. (5.92)

Combinando le equazioni (5.89)-(5.92) otteniamo il sistema di due equazioni algebriche non lineari

wb

R1

+ gb wb ,wc( ) !E1

R1

= 0

gc wb ,wc( ) +wc

R2

!E2

R2

= 0.

"

#

$$

%

$$

(5.93)

29


Figura 5.12 Circuito di polarizzazione per la configurazione di emettitore comune.

♦

Prima di estendere il metodo di Picard e il metodo di Newton-Raphson ai sistemi di equazioni, vogliamo fare qualche ulteriore considerazione sull’esistenza e unicità della soluzione dell’equazione (5.3). Una funzione f :D ! R

n" C ! R

n è strettamente monotona in D se

f x( ) ! f y( )"# $%T

x ! y( ) > 0 , !x " y #D . (5.94) Se f è strettamente monotona in D , allora essa definisce una corrispondenza biunivoca tra il dominio di definizione D e il codominio C . In conclusione, se C contiene lo zero e f è strettamente monotona, allora l’equazione (5.3) ha una ed una sola soluzione in D . La funzione f è strettamente monotona in D se la matrice jacobiana

J f x( )

J f x( ) !

"f1/ "x

1"f1/ "x

2... "f

1/ "xn

"f2/ "x

1"f2/ "x

2... "f

2/ "xn

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

"fn / "x1 "fn / "x2 ... "fn / "xn

(5.95)

è definita positiva per ogni x !D , [5.1]. Lasciamo al Lettore la dimostrazione di questa proprietà. Pertanto, se il codominio di f contiene lo zero e

J f x( ) è definita positiva in

D , la soluzione dell’equazione (5.3) esiste ed è unica. Esiste un altro importante risultato sull’esistenza e unicità della soluzione dell’equazione (5.3). Esso deriva dal teorema di Hadamard, [5.1]. Se f è continua e differenziabile in !n e J f

!1x( ) " # < $ per ogni x !"

n , allora f è un omeomorfismo

di !n in !n . Di conseguenza esiste ed è unica la soluzione dell’equazione (5.3). Se

J f x( ) diventa singolare in qualche punto di D , allora la soluzione può non essere unica.

30


5.3.1 Metodo di Picard

Come nel caso scalare, è sempre possibile riscrivere il sistema di equazioni algebriche non lineari (5.3) nella forma

u = ! u( )"

u1= !

1u

1,u

2,...,u

n( )u

2= !

2u

1,u

2,...,u

n( ) ...

un= !

nu

1,u

2,...,u

n( )

#

$

%%

&

%%

(5.96)

dove ! :Rn " R

n è una funzione di iterazione vettoriale opportuna. Una scelta possibile per la funzione di iterazione è

! x( ) = x + f x( ) . (5.97) Ovviamente, come nel caso scalare, altre scelte sono possibili, in seguito le considereremo. La soluzione u del sistema (5.3) è il punto fisso della funzione di iterazione ! . Come nel caso scalare, consideriamo la successione

x1( ),x

2( ),...,x

k( ),... (5.98)

generata da

x

k( )= ! x

k"1( )( ) (5.99) con x 0( ) assegnato. Sotto opportune condizioni sulla funzione di iterazione ! la successione (5.98) converge al punto fisso u . Assumiamo che il metodo di Picard converga. Posto

!x = xk+1( )

" xk( ) (5.100)

e

rk( )= ! x

k( )( ) " x k( ) (5.101) per arrestare l’iterazione possiamo utilizzare uno dei due criteri

!xN( )

" # , (5.102)

rN( )

! " , (5.103)

31


dove con x indichiamo una norma del vettore x (vedi Appendice A2), ad esempio, la norma euclidea

x2= x

i

2

i=1

n

! (5.104)

ed ! è un parametro dato. Se si sceglie la (5.97) come funzione di iterazione si ha

!x

N +1( )= r

N( ) . (5.105) Nel Programma 5.4 è riportata una implementazione in ambiente MATLAB del

metodo di Picard per un sistema di equazioni utilizzando come funzione di iterazione la (5.97) e come criterio di arresto il controllo della variazione della soluzione (5.102).

Le variabili di ingresso sono: “x0” è il dato iniziale , “toll” è l’errore ! e “nmax” è il numero massimo di iterazioni accettate (è impiegata la norma euclidea per valutare l’errore).

Le variabili di uscita sono: “zero” contiene la soluzione calcolata, “niter” il numero di iterazioni effettuate per determinarla secondo il criterio di arresto (5.102). La funzione di iterazione ! è costruita attraverso la function “fun(y)”.

Programma 5.4: Picard vettoriale function [zero,niter,resrel]=picard(x0,toll,nmax) % Metodo di Picard vettoriale. % Attenzione: y è un vettore colonna. y=x0; [fx]=fun(y); niter=0; diff=toll+1; while diff>=toll&niter<=nmax niter=niter+1; diff=fx; y=y+diff; diff=norm(diff); [fx]=fun(y); end; zero=y; return 5.3.2 Metodo di Newton-Raphson

In analogia con il caso scalare, sostituiamo la funzione f x( ) con il suo sviluppo in

serie di Taylor nell’intorno di x k( ) arrestato al termine lineare,

32


fi x( ) ! fi xk( )( ) +

"fi

"x jj=1

n

#x=x

k( ) x j $ x jk( )( ) i = 1,2,...,n . (5.106)

Le n relazioni (5.106) possono essere rappresentate sinteticamente in forma vettoriale in questo modo

f x( ) ! f x k( )( ) + J f x

k( )( ) x " x k( )( ) (5.107) dove

J f x( ) è la matrice jacobiana di f x( ) la cui espressione è data dalla (5.95). La legge

che governa l’iterazione di Newton-Raphson è determinata sostituendo il segno “! ” con il segno “= ” nella (5.106) e sostituendo l’espressione di f così ottenuta nell’equazione (5.3),

J f x

k( )( ) x k+1( )! x

k( )( ) = !f xk( )( ) . (5.108)

Ad ogni passo dell’iterazione bisogna risolvere un sistema di equazioni algebriche lineari nell’incognita x k+1( )

! xk( )( ) . Per questa ragione il costo computazionale per del metodo

di Newton-Raphson è più alto di quello per il metodo di Picard. Risolvendo simbolicamente il sistema (5.108) si ha

x

k+1( )= x

k( )! J f

!1x

k( )( ) f x k( )( ) . (5.109) Questa è la legge che governa l’iterazione di Newton-Raphson per un sistema di equazioni algebriche non lineari. Se si blocca il valore della matrice jacobiana al valore che essa ha assunto a una precedente iterazione si ha il metodo di Newton-Raphson modificato. In questo caso bisogna invertire una sola volta la matrice jacobiana. Osservazione

Se si sceglie come funzione di iterazione nel metodo di Picard la funzione

! x( ) = x " J f

"1x( ) f x( ) (5.110)

si ottiene il metodo di Newton-Raphson.

♦ L’iterazione di Newton-Raphson può essere arrestata utilizzando i criteri (5.102) e (5.103) che abbiamo introdotto per i metodo di Picard. Essendo

33


f x

k( )( ) = J f xk( )( )!x k+1( ) (5.111)

e

Ay ! A y (vedi Appendice 2), si ha

f x

k( )( ) ! J f xk( )( ) "x

k+1( ) . (5.112)

La (5.109) dà un maggiorante della norma del residuo in termini della norma dell’incremento tra due iterate consecutive. Programma 5.5: Newton-Raphson vettoriale function [zero,niter]=newton(x0,toll,nmax) % Metodo di Newton-Raphson vettoriale % Attenzione: y è un vettore colonna y=x0; [fx,dfx]=fun(y); niter=0; diff=toll+1; while diff>=toll&niter<=nmax niter=niter+1; diff=-dfx\fx'; y=y+diff; diff=norm(diff); [fx,dfx,fL]=fun(y); end; zero=y; return

Nel Programma 5.5 è riportata una implementazione in ambiente MATLAB del metodo di Newton-Raphson per un sistema di equazioni utilizzando come criterio di arresto il controllo della variazione della soluzione (5.102). La funzione f e la corrispondente matrice jacobiana sono costruite attraverso la chiamata alla function “fun(y)”. Per la soluzione del sistema (5.108) si utilizza la fattorizzazione LU.

Le variabili di ingresso sono: “x0” è il dato iniziale, “toll” è l’errore ! e “nmax” è il numero massimo di iterazioni accettate.

Le variabili di uscita sono: “zero” contiene la soluzione calcolata e “niter” il numero di iterazioni effettuate per determinarla secondo il criterio di arresto (5.102). Esempio 5.3

Risolviamo il sistema (5.93) sia con il metodo di Picard (Programma 5.4), sia con il metodo di Newton-Raphson (Programma 5.5). Si assumano i seguenti valori per i parametri: R

1= 2! , R

2= 1! , E

1= 1V , E

2= 0.5V , !

R= 0.8 , !

F= 0.9 ,

IES = ICS = 1pA , VT= 26mV . L’iterazione di Picard non converge mai alla soluzione,

comunque si scelga x 0( ) . L’iterazione di Newton-Raphson con x 0( )= 0 , toll = 10!3

34


converge dopo 38 iterazioni alla soluzione wb= 0.722V , w

c= 0.015V . In

corrispondenza di questi valori si ha f = 1.94 !9,7 µA .

Con altri valori dei parametri l’iterazione di Newton-Raphson con x 0( )= 0 può non

convergere. Ad esempio, con R1= 500! , R

2= 200k! , E

1= 0.5V , E

2= 10V dopo tre

iterazioni si ha NaN (in MATLAB “NaN” è un non numero, cioè una forma indeterminata del tipo 0 / 0 o ! /! , [5.2]; nel caso in esame è del tipo ! /! ). Questo problema di convergenza può essere risolto con il metodo della continuazione o con il metodo di Newton-Raphson modificato.

♦

5.4 Il problema della convergenza

Si assuma che l’equazione (5.3) abbia una ed una sola soluzione u in D !Rn . Sia

x

0( ),x

1( ),x

2( ),...,x

k( ),... (5.113)

la successione generata risolvendo l’equazione (5.3) attraverso uno dei due metodi iterativi descritti nel precedente paragrafo (Picard, Newton-Raphson). Come nel caso scalare affronteremo le seguenti questioni:

i) la successione converge alla soluzione u? ii) con quale velocità l’errore decresce? iii) che errore si commette arrestando l’iterazione dopo un numero fissato di passi?

Analogamente al caso scalare, introduciamo l’errore assoluto e k( ) al passo k

ek( )! x

k( )" u . (5.114)

L’incremento !x k( ) è definito nella (5.100) e il residuo r k( ) è definito nella (5.101). Come nel caso scalare, se

limk!"

ek( )= 0 (5.115)

l’iterazione converge alla soluzione u . In queste condizioni si ha immediatamente che

limk!"

rk( )= 0 (5.116)

e lim

k!"

#xk( )= 0 . (5.117)

35


Si consideri il rapporto tra la norma dell’errore al passo k e la norma dell’errore al passo k !1 ,

!k"

ek( )

ek#1( )

. (5.118)

Si assuma che l’iterazione converga alla soluzione u e che esista un intero p e una costante positiva C tali che

limk!"

#k = C ek( ) p$1( )

. (5.119)

Il parametro p non è altro che l’ordine di convergenza del metodo iterativo e C la costante asintotica di errore. La legge che governa il comportamento asintotico della norma dell’errore è, allora,

eN +1( )

! C eN( ) p

. (5.120)

Come nel caso scalare, è importante cercare le relazioni che esistono tra la norma dell’errore a un certo passo e la norma dell’incremento !x o del residuo r . 5.4.1 Iterazione di Picard

Analizziamo ora le proprietà della successione (5.113) generate dall’iterazione di punto fisso

xk= ! x

k"1( ) (5.121) nell’ipotesi che esiste uno ed un solo valore u nel dominio limitato e convesso 9 D ! R

n tale che

u = ! u( ) (5.122) e che ! sia differenziabile.

Il teorema della media nella forma illustrata nel § 5.2 non vale per le funzioni vettoriali di variabili vettoriali10, [5.1]. Prima di iniziare lo studio delle proprietà di

9 Un dominio D si dice convesso se per ogni x, y !D , tutti i punti del “segmento “ definito da

s = x + t y ! x( ) per 0 ! t ! 1 appartengono a D .

36


convergenza del metodo di Picard, enunceremo un’estensione. Per ogni x,y !D esiste un 0 ! t ! 1 tale che [5.1]

! x( ) " ! y( ) # J! $( ) x " y , (5.123)

dove ! = x + t y " x( ) ;

J! indica la norma della matrice jacobiana

J! (vedi Appendice

2). Posto

L = sup

x!D

J" x( ) (5.124)

dalla (5.123) si ha

! x( ) " ! y( ) # L x " y (5.125) per ogni x,y !D . Essendo ! differenziabile in D si ha che la costante L è limitata, quindi, come nel caso scalare, ! è lipschitziana 11. Come nel caso scalare, diciamo che la funzione di iterazione ! è L -contrattiva se L < 1 . 5.4.1.1 Convergenza

Sottraendo membro a membro le equazioni (5.121) e (5.122) otteniamo

ek( )= ! x

k"1( )( ) " ! u( ) . (5.126) Utilizzando, ora, la (5.125) dalla (5.126) otteniamo

ek( )

! L ek"1( ) . (5.127)

Iterando k volte la (5.127) si ha

ek( )

! Lke0( ) . (5.128)

Se L < 1abbiamo che lim

k!"

Lk= 0 e

10 Esso continua a valere per le singole componenti della funzione ! : per ogni x, y !D esiste un

0 ! ti! 1 tale che !i x( ) "!i y( ) =

#!i#x j

x=sixi " yi( )

j=1

n

$ dove si= x + t

iy ! x( ) , però t

i varia al variare di i

i = 1, 2, ...,n( ) . 11 Si dice che la funzione f : Rn

! Rn è lipschitziana in D se esiste una costante finita L tale che per

ogni x, y !D si ha che f x( ) ! f y( ) " L x ! y .

37


limk!"

ek( )

= 0 . (5.129)

Pertanto, come nel caso scalare, se la funzione di iterazione ! è L -contrattiva nel dominio D allora l’iterazione di Picard (5.121) converge al punto fisso u .

La funzione ! è L -contrattiva se

J! x( ) < 1 "x #D . (5.130)

Per ogni ! positivo piccolo a piacere esiste una norma

J! x( )

" tale che [2]

J! "

# $ J!( ) + % , (5.131)

dove

! J"( ) è il raggio spettrale della matrice

J! , che per definizione è

! J"( ) = max

i

#i

J"( ) (5.132)

e

!i

J"( ) è il generico autovalore di J! . Allora se

! J" x( )( ) < 1 #x $D (5.133)

certamente esiste una norma per la quale la (5.130) è verificata e, quindi l’iterazione di punto fisso (5.121) converge.

Con la funzione di iterazione (5.97) si ha

J! = I + J f . (5.134)

Allora, gli autovettori di

J f sono anche autovettori di

J! e per gli autovalori si ha

!i J"( ) = 1+ !i J f( ) . (5.135)

In termini dello jacobiano della funzione f la condizione (5.133) diventa

max

i1+ !i J f x( )( ) < 1 "x #D . (5.136)

Questa condizione può essere verificata solo se gli autovalori di

J f sono tutti minori di

zero e in modulo minori di uno,

38


!1 < "i J f x( )( ) < 0 #x $D e i = 1,2,...,n . (5.137)

Con la funzione di iterazione

! x( ) = x " f x( ) (5.138) la condizione (5.133) diventa

0 < !i J f x( )( ) < 1 "x #D e i = 1,2,...,n . (5.139)

Osservazione

Come nel caso scalare, il metodo di Picard può essere impiegato solo se f verifica una delle due condizioni (5.137), (5.139). Questo vincolo rende praticamente inapplicabile il metodo di Picard per risolvere, ad esempio, un sistema del tipo (5.93).

Questo problema può essere risolto come nel caso scalare, osservando che esiste una riformulazione più generale dell’equazione (5.3) come problema di punto fisso. L’equazione (5.3) può essere riscritta come

x = x + Af x( ) , (5.140)

dove A è matrice invertibile, per il momento arbitrario. In questo caso la funzione di iterazione è data da

! x( ) = x + Af x( ) (5.141)

e

J! = I + AJ f . (5.142)

Si può pensare di scegliere A in modo tale che

max

i!i I + AJ f( )( ) < 1 . (5.143)

Esercizio 5.16

Verificare che la funzione f del problema (5.93) non soddisfa né la condizione (5.137), né la condizione (5.139); si assumano per i parametri del circuito gli stessi valori utilizzati nell’Esempio 5.3. Determinare una matrice A per cui sia verificata la condizione (5.143) e risolvere il sistema il (5.93) con il metodo di Picard.

♦

39


5.4.1.2 Velocità di Convergenza

Assumiamo che l’iterazione di Picard (5.121) converga al punto fisso u . Vogliamo studiarne la velocità di convergenza quando x

k si trova in un intorno di u . Sviluppiamo

in serie di Taylor la funzione !ix1, x

2,..., x

n( ) nell’intorno di u1,u

2,...,u

n( ) e arrestiamo lo sviluppo al termine quadratico. Abbiamo per l = 1,2,...,n

!l x1, x

2,..., xn( ) " !l u1

,u2,...,un( ) +

#!l#x j

x=u

j=1

n

$ x j % uj( ) +

+1

2

#2!l#xi#x jj=1

n

$ x=u

i=1

n

$ xi % ui( ) x j % uj( ). (5.144)

Le n relazioni (5.144) possono essere riscritte in forma vettoriale nel seguente modo

! x( ) " ! u( ) + J! u( ) x # u( ) +

1

2x # u( )

T

H! x # u( ) (5.145)

dove

H! x( ) è la matrice hessiana di ! x( )

H! =

"2!

"x1"x

1

"2!

"x1"x

2

..."2!

"xn"x

1

"2!

"x1"x

2

"2!

"x2"x

2

..."2!

"xn"x

2

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

"2!

"x1"x

n

"2!

"x2"x

n

..."2!

"xn"x

n

. (5.146)

Dalla relazione (5.126) e dall’espansione (5.145) segue che

ek( ) ! J" u( )e k#1( )

+1

2ek#1( )T

H"ek#1( ) . (5.147)

Siccome e k!1( )

" 0 per k!" , per k sufficientemente grande possiamo trascurare nella (5.147) i termini quadratici, quindi

ek( ) ! J" u( )e k#1( ) . (5.148)

Dalla (5.148) si ottiene

40


!k"ek( )

ek#1( )

$ J% u( ) . (5.149)

Dunque, per k sufficientemente grande la norma dell’errore al passo k è inferiore alla norma dell’errore al passo precedente moltiplicata per la norma della matrice

J! u( ) . Di

conseguenza, l’errore decresce almeno linearmente in un intorno sufficientemente piccolo del punto fisso u . In realtà, come ora mostreremo, l’errore decresce linearmente e non più velocemente.

Indichiamo con !1

J"( ),!2 J"( ),...,!n J"( ) gli autovalori della matrice J! u( ) e con

a1,a

2,...,a

n i corrispondenti autovettori. Se gli autovettori di

J! u( ) sono linearmente

indipendenti, è sempre possibile rappresentare J! u( ) nel modo seguente (vedi

Appendice 1)

J! u( ) = A!"!A

!

#1 (5.150) dove

!" = diag #1,#

2,...,#n( ) (5.151)

e

A! = a1,a2 ,...,an . (5.152)

La matrice

A! è invertibile perché gli autovettori di

J! u( ) sono linearmente

indipendenti. Sostituendo la (5.150) nella (5.148) abbiamo

!ek( ) ! "# u( ) !e k$1( ) (5.153)

dove

!ek( )= A!

"1ek( ) . (5.154)

Essendo !" diagonale si ha che

!ei

k( )

!ei

k!1( )" #

iJ$( ) per i = 1,2,...,n . (5.155)

In conclusione, la velocità di convergenza del metodo di Picard è di ordine uno. 5.4.1.3 Stima dell’errore

Come nel caso scalare, cerchiamo una relazione tra la norma dell’errore al passo k , ek( ) , e la norma dell’incremento al passo k +1 , !x

k+1( ) nell’intorno del punto fisso. Essendo

41


x

k( )! u = x

k+1( )! u( ) ! x

k+1( )! x

k( )( ) , (5.156) dalle (5.121) e (5.122) abbiamo che

ek( )= ! x k( )( ) " ! u( )#$

%& " 'x k+1( ) , (5.157)

da cui, applicando la (5.145) si ottiene

ek( ) ! J" u( )e k( )

+1

2ek( )T

H"ek( ) # $x

k+1( ) . (5.158)

Siccome e k( )

! 0 per k!" , per k sufficientemente grande possiamo trascurare nella (5.158) i termini quadratici, quindi si ha

ek( ) ! J" u( )e k( ) # $x

k+1( ) . (5.159) Allora, per k sufficientemente grande si ha

ek( ) !

1

1" J# u( )$x

k+1( ) (5.160)

(abbiamo usato le seguenti proprietà di una norma ! : a + b ! a + b e

La ! L a ,

vedi Appendice 2). Questa relazione ci consente di stimare il limite superiore per l’errore eN( ) per un dato !x

N +1( ) . 5.4.2 Iterazione di Newton-Raphson

La scelta più generale di funzione di iterazione per risolvere l’equazione (5.4) è data da

! x( ) = x + L x( ) f x( ) (5.161)

dove, ora,

L = L x( ) è una matrice dipendente da x non singolare, da scegliere in modo

tale che J! x( ) < 1 almeno in un intorno della soluzione u . Determiniamo la matrice

jacobiana di ! . Abbiamo

!"!x

= I +!

!xL x( ) f x( )#$ %& = I + L x( ) J f x( ) + B x( ) f x( ){ } (5.162)

42


dove

B x( ) f x( ){ } è la matrice n ! n con elementi

bhkf{ } dati da

bhkf{ } = a

hk

T f , (5.163)

ahk =!lh1

!xk,!lh2

!xk,...,

!lhp

!xk

T

, (5.164)

e l

hk sono gli elementi di L . Analogamente al caso scalare conviene scegliere

L = !J f

!1 (5.165) purché

J f sia invertibile. Con questa scelta l’iterazione di Picard si riduce proprio

all’iterazione di Newton-Raphson e

J! x( ) = B x( ) f x( ){ } . (5.166)

E’ immediato allora che

J! u( ) = B u( ) f u( ){ } = 0 . Consideriamo la norma del modulo di

J! x( ) ,

J! x( )

1

. Abbiamo

J! x( )1

= maxk

ahk

Tf

h=1

N

" # f2max

k

ahk 2

h=1

N

" . (5.167)

Allora, esiste un intorno J del punto fisso u nel quale

J! x( )

1

< 1 . Scegliendo x0( )!J

l’iterazione di Newton-Raphson converge alla soluzione u . Esempio 5.4

Nell’Esempio 5.3 abbiamo visto che l’iterazione di Newton-Raphson non converge quando si risolve il sistema (5.93) con questi valori dei parametri R

1= 500! ,

R2= 200k! , E

1= 0.5V , E

2= 10V , !

R= 0.8 , !

F= 0.9 , IES = ICS = 1pA , V

T= 26mV e

con x 0( )= 0 . Alla terza iterazione la soluzione è “NaN”.

Noi sappiamo che il metodo di Newton-Raphson converge se il dato iniziale x 0( ) è sufficientemente vicino alla soluzione. Possiamo, allora, applicare il metodo della continuazione per risolvere questo problema. Il sistema di equazioni in esame ha la soluzione w

b= w

c= 0 per E

1= E

2= 0 .

Si consideri la soluzione del sistema (5.93) per

43


E1= h!E

1, E

2= h!E

2 con h = 1,2,...,M (5.168)

dove M , !E

1 e !E

2 sono scelti in modo tale che M!E

1= 0.5V e M!E

2= 10V . Il

problema h ! esimo, per h = 1,2,...,M , può essere risolto con il metodo di Newton-Raphson assumendo come dato iniziale x

h

0( ) la soluzione del problema h !1( ) ! esimo, uh!1

. Se !E1 e !E

2 sono sufficientemente piccoli l’iterazione di Newton-Raphson con

xh

0( )= u

h!1 converge alla soluzione u

h per h = 1,2,...,M .

♦ Esercizio 5.17

Risolvere con il Metodo di Newton-Raphson il sistema (5.93) per R1= 500! ,

R2= 200k! , E

1= 0.5V , E

2= 10V , !

R= 0.8 , !

F= 0.9 , IES = ICS = 1pA , V

T= 26mV

applicando la procedura della continuazione descritta nell’Esempio 5.4. ♦

Ora mostreremo che l’ordine di convergenza è uguale a 2. Siccome nel caso in esame

J! u( ) = B u( ) f u( ){ } = 0 , dalla (5.147) si ha che

ek( ) !

1

2ek"1( )T

H# u( )e k"1( ) . (5.169)

Essendo

ek!1( )T

H" u( )e k!1( )= e

k!1( )TH1u( )e k!1( )

,ek!1( )T

H2u( )e k!1( )

,...,ek!1( )T

Hnu( )e k!1( ) T (5.170)

dove

H1,H

2,...,H

n sono gli le matrice hessiane di !

1,!

2,...,!

n, rispettivamente, si ha che

(vedi Appendice 2)

ek!1( )T

H" u( )e k!1( )

1

= ek!1( )T

Hiu( )e k!1( )

i=1

n

# $ e k!1( )

2

2

Hiu( )

2

i=1

n

#

$ e k!1( )

1

2

Hiu( )

2

i=1

n

# .

(5.171)

(Ricordiamo che !1 è la norma del modulo, x

1= x

i

i=1

n

! .) Combinando la (5.169) e la

(5.171) si ha

44


ek( )

1

!1

2ek"1( )

1

2

Hiu( )

2

i=1

n

# . (5.172)

Per l’equivalenza tra le norme esistono due coefficienti positivi c

ms( ) e c

Ms( ) tali che

c

ms( ) e k( )

s

! ek( )

1

(5.173)

e c

Ms( ) e k( )

s

! ek( )

1

(5.174)

per ogni s . Combinando le (5.172)-(5.174) abbiamo

ek( )

p!cM2p( )

2cm p( )Hi u( )

2

i=1

n

" ek#1( )

p

2

, (5.175)

quindi

!k =ek( )

p

ek"1( )

p

#cM2p( )

2cm p( )Hi u( )

2

i=1

n

$%

&'

(

)* e

k"1( )

p. (5.176)

Pertanto, l’ordine di convergenza del metodo di Newton-Raphson è almeno uguale a 2. La costante asintotica di errore è da

C =cM2p( )

2cm p( )Hi u( )

2

i=1

n

! . (5.177)

Referenze

[5.1] J.M. Ortega, W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables, Academic Press, New York, 1969. [5.2] A. Quarteroni, F. Saleri, Introduzione al Calcolo Scientifico, Springer, 2002. [5.3] V. Comincioli, Analisi numerica: metodi ed applicazioni, McGraw-Hill, Milano 1990.

5. Equazioni algebriche non lineari - Elettrotecnica · Solo in qualche caso è possibile risolvere...

Documents

Transcript of 5. Equazioni algebriche non lineari - Elettrotecnica · Solo in qualche caso è possibile risolvere...