Metodi iterativi per la soluzione di sistemi nonlineari

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica

28 aprile 2014

Introduzione

Il problema della soluzione di sistemi nonlineari è di importanza fondamentale in vari aspetti della modellistica matematica come ad esempio nell’ambito dell’integrazione di equazioni differenziali nonlineari con metodi impliciti.

Data una funzione F : ℝ^{n} → R^{n} si tratta di trovare gli x^{∗} tali che F (x^{∗}) = 0 (cf. [6, p.254]).

Definizione

Se d è una distanza, e (X , d) uno spazio metrico (come ad esempio (ℝ, d) con d(x, y) = |x − y|), M un sottospazio non vuoto di X allora F : M → M è L-contrattiva in M se e solo se L < 1 e

d(F (x), F (y)) ≤ Ld(x, y), ∀x, y ∈ M.

Teorema di Banach (o di punto fisso)

Banach provò nel 1922 il seguente teorema (cf. [3, p.432]).

Teorema

Sia (X , d) uno spazio metrico completo ed M un sottospazio non vuoto e chiuso di X . Se T : M → M è L-contrattiva allora

1. l’equazione x = T (x) ha una ed una sola soluzione α;

2. la successione x_{k+1} = T (x_{k}) (detta di Banach o di punto fisso) converge alla soluzione α per ogni scelta del punto iniziale x_{0};

3. per k = 0, 1, 2, … si ha

d(x_{k}, α) ≤ L^{k}(1 − L)^{−1}d(x_{0}, x_{1}), a priori

d(x_{k+1}, α) ≤ L(1 − L)^{−1}d(x_{k+1}, x_{k}), a posteriori

4. per k = 0, 1, 2, … si ha

d(x_{k+1}, α) ≤ Ld(x_{k}, α)
Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  3/ 41

Si vede inoltre che

Teorema

Condizione sufficiente affinchè T sia una contrazione ∥ · ∥_{∞} è che sia di classe C^{1}(K), con K chiusura di un aperto convesso e

sup_{x∈K} ∥JT (x)∥_{∞} < 1

dove JT è la matrice jacobiana di T .

\underline{Sketch della dimostrazione}: applicazione del teorema del valor medio in più variabili.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  4/ 41

Teorema (Convergenza locale)

Se T è di classe C^{1}(Ω) con Ω intorno del punto fisso x^{∗} e sup_{x∈Ω} ∥JT (x)∥_{∞} < 1, allora il metodo delle approssimazioni successive converge localmente a x^{∗}.

\underline{Sketch della dimostrazione}: si prenda come K una opportuna palla chiusa centrata in x^{∗} dove sup_{x∈K} ∥JT (x)∥_{∞} < 1.

Teorema (Convergenza locale)

Sia B_{∞}[x_{0}, r] la palla chiusa di centro x_{0} e raggio r in ∥ · ∥_{∞}. Se θ = max_{x∈B_{∞}[x_{0},r]} ∥JT (x)∥_{∞} < 1, allora il metodo delle approssimazioni successive converge quando ∥x_{1} − x_{0}∥ ≤ r(1 − θ).

\underline{Sketch della dimostrazione}: si prenda K = B_{∞}[x_{0}, r] e si verifichi che T (K) ⊆ K.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  5/ 41

Per quanto riguarda la velocità di convergenza del metodo di punto fisso, nelle ipotesi del teorema delle contrazioni, si ha che

▶ la convergenza è almeno lineare visto che

∥x^{∗} − x_{n+1}∥ ≤ L∥x^{∗} − x_{n}∥

con L < 1. ▶ se T ∈ C^{2}(Ω) con Ω intorno del punto fisso x^{∗} e JT (x^{∗}) = 0,

la convergenza diventa localmente almeno quadratica, ovvero

∥x^{∗} − x_{n+1}∥ ≤ c∥x^{∗} − x_{n}∥^{2}

con una opportuna costante c per n abbastanza grande.

Esempio

Supponiamo di dover risolvere il sistema nonlineare (cf. [3, p.449])

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}  (1)

Chiaramente è un sistema nonlineare avente due equazioni e due incognite. Per prima cosa ci si chiede quante soluzioni abbia questo problema.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  7/ 41

