Sistemi di equazioni nonlineari

1 giugno 2008

1 Un esempio introduttivo sul metodo di punto fisso

Il problema della soluzione di sistemi nonlineari è di importanza fondamentale in vari aspetti della modellistica matematica. Ad esempio è fondamentale nella risoluzione di equazioni differenziali nonlineari con metodi impliciti.

Data una funzione : �^{}^{}si tratta di trovare gli ^{}tali che ( ^{}) 0 (cf. [6, p.254]). Vediamone un esempio. Supponiamo di dover risolvere il sistema nonlineare (cf. [3, p.449])

\begin{matrix}\frac{1}{2} cos( ) 0 \\ \frac{1}{2} sin( ) 0 \end{matrix}  (1)

Per prima cosa ci si chiede quante soluzioni abbia questo problema. Dal punto di vista pratico, sostituendo

\frac{1}{2} sin( )

nella prima equazione, otteniamo

\frac{1}{2} cos \frac{1}{2} sin( )

La funzione

( ) \frac{1}{2} cos \frac{1}{2} sin( )

è tale che (1)( ) ( 1/4)sin((1/2) sin( ))cos( )

Per convincersene, usando il calcolo simbolico in Matlab

>> syms x
>> g=(1/2)*cos((1/2)*sin(x))
g =
1/2*cos(1/2*sin(x))
>>
>> % DERIVATA PRIMA.
>>
>> diff(g)

Sistemi nonlineari

ans =
-1/4*sin(1/2*sin(x))*cos(x)
>>
>> % DERIVATA SECONDA.
>>
>> diff(g,2)
ans =
-1/8*cos(1/2*sin(x))*cos(x)^2+1/4*sin(1/2*sin(x))*sin(x)
>>

ed essendo

0 sin \frac{( )}{2} cos( ) 1

si ha che

(1)( ) \frac{1}{4} �

Dalla formula di Taylor, centrata in un arbitrario _{0}

( ) ( _{0}) ^{(1)}( )( _{0})

e quindi portando ( _{0}) a primo membro e calcolando il modulo ad ambo i membri

( ) ( _{0}) ^{(1)}( ) ( _{0}) \frac{1}{4} ( _{0})

Ricordiamo ora che se è una distanza, e ( ) uno spazio metrico (come ad esempio (� ) con ( ) ), un sottospazio non vuoto di allora : è contrattiva in se e solo se

( ) ( )

Per quanto visto abbiamo provato che

( ) \frac{1}{2} cos \frac{1}{2} sin( )

è contrattiva in � per 1/4. Banach provò nel 1922 il seguente (cf. [3, p.432])

Teorema 1.1 Sia ( ) uno spazio metrico completo ed un sottospazio non vuoto e chiuso di . Se : è contrattiva allora

1. l’equazione ( ) ha una ed una sola soluzione ;

2. la successione _{1} ( _{}) (detta di Banach o di punto fisso) converge alla soluzione per ogni scelta del punto iniziale _{0};

3. per 0 1 2 si ha

( _{}) ^{}(1 )^{1} ( _{0} _{1}) a priori
( _{1} ) (1 )^{1} ( _{1} _{}) a posteriori

4. per 0 1 2 si ha

( _{1} ) ( _{})

Dal teorema di Banach deduciamo che

( ) \frac{1}{2} cos \frac{1}{2} sin( )

ha una ed una sola soluzione e quindi pure

\begin{matrix}\frac{1}{2} cos( ) 0 \\ \frac{1}{2} sin( ) 0 \end{matrix}  (2)

Vediamo i dettagli in Matlab. Scriviamo il codice punto_fisso

function xhist=punto_fisso(x0, maxit, tol, g)
% INPUT.
% x0      : PUNTO INIZIALE.
% maxit   : NUMERO MASSIMO DI ITERAZIONI.
% tol     : TOLLERANZA CRITERIO DI ARRESTO.
% g       : EQUAZIONE DA RISOLVERE x=g(x)
x_old=x0;
xhist=[x_old];
for index=1:maxit
    x_new=feval(g,x_old);
    if norm(x_new-x_old,inf) < tol
       return
   else
      xhist=[xhist; x_new];
      x_old=x_new;
   end
end
fprintf('\n \t [WARNING]: RAGGIUNTO MASSIMO NUMERO DI ITERAZIONI ');

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