Esempio

Dal punto di vista pratico, sostituendo

y = \frac{1}{2} sin (x)

nella prima equazione, otteniamo

x = \frac{1}{2} cos (\frac{1}{2} sin (x)).

La funzione

g(x) = \frac{1}{2} cos (\frac{1}{2} sin (x))

è tale che

g^{(1)}(x) = (−1/4) sin ((1/2) · sin (x)) cos (x).

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  8/ 41

Esempio

Per convincersene, usando il calcolo simbolico in Matlab

>> syms x

>> g =(1/2)∗ cos ( ( 1 / 2 ) ∗ s i n ( x ) )

 g =

1/2∗ cos (1/2∗ s i n ( x ) )

>>

>> % DERIVATA PRIMA .

>>

>> d i f f ( g )
ans =
−1/4∗ s i n (1/2∗ s i n ( x ) ) ∗ cos ( x )

>>

>> % DERIVATA SECONDA.

>>

>> d i f f ( g , 2 )
ans =
−1/8∗ cos (1/2∗ s i n ( x ) ) ∗ cos ( x ) ˆ 2+1/4∗ s i n (1/2∗ s i n ( x ) ) ∗ s i n ( x )

>>

ed essendo

0 ≤ | sin (\frac{sin(x)}{2}) · cos (x)| ≤ 1

si ha che

|g^{(1)}(x)| = |(−1/4) sin ((1/2) · sin (x)) cos (x)| ≤ \frac{1}{4}, ∀x ∈ ℝ.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  9/ 41

Esempio

Dalla formula di Taylor, centrata in un arbitrario y

g(x) = g(y) + g^{(1)}(ξ)(x − y)

e quindi portando g(y) a primo membro e calcolando il modulo ad ambo i membri

|g(x) − g(y)| = |g^{(1)}(ξ)||(x − y)| ≤ \frac{1}{4} |x − y|.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  10/ 41

Esempio

Per quanto detto, se d è una distanza, e (X , d) uno spazio metrico (come ad esempio (ℝ, d) con d(x, y) = |x − y|), M un sottospazio non vuoto di X allora F : M → M è L-contrattiva in M se e solo se L < 1 e

d(F (x) − F (y)) ≤ Ld(x, y), ∀x, y ∈ M.

Per quanto visto, da

|g(x) − g(x_{0})| = |g^{(1)}(ξ)||(x − x_{0})| ≤ \frac{1}{4} |(x − x_{0})|

abbiamo provato che

g(x) = \frac{1}{2} cos (\frac{1}{2} sin (x))

è L-contrattiva in M = ℝ per L = 1/4.

Teorema di Banach

Dal teorema di Banach, in virtù della L-contrattività , deduciamo che

g(x) = \frac{1}{2} cos (\frac{1}{2} sin (x))

ha una ed una sola soluzione e quindi pure

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}  (2)
Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  12/ 41

Punto fisso in Matlab

Calcoliamo la successione di punto fisso in Matlab. Scriviamo il

codice punto fisso

 f u n c t i o n xhist=punto_fisso ( x0 , maxit , tol , g )

% INPUT .
% x0  : PUNTO INIZIALE .
% maxit  : NUMERO MASSIMO DI ITERAZIONI .
% t o l  : TOLLERANZA CRITERIO DI ARRESTO.
% g  : EQUAZIONE DA RISOLVERE x=g ( x )

x_old=x0 ;
xhist =[ x_old ] ;
 f o r index =1: maxit
     x_new=f e v a l ( g , x_old ) ;
      i f norm ( x_new−x_old , inf ) < tol

r e t u r n e l s e

        xhist =[ xhist ; x_new ] ;
        x_old=x_new ;

end end

f p r i n t f ( ’\n \t [WARNING ] : RAGGIUNTO MASSIMO NUMERO DI ITERAZIONI ’ ) ;

Applichiamo la successione descritta nel teorema di Banach per risolvere il problema x = g(x),

>> g=inline ( ’ (1/2) ∗ cos ( ( 1 / 2 ) ∗ s i n ( x ) ) ’ ) ;
>> xhist=punto_fisso ( 0 , 50 , 10ˆ ( −15) , g ) ;
>> format long

>> xhist
xhist =

0

0.50000000000000 0.48570310513698 0.48644107261515 0.48640331614624 0.48640524877124 0.48640514984910 0.48640515491247 0.48640515465330 0.48640515466657 0.48640515466589 0.48640515466592 0.48640515466592

>> f e v a l ( g , 0 . 4 8 6 4 0 5 1 5 4 6 6 5 9 2 )
ans =

0.48640515466592 >>

Nelle ultime due righe dell’esperimento abbiamo verificato che effettivamente per t = 0.48640515466592, t ≈ g(t).

Punto fisso in Matlab

Vediamo se anche il sistema nonlineare

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}  (3)

è risolto correttamente.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  15/ 41

Punto fisso in Matlab

▶ Se x = 0.48640515466592, allora dalla seconda equazione

y = \frac{1}{2} sin (0.48640515466592) ≈ 0.23372550195872.

▶ Verifichiamo che anche la prima equazione sia risolta. In

effetti si ha x = 0.48640515466592 e

\frac{1}{2} sin (x) ≈ 0.23372550195872

che coincide proprio con y.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  16/ 41

In Matlab/Octave

>> format long ;
>> x =0.48640515466592;
>> y=0.5∗ s i n ( x ) ;
>> y

y =

0.23372550195872

>> 0.5∗ cos ( y )
ans =

0.48640515466592 >>

Dal punto di vista pratico abbiamo risolto un sistema nonlineare in due variabili e due incognite come un’equazione lineare in un’incognita. Evidentemente ciò non è sempre possibile.

Punto fisso in Matlab

Riproponiamo quanto visto nell’esempio precedente

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}

in modo differente.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  18/ 41

Siano v = (v_{i})_{i} = [x; y] e

T_{1}(v) = \frac{1}{2} cos (v_{2})  (4)
T_{2}(v) = \frac{1}{2} sin (v_{1}).

Per T (v) = [T_{1}(v); T_{2}(v)] abbiamo v = T (v). Mostriamo che T : ℝ^{2} → ℝ^{2} è contrattiva.

Dalla formula di Taylor in più variabili (cf. [4, p.192]), T (v) = T (v_{0}) + (T′(ξ)) · (v − v_{0}), ξ ∈ S

in cui S il segmento che congiunge v con v_{0} e T′ la matrice Jacobiana definita per componenti come

(T′)_{j,k}(v) = \frac{∂T_{j}(v)}{∂v_{k}}.

Punto fisso in Matlab

Nel nostro esempio quindi, T_{1}(v) = \frac{1}{2} cos (v_{2}), T_{2}(v) = \frac{1}{2} sin (v_{1}),

T′(v) = ( \begin{matrix}0 & (−1/2) sin (v_{1}) \\ (+1/2) cos (v_{2}) & 0 \end{matrix} )

Essendo la norma infinito di una matrice A definita da

∥A∥_{∞} = max_{j} ∑_{i}|a_{i,j}|

si ha

∥(T′)(v)∥_{∞} = max_{j} ∑_{i}| ((T′)(v))_{i,j} |
= max (|(−1/2) sin (v_{1})|, |(+1/2) cos (v_{2})|) ≤ \frac{1}{2}
Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  20/ 41

Dalla formula di Taylor si ha così che

T (v) = T (v_{0}) + (T′(ξ)) · (v − v_{0}), ξ ∈ S

e calcolando la norma infinito di ambo i membri

∥T (v) − T (v_{0})∥_{∞} ≤ ∥(T′(ξ))∥_{∞}∥(v − v_{0})∥_{∞}, ξ ∈ S

da cui T : ℝ^{2} → ℝ^{2} è una L-contrazione con L = 1/2.

Possiamo nuovamente applicare il metodo di Banach (detto anche di punto fisso), questa volta però non per risolvere l’equazione x = g(x), bensì direttamente il sistema nonlineare v = T (v).

Punto fisso in Matlab

Scriviamo il file g1

f u n c t i o n res=g1 ( v )
res ( 1 , 1 ) =0.5∗ cos ( v ( 2 ) ) ;
res ( 2 , 1 ) =0.5∗ s i n ( v ( 1 ) ) ;