>> g=inline('(1/2)*cos((1/2)*sin(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
>> feval(g,0.48640515466592)
ans =
   0.48640515466592
>>

Nelle ultime due righe dell’esperimento abbiamo verificato che effettivamente per 0 48640515466592, ( ). Vediamo se anche il sistema nonlineare

\begin{matrix}\frac{1}{2} cos( ) 0 \\ \frac{1}{2} sin( ) 0 \end{matrix}  (3)

è risolto correttamente.

• Se 0 48640515466592, allora dalla seconda equazione

\frac{1}{2} sin(0 48640515466592) 0 23372550195872

• Verifichiamo che anche la prima equazione sia risolta. In effetti si ha

0 48640515466592 e

\frac{1}{2} sin( ) 0 23372550195872

che coincide proprio con .

In Matlab/Octave

>> format long;
>> x=0.48640515466592;
>> y=0.5*sin(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. Riproponiamo quanto visto nell’esempio precedente

\begin{matrix}\frac{1}{2} cos( ) 0 \\ \frac{1}{2} sin( ) 0 \end{matrix}

in modo differente.

Siano v ( _{})_{}[ ; ]

1(v)  \frac{1}{2} cos( _{2}) 0  (4)
2(v)  \frac{1}{2} sin( _{1}) 0

per (v) [ _{1}(v); _{1}(v)] abbiamo v (v). Mostriamo che : �^{2} �^{2} è contrattiva. Dalla formula di Taylor in più variabili (cf. [4, p.192]),

(v) (v_{0}) ( ^{}( )) (v v_{0})

essendo il segmento che congiunge v con v_{0} e ^{}la matrice Jacobiana definita per componenti come

( ^{})_{}(v) \frac{_{}(v)}{_{}}

Nel nostro esempio quindi,

( ^{})_{}(v) \begin{matrix}0 & ( 1/2)sin( _{1}) \\ ( 1/2)cos( _{2}) & 0 \end{matrix}

Essendo la norma infinito di una matrice definita da

max_{}_{}_{}

si ha

( ^{})(v) _{}max_{}_{}( ^{})(v)_{}(5)
max( ( 1/2)sin( _{1}) ( 1/2)cos( _{2}) )\frac{1}{2}  (6)

Dalla formula di Taylor si ha così che

(v) (v_{0}) ( ^{}( )) (v v_{0})

e calcolando la norma infinito di ambo i membri 0

(v) (v_{0}) _{}( ^{}( )) _{}(v v_{0}) _{}

da cui : �^{2} �^{2} è una -contrazione con 1/2.

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

Scriviamo il file g1

function res=g1(v)
res(1,1)=0.5*cos(v(2));
res(2,1)=0.5*sin(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.48640524877124   0.23372554355416
0.48640514984910   0.23372554355416
0.48640514984910   0.23372549982964
0.48640515491247   0.23372549982964
0.48640515491247   0.23372550206770
0.48640515465330   0.23372550206770
0.48640515465330   0.23372550195314
0.48640515466657   0.23372550195314
0.48640515466657   0.23372550195901
0.48640515466589   0.23372550195901
0.48640515466589   0.23372550195871
0.48640515466592   0.23372550195871
0.48640515466592   0.23372550195872
0.48640515466592   0.23372550195872
>>

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

1. Gli errori assoluti compiuti, posto

[0 486405154665920 23372550195872]

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

ferma

( _{1} ) ( _{1} _{}) 0 1

abbiamo

>> format short e
>>
>> % CALCOLIAMO ERRORI ASSOLUTI.
>> xhist_err=abs(xhist-repmat([0.48640515466592   0.23372550195872],24,1));
>> xhist_abserr=sqrt((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=sqrt( (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.3412e-001  5.0000e-001
1.4855e-002  2.3971e-001
6.0283e-003  1.4297e-002
7.6760e-004  6.2976e-003
3.1244e-004  7.3797e-004
3.9270e-005  3.2625e-004
1.5982e-005  3.7756e-005
2.0101e-006  1.6689e-005
8.1807e-007  1.9326e-006
1.0289e-007  8.5424e-007
4.1873e-008  9.8922e-008
5.2664e-009  4.3725e-008
2.1433e-009  5.0634e-009
2.6956e-010  2.2381e-009
1.0971e-010  2.5917e-010
1.3796e-011  1.1456e-010
5.6147e-012  1.3266e-011
7.0775e-013  5.8636e-012
2.8805e-013  6.7901e-013
3.4630e-014  3.0012e-013
1.4144e-014  3.4750e-014
3.3766e-015  1.5377e-014
1.9767e-015  1.7764e-015
>>

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, ragion per cui il risultato non deve essere ivi considerato.

Osservazioni. Si noti che

• in generale un’equazione si scrive come ( ) 0, mentre nella formulazione di punto fisso si descrive come ( ); evidentemente non è un problema, in quanto se ( ) 0 allora ( ) per ( ) : ( ) ; • l’esempio mostrato nella sezione è particolarmente semplice, in quanto può essere ricondotto ad un problema unidimensionale; ciononostante abbiamo mostrato che può essere risolto da un metodo che lavora in più variabili.

2 Metodo di Newton

Dalla formula di Taylor (multivariata), se la funzione : �^{}�^{}è sufficientemente regolare, allora per v_{k 1} sufficientemente vicino a v_{k}

(v_{k 1}) (v_{}) ( ^{}(v_{})) (v_{k 1} v_{})  (7)

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

Da (7), supposto che la matrice ^{}(v_{0}) sia invertibile, e che (v_{k 1}) 0, cioè v_{k 1} sia una approssimazione di uno zero di , abbiamo

0 (v_{k 1}) (v_{}) ( ^{}(v_{})) (v_{k 1} v_{})  (8)

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

v_{k 1} : v_{k} ( ^{}(v_{}))^{1} (v_{})  (9)

Talvolta, la successione di Newton viene descritta come

\begin{matrix}^{}(v_{})) _{1}  (v_{}) \\ v_{k 1} : v_{k} _{1} \end{matrix}  (10)

sottolineando che l’inversa di ^{}(v_{})) non deve necessariamente essere calcolata, bensì bisogna risolvere un sistema lineare la cui soluzione _{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]. Vediamo di rifare l’esempio della sezione precedente, cioè

\begin{matrix}\frac{1}{2} cos( ) 0 \\ \frac{1}{2} sin( ) 0 \end{matrix}  (11)

mediante il codice implementante il metodo di Newton

function [v_hist, step_hist, residual_hist]=newton_sistemi(v0,maxit,tol)
% INPUT:
% v0      : APPROSSIMAZIONE INIZIALE.
% maxit   : NUMERO MASSIMO DI ITERAZIONI.
% tol     : TOLLERANZA DEL CRITERIO D'ARRESTO.

% F       : FUNZIONE.
% DF      : FUNZIONE DA CUI VIENE CALCOLATA LA JACOBIANA DI F.
% OUTPUT:
% v_hist : VETTORE CONTENENTE LE APPROSSIMAZIONI DELLA SOLUZIONE.
% step_hist    : VETTORE DEGLI STEP.
% residual_hist: VETTORE DEI RESIDUI.
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=[];
for 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];
    if max(loc_step,residual) < tol
        return;
    else
        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
fprintf('\n \t [WARNING]: NUMERO MASSIMO DI ITERAZIONI RAGGIUNTO.')

function res=F(v)

res(1,1)=0.5*cos(v(2))-v(1);
res(2,1)=0.5*sin(v(1))-v(2);

function res=DF(v)
% res(1,1)=0.5*cos(v(2));
% res(2,1)=0.5*sin(v(1));
res=[-1 -0.5*sin(v(2)); 0.5*cos(v(1)) -1];

Quindi scriviamo il programma principale driver_newton che setta i parametri

della procedura newton_sistemi
clear;
v0=[0; 0]; maxit=50; tol=10^(-10);
[v_hist, step_hist, residual_hist]=newton_sistemi(v0,maxit,tol);

Otteniamo così

>> driver_newton
>> v_hist
v_hist =
            0            0
  5.0000e-001  2.5000e-001
  4.8646e-001  2.3377e-001
  4.8641e-001  2.3373e-001
  4.8641e-001  2.3373e-001
  4.8641e-001  2.3373e-001
>> step_hist
step_hist =
  5.5902e-001
  2.1132e-002
  7.5293e-005
  7.7616e-010
            0
>> residual_hist
residual_hist =
  5.0000e-001
  1.8640e-002
  6.7480e-005
  6.7927e-010
            0
            0
>>

Dal vettore degli step _{1} _{}e da quello dei residui ( _{}) 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(size(v_hist,1),:);
>> err_vett=v_hist-repmat(alpha,size(v_hist,1),1);
>> abs_err_vett=sqrt( (err_vett(:,1)).^2 + (err_vett(:,2)).^2 )
abs_err_vett =
  5.3965e-001
  2.1206e-002
  7.5294e-005
  7.7616e-010
            0
            0
>>

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

err_vett=v_hist-repmat(alpha,size(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
>>

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

>> help 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.
>>

References

[1] K. Atkinson, An Introduction to Numerical Analysis, Wiley, (1989).

[2] K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer Verlag, (2001).

[3] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[4] G. Gilardi, Analisi due, seconda edizione, Mc Graw-Hill, 1996.

[5] J.H. Mathews e K.D. Fink, Numerical methods using Matlab, terza edizione, Prentice-Hall, 1999.

[6] A. Quarteroni, F. Saleri Introduzione al Calcolo Scientifico, Springer, (2002).

[7] The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler .

[8] G. Rodriguez, Algoritmi numerici, Pitagora editrice, 2008.