e quindi dalla shell di Matlab/Octave

>> format long ;
>> x0 =[0; 0 ] ;
>> maxit =50;
>> tol=10ˆ (−15) ;

>> xhist=punto_fisso ( x0 , maxit , tol , ’ g1 ’ )
xhist =

0  0
0.50000000000000  0
0.50000000000000  0.23971276930210
0.48570310513698  0.23971276930210
0.48570310513698  0.23341513183103
0.48644107261515  0.23341513183103
0.48644107261515  0.23374137788239
0.48640331614624  0.23374137788239
0.48640331614624  0.23372468931518
0.48640524877124  0.23372468931518

…

0.48640515466657  0.23372550195901
0.48640515466589  0.23372550195901
0.48640515466589  0.23372550195871
0.48640515466592  0.23372550195871
0.48640515466592  0.23372550195872
0.48640515466592  0.23372550195872

>>

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  22/ 41

Vediamo quanto accurato è il teorema di Banach nelle sue stime. Il metodo esegue 24 iterazioni. Calcoliamo

1. Gli errori assoluti compiuti, posto

α = [0.48640515466592 0.23372550195872].

2. Essendo L = 1/2 dal teorema di Banach si vede che la stima a posteriori

d(x_{k+1}, α) ≤ L(1 − L)^{−1}d(x_{k+1}, x_{k}),

afferma

d(x_{k+1}, α) ≤ d(x_{k+1}, x_{k}), k = 0, 1, …

Punto fisso in Matlab

abbiamo

>> format short e
>> % CALCOLIAMO ERRORI ASSOLUTI .
>> xhist_err=abs ( xhist−repmat ( [ 0 . 4 8 6 4 0 5 1 5 4 6 6 5 9 2  0 . 2 3 3 7 2 5 5 0 1 9 5 8 7 2 ] , 2 4 , 1 ) ) ;
>> xhist_abserr=s q r t ( ( xhist_err ( : , 1 ) ) .ˆ 2+( xhist_err ( : , 2 ) ) . ˆ 2 ) ;

>>

>> % CALCOLIAMO STIMA A POSTERIORI .
>> N=24;

>> xplus=xhist ( 2 : N , : ) ; xminus=xhist ( 1 : N −1 ,:) ;
>> xdiff=xplus−xminus ;

>> xerr_est=s q r t ( ( xdiff ( : , 1 ) ) . ˆ 2 + ( xdiff ( : , 2 ) ) . ˆ 2 )

>>

>> % PARAGONIAMO GLI ERRORI A PARTIRE DA k=1. OSSERVIAMO CHE
>> % GLI INDICI IN MATLAB PARTONO DA 1 .

>> [ xhist_abserr ( 2 : N , : ) xerr_est ]

ans =
2.3412 e−001 5.0000 e−001
1.4855 e−002 2.3971 e−001
6.0283 e−003 1.4297 e−002
7.6760 e−004 6.2976 e−003
3.1244 e−004 7.3797 e−004
3.9270 e−005 3.2625 e−004
1.5982 e−005 3.7756 e−005

…

2.8805 e−013 6.7901 e−013
3.4630 e−014 3.0012 e−013
1.4144 e−014 3.4750 e−014
3.3766 e−015 1.5377 e−014
1.9767 e−015 1.7764 e−015

>>

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  24/ 41

▶ Si osservi che, eccetto nell’ultima riga, la stima dall’alto è

effettivamente realizzata. ▶ Nell’ultima riga, abbiamo approssimato la soluzione α a 14

cifre dopo la virgola, ed è la ragione dell’anomalia numerica.

Nota.

▶ in generale un’equazione si scrive come F (x) = 0, mentre nella formulazione di punto fisso si descrive come x = T (x); evidentemente non è un problema, in quanto se F (x) = 0 allora x = T (x) per T (x) := F (x) + x;

▶ l’esempio mostrato nella sezione è particolarmente semplice, in quanto può essere ricondotto ad un problema unidimensionale;

ciò nonostante, abbiamo mostrato che può essere risolto da un metodo che lavora in più variabili.

Metodo di Newton

Supponiamo di dover risolvere il problema F (v) = 0. Dalla formula di Taylor (multivariata), se la funzione F : M ⊆ ℝ^{n} → ℝ^{n} è sufficientemente regolare, allora per v_{k+1} sufficientemente vicino a v_{k}

F (v_{k+1}) ≈ F (v_{k}) + (F′(v_{k})) · (v_{k+1} − v_{k}) .  (5)

Ricordiamo che F′(v_{k}) è una matrice, e il successivo prodotto · è di tipo matrice per vettore.

Da (5), supposto che la matrice F′(v_{0}) sia invertibile, e che F (v_{k+1}) ≈ 0, cioè v_{k+1} sia una approssimazione di uno zero di F , abbiamo

0 ≈ F (v_{k+1}) ≈ F (v_{k}) + (F′(v_{k})) · (v_{k+1} − v_{k})  (6)

da cui ha senso definire la successione (del metodo di Newton)

v_{k+1} := v_{k} − (F′(v_{k}))^{−1} · F (v_{k}).  (7)

Talvolta, la successione di Newton viene descritta come

\{ \begin{matrix}F′(v_{k})) · h_{k+1} = −F (v_{k}), \\ v_{k+1} := v_{k} + h_{k+1} \end{matrix}  (8)

sottolineando che l’inversa di F′(v_{k})) non deve necessariamente essere calcolata, bensì bisogna risolvere un sistema lineare la cui soluzione h_{k+1} è la correzione da apportare a v_{k} per ottenere v_{k+1} (cf. [5, p.174] per un esempio in ℝ^{2}).

L’analisi di convergenza del metodo di Newton in dimensione maggiore di 1 e più in generale in spazi di Banach è un attivo e difficile argomento di ricerca. A tal proposito si consulti [2, p.154].

Per quanto riguarda la velocità di convergenza del metodo di Newton

Teorema

Se f ∈ C^{2}(K) con K chiusura di un aperto convesso e limitato contenente x^{∗}, in cui la Jacobiana di f è invertibile, e supposto che le iterazioni x_{n} siano tutte in K, posto e_{n} = ∥x^{∗} − x_{n}∥_{2} vale la seguente stima (convergenza almeno quadratica)

e_{n+1} ≤ ce_{n}^{2}, n ≥ 0

con

c = \frac{\sqrt{m}}{2} max_{x∈K}∥(Jf (x))^{−1}∥_{2} max_{1≤i≤m} max_{x∈K}∥Hf_{i}(x)∥_{2}

in cui H è la matrice Hessiana cioè (Hf )_{i,j} = \frac{∂^{2}f}{∂x_{i}∂x_{j}} , e Hf_{i} la matrice hessiana della componente f_{i}.

Per quanto riguarda la convergenza locale del metodo di Newton

Teorema

Se f ∈ C^{2}(K) e Jf (x) è invertibile in K = B_{2}[x^{∗}, r] (dove x^{∗} è la soluzione del sistema, f (x) = 0), per

c = \frac{\sqrt{m}}{2} max_{x∈K}∥(Jf (x))^{−1}∥_{2} max_{1≤i≤m} max_{x∈K}∥Hf_{i}(x)∥_{2}

scelto x_{0} tale che e_{0} < min{1/c, r}, il metodo di Newton è convergente e vale la seguente stima dell?errore

ce_{n} ≤ (c e_{0})^{2^{n}}, n ≥ 0.

\underline{Sketch della dimostrazione}: per induzione e_{n+1} ≤ (ce_{n})e_{n} < e_{n} e quindi x_{n+1} ∈ B_{2}[x^{∗}, r].

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  31/ 41

Nelle ipotesi di convergenza locale, la stima a posteriori dell’errore con lo step ∥x_{n+1} − x_{n}∥ è una buona stima (almeno per n abbastanza grande), cioè

e_{n} = ∥x^{∗} − x_{n}∥ ≈ ∥x_{n+1} − x_{n}∥.

\underline{Sketch della dimostrazione}: si osservi che f è localmente invertibile e Jf^{−1}(f (x)) = (Jf (f (x)))^{−1}, da cui

x^{∗} − x_{n} = f^{−1}(f (x^{∗})) − f^{−1}(f (x_{n})) ≈ Jf^{−1}(f (x_{n}))(f (x^{∗}) − f (x_{n})).

Nota.

Il metodo di Newton corrisponde ad un’iterazione di punto fisso con

ϕ(x) = x − (Jf (x))^{−1}f (x),

da cui si deduce che se f ∈ C^{2}(Ω), con Ω intorno di x^{∗} allora la convergenza è localmente almeno quadratica perchè Jϕ(x^{∗}) = 0.

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  32/ 41

Metodo di Newton in Matlab

Consideriamo l’esempio precedente, cioè

\{ \begin{matrix}x − \frac{1}{2} cos (y) = 0 \\ y − \frac{1}{2} sin (x) = 0 \end{matrix}  (9)

e implementiamo il metodo di Newton

 f u n c t i o n [ v_hist , step_hist , residual_hist ]= newton_sistemi ( v0 , maxit , tol )

% INPUT :
% v0  : APPROSSIMAZIONE INIZIALE .
% maxit  : NUMERO MASSIMO DI ITERAZIONI .
% t o l  : TOLLERANZA DEL CRITERIO D’ARRESTO.
% F  : FUNZIONE .
% DF  : FUNZIONE DA CUI VIENE CALCOLATA LA JACOBIANA DI F .
% OUTPUT:
% v h i s t : VETTORE CONTENENTE LE APPROSSIMAZIONI DELLA SOLUZIONE .
% s t e p h i s t  : VETTORE DEGLI STEP .
% r e s i d u a l h i s t : VETTORE DEI RESIDUI .

Metodo di Newton in Matlab

v_old=v0 ;
v_hist =[v0 ’ ] ;

JF=DF ( v_old ) ; % CALCOLIAMO LA MATRICE JACOBIANA .
Fv=F ( v_old ) ; % VALUTIAMO LA FUNZIONE .

residual=norm ( Fv ) ;
residual_hist =[ residual ] ;

h=−JF\Fv ;  % CALCOLIAMO LA CORREZIONE .
v_new=v_old+h ; % NUOVA APPROSSIMAZIONE .

v_hist =[ v_hist ; v_new ’ ] ;
step_hist = [ ] ;
f o r index =1: maxit

Fv=F ( v_new ) ; % VALUTIAMO LA FUNZIONE .
residual=norm ( Fv , 2 ) ;  % RESIDUO .

     residual_hist =[ residual_hist ; residual ] ;
     loc_step=norm ( v_new−v_old , 2 ) ;  % STEP .
     step_hist =[ step_hist ; loc_step ] ;
     i f max( loc_step , residual ) < tol

r e t u r n ;

e l s e

          v_old=v_new ;

JF=DF ( v_old ) ; % CALCOLIAMO LA MATRICE JACOBIANA .
h=−JF\Fv ;  % CALCOLIAMO LA CORREZIONE .
v_new=v_old+h ; % NUOVA APPROSSIMAZIONE .

          v_hist =[ v_hist ; v_new ’ ] ;

end end

f p r i n t f ( ’\n \t [WARNING ] : NUMERO MASSIMO DI ITERAZIONI RAGGIUNTO . ’ )

Metodo di Newton in Matlab

f u n c t i o n res=F ( v )
res ( 1 , 1 ) =0.5∗ cos ( v ( 2 ) )−v ( 1 ) ;
res ( 2 , 1 ) =0.5∗ s i n ( v ( 1 ) )−v ( 2 ) ;
f u n c t i o n res=DF ( v )
% r e s ( 1 , 1 ) =0.5∗ cos ( v ( 2 ) ) ;
% r e s ( 2 , 1 ) =0.5∗ s i n ( v ( 1 ) ) ;
res=[−1 −0.5∗ s i n ( v ( 2 ) ) ; 0.5∗ cos ( v ( 1 ) ) −1];

Quindi scriviamo il programma principale driver newton che setta i parametri della procedura newton sistemi

c l e a r ;
v0 =[0; 0 ] ; maxit =50; tol=10ˆ (−10) ;

[ v_hist , step_hist , residual_hist ]= newton_sistemi ( v0 , maxit , tol ) ;

Metodo di Newton

Otteniamo così

>> driver_newton
>> v_hist
v_hist =

0  0
5.0000 e−001 2.5000 e−001
4.8646 e−001 2.3377 e−001
4.8641 e−001 2.3373 e−001
4.8641 e−001 2.3373 e−001
4.8641 e−001 2.3373 e−001

>> step_hist
step_hist =

5.5902 e−001
2.1132 e−002
7.5293 e−005
7.7616 e−010

0

>> residual_hist
residual_hist =

5.0000 e−001
1.8640 e−002
6.7480 e−005
6.7927 e−010

0 0

>>

▶ Dal vettore degli step |x_{k+1} − x_{k}| e da quello dei residui |F (x_{k})| si capisce che il metodo ha convergenza superlineare (in realtà è quadratica).

▶ Se consideriamo quale soluzione il risultato fornito dall’ultima iterazione (il residuo è 0!), possiamo calcolare gli errori assoluti ad ogni iterazione, convincendoci una volta di più della convergenza superlineare del metodo.

>> alpha=v_hist ( s i z e ( v_hist , 1 ) , : ) ;
>> err_vett=v_hist−repmat ( alpha , s i z e ( v_hist , 1 ) , 1 ) ;

>> abs_err_vett=s q r t ( ( err_vett ( : , 1 ) ) . ˆ 2 + ( err_vett ( : , 2 ) ) . ˆ 2 )

 abs_err_vett =

5.3965 e−001
2.1206 e−002
7.5294 e−005
7.7616 e−010

0 0

>>

Sul comando Matlab repmat

Osservazione. Il comando Matlab repmat è di uso molto comune. Nel nostro caso volevamo costruire, alla riga

err_vett=v_hist−repmat ( alpha , s i z e ( v_hist , 1 ) , 1 ) ;

la differenza tra la matrice v hist e la matrice in cui ogni riga è costituita dal vettore riga soluzione (α_{1} α_{2}). Per avere le dimensioni corrette per poter eseguire la differenza abbiamo dovuto replicare il vettore riga α tante volte quanto il numero di righe di v hist, cosa eseguita appunto dal comando repmat. Vediamo un esempio di repmat.

>> v=[1 2 ] ;
>> repmat ( v , 3 , 1 )
ans =
1  2
1  2
1  2
>> repmat ( v , 3 , 2 )
ans =
1  2  1  2
1  2  1  2
1  2  1  2

>>

Alvise Sommariva  Metodi iter. per la soluzione di sistemi nonlineari  39/ 41

Sul comando Matlab repmat

In pratica, repmat(v,m,n) costruisce una matrice m × n in cui ogni componente è uguale al vettore v, come si può evincere da

>> h e l p repmat
REPMAT  Replicate  and  tile an  array .
B = repmat ( A , M , N ) creates a large  matrix B consisting  of an M−by−N
tiling  of  copies  of A .
B = REPMAT ( A , [ M N ] ) accomplishes  the  same  result  as  repmat ( A , M , N ) .
B = REPMAT ( A , [ M N P . . . ] ) tiles  the  array A to  produce a
M−by−N−by−P−by −. . . block  array . A can be N−D .
REPMAT ( A , M , N ) when A is a scalar  is  commonly  used to  produce
an M−by−N matrix  filled  with A ’ s value . This  can be much  faster
than A∗ONES ( M , N ) when M and /or N are  large .
Example :

repmat ( magic ( 2 ) , 2 , 3 )
repmat (NaN, 2 , 3 )

See  also  MESHGRID .

>>

Bibliografia

K. Atkinson, An Introduction to Numerical Analysis, Wiley, (1989). K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer Verlag, (2001). V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990. G. Gilardi, Analisi due, seconda edizione, Mc Graw-Hill, 1996. J.H. Mathews e K.D. Fink, Numerical methods using Matlab, terza edizione, Prentice-Hall, 1999. A. Quarteroni, F. Saleri Introduzione al Calcolo Scientifico, Springer, (2002). The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler. G. Rodriguez, Algoritmi numerici, Pitagora editrice, 2008.