Appunti dal corso di calcolo numerico

Carli Samuele
Matricola: 4036768
E-mail: winsucks@tin.it

Gennaio 2007

Indice

1 Errori ed aritmetica finita
 1.1 Valutazione dell’errore
  1.1.1 L’errore assoluto
  1.1.2 L’errore relativo
 1.2 Sorgenti di errore
  1.2.1 Errori di discretizzazione
  1.2.2 Errori di convergenza
  1.2.3 Errori di round-off
 1.3 Rappresentazione dei numeri reali
  1.3.1 I numeri di macchina
  1.3.2 Overflow e underflow
 1.4 Lo standard IEEE 754
 1.5 Condizionamento di un problema
  1.5.1 Condizionamento delle operazioni elementari
2 Radici di una equazione
 2.1 Il metodo di bisezione
  2.1.1 Criteri di arresto:
  2.1.2 Ordine di convergenza
 2.2 Il metodo di Newton
  2.2.1 Convergenza locale
  2.2.2 Criteri di arresto
  2.2.3 Ripristino della convergenza quadratica in caso di radici multiple
 2.3 Metodi quasi-Newton
  2.3.1 Metodo delle corde
  2.3.2 Metodo delle secanti
3 Risoluzione di sistemi lineari
 3.1 Matrici quadrate
  3.1.1 Matrici diagonali
  3.1.2 Matrici triangolari
  3.1.3 Matrici ortogonali
  3.1.4 Propriet
 3.2 Metodi di fattorizzazione
 3.3 Fattorizzazione LU
  3.3.1 Costo computazionale
 3.4 Matrici a diagonale dominante
 3.5 Fattorizzazione LDLT
  3.5.1 Costo computazionale
 3.6 Pivoting
 3.7 Condizionamento del problema
 3.8 Sistemi lineari sovradeterminati
  3.8.1 Esistenza della fattorizzazione QR
4 Approssimazione di funzioni
 4.1 Interpolazione polinomiale
 4.2 Forma di Lagrange e forma di Newton
 4.3 Interpolazione di Hermite
 4.4 Errore nell’interpolazione
 4.5 Condizionamento del problema
 4.6 Ascisse di Chebyshev
 4.7 Funzioni spline
  4.7.1 Spline cubiche
  4.7.2 Calcolo di una spline cubica
 4.8 Approssimazione polinomiale ai minimi quadrati
5 Formule di quadratura
 5.1 Formule di Newton-Cotes
  5.1.1 Formula dei trapezi
  5.1.2 Formula di Simpson
  5.1.3 Condizionamento del problema
 5.2 Errore e formule composite
  5.2.1 Formule composite
 5.3 Formule adattative
  5.3.1 Formula dei trapezi
  5.3.2 Formula di Simpson
6 Altre implementazioni
 6.1 Modifica al metodo delle secanti
 6.2 Fattorizzazione LU
  6.2.1 Algoritmo ottimizzato
  6.2.2 Ottenere i fattori L ed U
  6.2.3 Risolvere il sistema
  6.2.4 Ottimizzazione per sistemi tridiagonali
 6.3 Fattorizzazione LDLT
  6.3.1 Algoritmo ottimizzato
  6.3.2 Trasformare una matrice sdp in un vettore utile all’utilizzo dell’algoritmo ottimizzato
  6.3.3 Piccolo miglioramento
  6.3.4 Ottenere i fattori L e D
  6.3.5 Risolvere il sistema
 6.4 Fattorizzazione LU con pivoting parziale
  6.4.1 Algoritmo ottimizzato
  6.4.2 Risolvere il sistema
 6.5 Fattorizzazione QR
  6.5.1 Ottenere i fattori Qed R
  6.5.2 Risolvere il sistema
 6.6 Algoritmo di Horner
  6.6.1 Algoritmo originale
  6.6.2 Algoritmo greneralizzato
  6.6.3 Polinomio interpolante di Newton
  6.6.4 Polinomio interpolante di Hermite
 6.7 Spline
  6.7.1 Spline lineare
  6.7.2 Spline cubica naturale
  6.7.3 Spline cubica not-a-knot
  6.7.4 Spline cubica completa
  6.7.5 Spline cubica periodica
  6.7.6 Approssimazione ai mini quadrati
 6.8 Formule di quadratura
  6.8.1 Formula dei trapezi adattiva non ricorsiva
  6.8.2 Formula di Simpson adattiva non ricorsiva

Capitolo 1
Errori ed aritmetica finita

1.1 Valutazione dell’errore

1.1.1 L’errore assoluto

Si supponga di avere a disposizione un metodo numerico che dia come risultato un numero ˜x approssimazione del risultato esatto x . possibile definire la grandezza errore assoluto come:

x = ˜x x

o in altra forma:

˜x = x + x.

Questa grandezza, bench fornisca un’indicazione precisa del valore dell’errore commesso dal metodo matematico preso in esame, non permette di valutare quanta influenza abbia l’errore sul risultato ottenuto. Ad esempio, un errore assoluto x = O(105) potrebbe avere un peso accettabile se il nostro risultato corretto x fosse dell’ordine di grandezza O(1) o superiore, ma non si avrebbe un risultato valutabile se x fosse anch’esso O(105).

1.1.2 L’errore relativo

Volendo conoscere quanto un errore influenzi il risultato quando x0, si definisce l’errore relativo:

εx △x-
 x = ˜x−-x-
 x

da cui:

˜x = x(1 + εx), e quindi ˜x-
x = 1 + εx

ovvero l’errore relativo deve essere comparato a 1: un errore relativo vicino a zero indicher che il risultato approssimato molto vicino al risultato esatto,mentre un errore relativo uguale a 1 indicher la totale perdita di informazione.

Esempi

Sia x = π 3.1415 = x˜. L’errore relativo quindi: εx = ˜x−x
 x 3.1415−3.141592653-
  3.141592653 = 0.000029492. Notiamo che la formula: log10(εx) restituisce all’incirca il numero di cifre significative corrette all’interno di ˜x: in questo caso il risultato del calcolo 4.53, che abbastanza vicino alla realt di 5 cifre significative corrette.

1.2 Sorgenti di errore

1.2.1 Errori di discretizzazione

I calcolatori elettronici non sono in grado di lavorare su un insieme numerico continuo, ma sono limitati a operare su un insieme numerico finito e discreto. In generale, per, i problemi matematici sono naturalmente pensati su un insieme continuo, ed necessario sostituirli con opportuni problemi discreti che il calcolatore in grado di gestire. Ad esempio, si supponga di voler calcolare un’approssimazione della derivata prima di una funzione derivabile f in un punto assegnato x0:

f(x0) f(x0 +-h)-− f-(x0)
       h con h > 0.

L’incremento h, per quanto piccolo, una quantit finita.
Consideriamo l’espansione in serie di Taylor con resto al secondo ordine:

f(x0 + h) = f(x0) + hf(x0) + h2
 2!f′′(ψ) con ψ (x0,x0 + h)

sostituendo si ha:

f(x0) f(x0 +-h)-− f-(x0)
       h = f(x0) + h2f′′(ψ)

dove il termine h
2f′′(ψ) quantifica l’errore di discretizzazione, che in questo caso un O(h).
Consideriamo invece la formula:

f(x0) = f(x0 +-h)-− f-(x0-−-h)
        2h con h > 0

e lo sviluppo in serie di Taylor di f(x) di punto inziale x0 con resto al terzo ordine si ottiene:

f(x0) = (f(x )+ hf′(x )+ h2f′′(x )+ h3f′′′(ψ))− (f(x )− hf′(x )+ h2f′′(x )− h3f′′′(η))
----0--------0---2!----0---3!------------0-------0----2!----0----3!------
                                   2h con ψ (x0,x0 + h) e η (x0 h,x0)

ovvero:

f(x0) = f(x0) + h2
12-(f′′′(ψ) f′′′(η))

In questo caso l’errore commesso O(h2), minore rispetto al caso precedente.
Questo esempio mostra come, approssimando un problema continuo a un problema discreto in due modi diversi, vi possa essere una grossa differenza nell’errore di discretizzazione commesso.

1.2.2 Errori di convergenza

Spesso per risolvere un problema possibile definire un metodo numerico di tipo iterativo che approssimi il risultato esatto, qui indicato con x, attraverso una successione di risultati intermedi {xn}, ricavabili mediante la funzione di iterazione del metodo:

xn+1 = Φ(xn), n = 0,1,2...,

Ovviamente, affinch il calcolo iterativo possa avere inizio, necessario definire un x0 approssimazione iniziale su cui eseguire il metodo la prima volta.
Condizione minimale che rende il metodo iterativo utile che sia verificata la condizione di consistenza: Φ(x) = x.

Un metodo iterativo si potr dire convergente se:

 lim xn = x
n→ ∞

tuttavia, il risultato corretto disponibile solo dopo infinite iterazioni. Ovviamente al di fuori del mondo matematico si sar costretti ad eseguire un numero finito di iterazioni, rendendo quindi necessario definire un opportuno criterio di arresto che all’n-esimo ciclo arresti la successione e determini l’errore:

x = xn x

detto errore assoluto di convergenza del metodo Φ(xn).

Teorema

Il metodo iterativo appena descritto, convergente a x, deve verificare la condizione di consistenza; ovvero la soluzione cercata deve essere un punto fisso per la funzione di iterazione che definisce il metodo.
Dimostrazione:
Supponiamo che non valga la condizione di consistenza, ovvero Φ(x)x. Allora se all’n-esimo passo dell’iterazione xn = x, xn+1x e non sarebbe vero che limn→∞xn = x, ovvero la funzione non sarebbe convergente e non sarebbe possibile arrivare a un risultato nemmeno dopo infinite iterazioni.

Esempio

Il metodo iterativo

      1 (      2)
xn+1 = 2 xn + xn- ,  n = 0,1,2..., x0 = 2

definisce un metodo numerico convergente per calcolare √ -
  2.
Dimostrazione:
√ -
  2 un punto fisso per la funzione, infatti:

               √ -
1 (√ -   2 )     2   1     2   √ -
2    2+ √2-  = -2-+ √2- = √2-=   2.

e quindi vale la condizione di consistenza.
La funzione converge a √ -
  2 e a ogni passo dell’iterazione la distanza dal valore corretto diminuisce, ovvero:

 (       )          (           )
1      -2-   √ -  1         --2--   √ -
2 xn + xn  −   2 < 2 xn+1 + xn+1  −   2.

Dal momento che consideriamo x0 = 2, la funzione si avvicina a √ -
  2 per valori positivi. Quindi:

                  (                       )
 (      )             (     )
1  x+ 2-  − √2 < 1|| 1  x + 2- + --(-2---)-|| − √2,
2     x          2( 2      x    1      2- )
                                2  x + x

          (             )             (   2  2    )
x2 +-2  1  x2-+-2  --2-     x2-+-2   1  (x+22x)-+-4x-
 2x   < 2    2x  +  x2+2- ,    2x  <  2   2x(x2+2)   ,
                     2x                       2x

 2      (x2+2)2+8x2
x-+-2-< ----2x----,  (x2 + 2)2 < (x2 + 2)2 + 8x2, 8x2 > 0, x > 0.
  2x     2(x2 + 2)

Ovvero per ogni x > 0 punto di innesco della sequenza, questa converge a √-
 2. Cvd.

Esempio

Il metodo iterativo:

       xnxn−1 + 2
xn+1 = --x-x-----, n = 1,2,..., x0 = 2,x1 = 1.5,
          n n−1

fornisce una successione di approssimazioni convergente a √2.
Ecco un esempio di successione, arrestata a un errore di convergenza 1012, considerando √2- = 1.414213562374:




n xn εconv



0 2 4.1 101



1 1.5 6.0 102



2 1.42 4.1 103



3 1.4137 3.6 104



4 1.41423 1.2 105



5 1.41421358 1.2 108



6 1.414213562373394.3 1013



1.2.3 Errori di round-off

Gli errori di round-off si verificano quando si cerca di rappresentare dei numeri che contengano intrinsecamente una quantit infinita di informazione (come ad esempio i numeri irrazionali) all’interno di una memoria dalle dimensioni finite che quindi obbliga a utilizzare un’approssimazione di questi.

Rappresentazione dei numeri interi

Un numero intero rappresentato, all’interno di un calcolatore, come una stringa del tipo:

α0α1....αn

per cui, assegnata una base b :

α0 ∈{+,−}, αi ∈{0,1,...,b 1}, i = 1,....,N.

A questa stringa viene fatto corrispondere il numero intero:

n = ({  ∑Ni=1 αibN −i, se α0 = +,

(  ∑N   αibN −i − bN, se α0 = − .
     i=1

In particolare nel caso di α0 positivo si ha che:

                                                       ∑N
0 ≤ α1bN−1 +...+ αN b0 ≤ (b− 1)bN−1 + ...+ (b− 1)b0 = (b− 1) bN−i
                                                       i=1

ovvero il massimo numero positivo rappresentabile :

      N              N −1
(b− 1)∑  bN−i = (b− 1) ∑ bi = (b− 1)bN-−-1-= bN − 1
      i=1              i=0           b − 1

e quindi l’intervallo di numeri positivi rappresentabile : [0,bN 1]
mentre quello di numeri negativi : [bN,1].
Mediante questo sistema possibile rappresentare senza errore tutti gli interi nell’insieme {            }
 − bN,...,bN − 1.

1.3 Rappresentazione dei numeri reali

Un numero reale rappresentato come una stringa del tipo:

α0α1...αmβ1β2...βs

dove α0 ∈{+,−}ij ∈{0,1,...,b 1} con i = 1...m e j = 1...s ed α10 per la normalizzazione.
La stringa suddetta rappresenta il numero reale:

     (∑m       )        ∑s
r =     αib1−i be−ν,e =   βjbs−j
      i=1               j=1

in cui lo shift ν fissato e chiameremo le quantit ρ = i=1nαib1i mantissa e μ = eν esponente del numero reale. La presenza dello shift serve a evitare di sprecare un bit dell’esponente per discriminarne il segno.
In particolare si ha che:

ν e bs −-1
b − 1(b 1) ν = bs 1 ν

ovvero e = (b 1) j=1sbsj ν.

Teorema

Secondo la rappresentazione appena definita, 1 ρ < b.
Dimostrazione:
Il minimo numero rappresentabile si ha nel caso α1 = 12 = 0...αm = 0: 1.◜◞◟◝
 000  m1.
Il massimo invece si avr nel caso α1 = ... = αm = (b 1): (b 1).(◟b−-1).◝..◜(b−-1)◞m1
In quest’ultimo caso si ha:

                      ∑m       1−i        m∑−1 −1 i        1− b−m
(b− 1).(◟b−-1).◝.◜.(b−-1)◞=    (b− 1)b   = (b− 1)    (b  ) = (b− 1)1-−-b−-1 =
            m−1       i=1                  i=0

             −m
= (b− 1)-1−-b----= b(1− b−m)
        (b− 1)b− 1

che sicuramente minore di b in quanto (1 bm) < 1. Cvd.

Il valore assoluto minimo e massimo tra i numeri di macchina rappresentabili sono rispettivamente:

r1 = b−ν
r2 = b(1− b−m)b((bs− 1)−ν) = (1− b−m )bs−ν .

1.3.1 I numeri di macchina

Definiamo adesso l’insieme dei numeri di macchina:

     (                                         )
     {              e    ∑m    1−i    ∑s    s−j}
M  = (x ∈ ℝ | x = ρb ,ρ =   αib   ,e =   βjb   ) ∪{0} .
                          i=1          j=1

In particolare si ha che:

M  ⊆ [− r2,r1]∪ {0}∪ [r1,r2] = I

ovvero i numeri di macchina appartengono a un sottoinsieme della retta reale che per definizione un intervallo denso, quindi infinito e non numerabile.

Teorema

M ha un numero finito di elementi.
Infatti, α2...αmβ1...βs possono assumere bs+m1 valori, α1 ne pu assumere (b 1) e α0 invece soltanto 2. Ovvero i numeri rappresentabili sono: 2(b 1)bs+m1 < 2bs+m.
In particolare |ρ| = α1α2...αm α1 1, ma anche:
|ρ| = α1α2...αm (b 1).◜-----◞◟----◝
(b− 1)...(b− 1)  m1 b. Cvd.

Funzione floating

In quanto M ha un numero finito di elementi e I ne ha un numero infinito, necessario definire una funzione fl(x), detta funzione floating, che associ ad ogni numero reale x I un corrispondente numero di macchina:

fl : x ∈ I → fl(x) ∈ M

Troncamento

Sia x = α1α2...αmαm+1...be I
Possiamo definire:

˜x = fl(x) =  α1α2...αmbe.

Arrotondamento

Sia x = α1α2...αmαm+1...be I
Possiamo definire

                            (
                            |{               b
˜x = fl(x) =  α α ...˜α be,α˜ =   αm, se αm+1 < 2    .
             1 2   m    m   |(  αm + 1, se αm+1 ≥ b
                                               2

Teorema

x I,x0:

                            {  1−m
                               b   T roncamento
˜x = fl(x) = x(1 + εx),|εx| ≤ u =   1b1−m Arrotondamento
                               2

dove con u indichiamo la precisione di macchina.
Dimostrazione (nel caso di troncamento):

      |x − ˜x|  (α1α2...αm αm+1...)be−ν − (α1α2...αm )be−ν
|εx| =--|x|--= --------(α1α2...αmαm+1...)be−ν-------- =

      m◜−◞◟1◝
    0.000 αm+1...    (αm+1 αm+2...)b− m   bb−m    1−m
=  α1α2...αm-αm+1...≤ --------1------- ≤ --1--= b    = u

e quindi:

x ∈ I ⇒ fl(x ) = x(1− εx),|εx| ≤ u

Osservazione: log10|εx|≃ numero di cifre corrette nella mantissa.
In particolare log10|εx|≤−log10|u| e m = 1 log10|u|
La quantit u prende il nome di precisione di macchina

Dimostrazione:
nel caso di arrotondamento:
Adesso dobbiamo distinguere tra i due casi, αm+1 < b
2 oppure αm+1 b
2.
Nel primo:

|ε | = |x-−-˜x|= (α1α2...αm-αm+1...)be−ν-− (α1α2...α˜m-)be−ν =
  x     |x|            (α1α2...αmαm+1...)be−ν

     ◜m◞−◟1◝
    0.000 αm+1...   (αm+1αm+2...)b−m    12bb− m   1 1−m
= α1α2...αm-αm+1... ≤ -------1--------≤  --1---= 2b

osservazione: (αm+1αm+2...) minore di b∕2 per ipotesi.

Nel secondo:

                                        e−ν             e−ν
|εx| = |x-−x˜|= |x−-˜x|=  (α1α2...αm-αm+1...)b---−-(α1-α2...α˜m-)b---=
       |x|      x              (α1α2...αm αm+1...)be−ν

       m−1
       ◜◞◟◝                       −m    1  −m
= |−-0.000 αm+1...|≤ (αm+1αm+2...)b-- ≤ 2bb---= 1b1−m
  α1α2...αm αm+1...          1             1     2

osservazione: (αm+1αm+2...) sicuramente minore di b∕2.

1.3.2 Overflow e underflow

Pu capitare di dover rappresentare un numero non contenuto in I, e in questo caso ci troveremmo in una condizione di errore. I casi possibili sono:

Possiamo adesso completare la definizione di fl(x):

       (|  0, x = 0
       |{  x(1+ ε ), |ε | ≤ u, r ≤ |x| ≤ r
fl(x) = | underfxlow, x0 ≤ |x| <1r      2
       |(  overf low, |x| ≥ r    1
                        2

1.4 Lo standard IEEE 754

Lo standard IEEE 754 stato creato per definire la rappresentazione dei numeri reali su calcolatore, allo scopo di uniformare i risultati dei calcoli eseguiti su piattaforme di calcolo differenti; cosa molto importante in vari ambiti scientifici infatti la riproducibilit del risultato. Questo standard stato definito e ottimizzato sulla base binaria.
Si sono scelte due rappresentazioni per i numeri reali, una a precisione singola (32 bit), e una a doppia precisione (64 bit).
Nelle due rappresentazioni i bit sono cos distribuiti:

In quanto i numeri rappresentati sono da intendersi come normalizzati, si stabilito che α1 non venga rappresentato, in quanto sempre uguale a 1: si ottiene un bit in pi dedicabile alla mantissa, che quindi in realt di 24 bit per la rappresentazione a singola precisione e di 54 per quella a doppia precisione. Indicheremo con f la parte frazionaria della mantissa, ovvero l’unica parte rappresentata in macchina: α2...αm.
Per poter gestire anche i casi di eccezione, stato deciso di dedicare il primo e l’ultimo numero rappresentabile nell’esponente a queste situazioni, ovvero:

per quanto riguarda la singola precisione e:

per quanto riguarda la doppia.
In questo modo si ottiene un range di numeri rappresentabili di circa 1038 1038 per i numeri in singola precisione e 10308 10308 per quelli in doppia.
La rappresentazione scelta dallo standard per arrotondamento, ma con una piccola differenza da quello esposto precedentemente:
fl(x) definito come il numero di macchina pi vicino a x, ma in caso di due numeri di macchina equidistanti, viene selezionato quello il cui ultimo bit della mantissa 0:

(
{ α˜m  = 0 se αm+1 = 1
(
  α˜m  = αm  se αm+1 = 0

Questo tipo di arrotondamento chiamato arrotondamento al pari (round to even).

Esempi

Il pi grande numero di macchina si avr nel caso:

α0 = +; α1 = ... = αN = (b 1); β1 = ... = βs = (b 1);

In particolare, secondo lo standard IEEE754:
m=53, s=11, ν=1023 se il numero normalizzato. Quindi:

      0  − 1   − 53  (((bs−1)−1)−1023)   (2046−1023)∑52  −i       −53 1023
n = ((b ).(b )...(b ))b             = b            b  = b(1− b  )b   =
                                              i=0

        −53 1023
= 2(1− 2   )2   = 1.7976931e +308

mentre il risultato della funzione realmax di Matlab restituisce: 1.7977e+308.
Il pi piccolo numero di macchina normalizzato positivo si ha invece nel caso:

α0 = +; α1 = 12 = ... = αN = 0; β1 = ... = βs1 = 0; βs = 1.

Quindi secondo lo standard:

n = (1.0◟.◝..◜0◞)b1−1023) = b−1022) = 2.2250739e− 308
       52

mentre il risultato della funzione realmin di Matlab restituisce: 2.2251e-308.
Il pi piccolo numero di macchina denormalizzato si ha quando:

n = (0.0◟.◝..◜0◞1)b0−1022 = b−51b−1022 = b−1073 = 9.8813129e− 324.
       51

La precisione di macchina :

u = 1b1−m = 1.110223e− 16
    2

1.5 Condizionamento di un problema

Operazioni elementari in aritmetica finita

importante soffermarsi sul comportamento delle operazioni elementari in aritmetica finita. Lavorando su numeri interi, se il risultato dell’operazione cade all’interno dell’insieme di rappresentabilit, le operazioni coincidono con quelle algebriche. Considerando invece i numeri interi, le operazioni saranno definite solo su numeri di macchina e dovranno avere per risultato ancora un numero di macchina; ovvero ad esempio nel caso dell’addizione si ha che x y = fl(fl(x) + fl(y)). In particolare, le operazioni in macchina sui numeri reali generalmente non godono di tutte le usuali propriet algebriche, come ad esempio la propriet associativa e quella distributiva.

Conversione tra tipi

La conversione di tipo da intero a reale non presenta particolari problemi, in quanto, a parit di bit di rappresentazione, un intero sempre rappresentabile in forma reale, introducendo al massimo un errore dell’ordine della precisione di macchina u: questo infatti l’errore commesso approssimando tramite la funzione fl(x). Convertire un reale in un intero invece pu causare problemi, dovuti alla grande differenza nel range di rappresentabilit: circa 1038 per i reali e circa 1010 per gli interi.

Condizionamento di un problema

Supponiamo di dover risolvere un problema rappresentabile nella forma x = f(x), dove x sono i dati in ingresso, y rappresenta il risultato f un metodo numerico. Quello che ci troveremo in realt a risolverein macchina un problema del tipo = ˜f (˜x); ovvero eseguiremo una funzione perturbata (con operazioni eseguite in aritmetica finita, con possibile errori di discretizzazione e/o convergenza) su un dato perturbato (errore di rappresentazione sempre presente, possibile origine sperimentale del dato). Per semplicit ci limiteremo a studiare il problema = f(˜x), ovvero considereremo il metodo numerico come eseguito in aritmetica esatta.
Consideriamo quindi le grandezze ˜x = x(1 + εx) e = y(1 + εy) e studiamo la relazione che intercorre tra l’errore relativo in ingresso e quello in uscita, considerando un’analisi al primo ordine:

y+ yεy = f(x)+ f ′(x)xεx + O(ε2)
                           x

ovvero:

     ||      ||
|εy|  ||f′(x)x-|||εx| ≡ k|εx|.
          y

Definizione

Il valore di amplificazione k, che indica quanto gli errori iniziali possano essere amplificati sul risultato finale, denominato numero di condizionamento del problema.

Significato dei valori di k:

1.5.1 Condizionamento delle operazioni elementari

Moltiplicazione: y = x1x2
y(1+ εy) = x1x2(1 +ε1)(1+ ε2) =⇒ εy = ε1+ ε2 + ε1ε2 =⇒ |εy| ≤ |ε1|+|ε2| ≤ 2εx

Divisione: y = x1
x2
y(1 + εy) = x1(1-+ε1) ⇒ εy = ε1 −-ε2
          x2(1 +ε2)        1+ ε2

Somma algebrica: y = x1 + x2
y(1+εy) = x1(1+ ε1)+ x2(1+ ε2) =⇒ yεy = x1ε1 +x2ε2 = ⇒ εy = |x1ε1 +-x2ε2| ≤
                                                            y

   |x1ε1|+ |x2ε2|  |x1|+|x2|    |x1|+ |x2|
≤  -----|y|---- ≤ ---|y|---εx =-|x1-+-x2|εx
                               ◟---◝◜---◞
                                   K

Mentre la moltiplicazione ben condizionata vediamo che nel caso della somma dobbiamo differenziare due casi:

Radice: y = √ --
  x
                       1 1
f′(x) = 1x− 12 =⇒  K  = |2x2| = 1
       2               x12    2

Esempi

Matlab - Errori di rappresentazione

Il codice Matlab:

x = 0; delta = 0.1;  
while x ~=1, x=x+delta, end

esegue un loop infinito. Questo perch il numero 0.1 non rappresentabile in base binaria come numero finito ma solo come numero periodico, e quindi subisce un errore di approssimazione. In questo modo sommando ripetutamente 0.1 in binario non si arriva mai esattamente ad 1 e quindi il ciclo non si ferma. Questo codice potrebbe essere corretto sostituendo il controllo x1 con x <= 1, ma questo non garantisce che venga eseguito esattamente lo stesso numero di cicli che sarebbe stato eseguito dal codice precedente se avesse funzionato.

Efficienza: Individuare un algoritmo efficiente per calcolare ∘ -------
  x2 + y2:

Calcolare ∘ -------
  x2 + y2 eseguendo √-------
 xx + yy non il migliore dei modi: infatti, in caso x = y = r1,fl(x) = fl(y) = 0,fl(√ -
  0) = 0, oppure se x vicino ad r2, x2 non pu essere rappresentato. Questo problema pu essere risolto in modo pi efficiente eseguendolo invece in questo modo:

                            ∘-----------
m = max (|x|,|y|),∘x2--+-y2 = m |-x|2 + | y-|2.
                              m      m

che permette di effettuare il calcolo anche nelle condizioni estreme sopra descritte.

La somma algebrica non gode della propriet associativa e distributiva:

Eseguendo in matlab il codice:

((eps/2+1)-1)*(2/eps)  
(eps/2+(1-1))*(2/eps)

La variabile eps rappresenta la precisione di macchina. Pi precisamente il pi largo spazio relativo tra due numeri adiacenti del sistema floating point della macchina. Questo numero ovviamente dipendente dalla macchina, ma in macchine che supportano l’aritmetica floating point IEEE 64 bit, vale approssimativamente 2.2204e-16. Le due istruzioni sopra descritte danno come risultato:

Ovvero non vale la propriet associativa. Nel caso invece:

(1e300-1e300)*1e300  
(1e300*1e300)-(1e300*1e300)

si hanno i seguenti risultati:

Ovvero nel primo caso l’operazione viene eseguita correttamente, mentre nel secondo si incorre in overflow e non possibile svolgere il calcolo e quindi non vale la propriet distributiva.

Cancellazione Numerica

Supponiamo di voler calcolare: y = 0.12345678 0.12341234 0.00004444, utilizzando una rappresentazione decimale con arrotondamento alla quarta cifra significativa.
Otteniamo = 1.235 101 1.234 101 = 1 104.
L’errore relativo su y :

ε = 4.444∗10−-5 −-1∗10−-4≃ 5.560∗-10−-5≃ 1.25
 y       4.444 ∗10−5       4.444∗ 10− 5

mentre essendo εx1 = 3.5 104 e εx2 = 1 104, si ha che

   |3.5 ∗10−4|+ |− 1∗ 10− 4|
k =---------−4-------−4-- ≃ 5.5 ∗103
     |3.5∗ 10  − 1 ∗10  |

Il problema risulta essere mal condizionato.

Capitolo 2
Radici di una equazione

In questo capitolo tratteremo la risoluzione dell’equazione:
f(x) = 0, x ∈ r

in cui f : una funzione assegnata; ovvero siamo interressati a trovare, se esiste, uno zero reale della funzione. In generale ci possiamo trovare in tre diverse situazioni:

In particolare ci limiteremo a lavorare su funzioni tali che, definito [a,b] intervallo in con a < b, queste godano di certe propriet:

Questo garantisce l’esistenza di una radice di f all’interno dell’intervallo [a,b], ovvero la funzione, passando da positiva a negativa all’interno di questo intervallo, deve necessariamente attraversare l’asse delle ascisse in almeno un punto.

2.1 Il metodo di bisezione

Come primo passo dobbiamo scegliere un punto interno all’intervallo che sia una miglior approssimazione della radice (x), in particolare scegliamo il punto medio: x1 = a+b
 2.
A questo punto possono verificarsi solo tre casi:

2.1.1 Criteri di arresto:

Con questo metodo si ha che |x x1|≤b−a
 2. Ragionando per induzione, indicando con [ai,bi] l’intervallo di confidenza e con xi = (ai + bi)2 l’approssimazione della radice all’intervallo i-esimo, si ha che:

|x − xi| ≤ b-−-a
          2i

Ovvero ad ogni iterazione si dimezza la distanza massima della serie dalla soluzione. In questo modo, si pu pensare di voler ottenere il calcolo di una approssimazione della radice che disti al pi tolx dalla radice,e si ha che questa si raggiunge in un numero di iterazioni:

imax ≡ [log2(b− a)− log2(tolx)]

Adesso sappiamo che sicuramente, dopo imax itarazioni, abbiamo raggiunto la precisione richiesta; ma possiamo ancora stabilire un criterio efficace per accorgersi se si arrivati sufficientemene vicini alla soluzione prima di imax iterazioni. In linea di massima, se f(x) = 0, per la continuit di f si deve avere che il valore che la funzione assume in prossimit di x debba assumere valori vicini a 0. In questo modo possiamo pensare ad un controllo di arresto del tipo |f(x)| < tolf dove tolf deve essere scelta in modo da avere |x x| < tolx. Supponendo f C(1), sviluppando in serie di Taylor:

f(x) ≈ f(x)+ f′(x)(x − x) = f′(x)(x− x) ovvero |x − x| ≈-f(x)-
                                                     |f ′(x)|

e quindi tolf ≈|f(x)|tolx.

Osservazione:

Tanto pi grande sar |f(x)|, tanto meno stringente dovr essere il criterio di arresto sulla funzione. Infatti si ha che:

         f(xi)− f (x)                1
|xi − x| ≈--|f′(x)|--- ovvero  k = |f′(x)|.

dove k rappresenta il numero di condizionamento del problema. Questo metodo sar ben condizionato in presenza di funzioni che hanno una derivata prima molto grande nei pressi della radice, ma mal condizionato quando questa si avvicina o tende a zero.

Implementazione

Considerando che una buona approssimazione di f(x) ricavabile, a costo molto basso, in questo modo:

       f(b )− f(a)
f′(x) ≈---i------i
          bi − ai

si pu pensare al seguente algoritmo per l’applicazione del metodo di bisezione:

Costo computazionale

Per quanto riguarda il costo di questo algoritmo in termini di operazioni floating-point, si vede facilmente dal codice che questo algoritmo ad ogni ciclo esegue esattamente 3 somme, 4 moltiplicazioni ed una valutazione di funzione, tutto questo al massimo imax volte. Tenendo conto che imax definito come log2(b a) log2(tolx), si pu affermare che questo algoritmo ha un costo di andamento grossomodo logaritmico all’allargarsi dell’intervallo di confidenza. Il costo in termini di memoria anch’esso molto contenuto, si utilizzano infatto soltanto poche variabili.


PIC

Figura 2.1: Esempio di funzionamento del metodo di bisezione che mostra in che punti viene valutata la funzione (x 2)3 + 2 con l’intervallo di partenza [0,4]; gli estremi 0 e 4 sono stati esclusi dal grafico per migliorare la visualizzazione.

PIC

Figura 2.2: Esempio di funzionamento del metodo di bisezione che mostra come, nel caso la funzione abbia pi di uno zero all’interno dell’intervallo di confidenza, questo converga verso uno di essi. In questo caso la funzione y = sin(1
x) con l’intervallo di partenza [103,101].


2.1.2 Ordine di convergenza

Definito come ei = xi x l’errore commesso al passo i-esimo nell’approssimazione fornita da un dato metodo numerico, questo medoto convergente se:

lim  ei = 0
i→ ∞

Il metodo ha ordine di convergenza p se p il pi grande numero positivo per cui:

    |ei+1|
li→im∞ -|e-|p-= c < ∞
      i

dove c e detta costante asintotica dell’errore. In particolare si avranno metodi che convergono linearmente quando p = 1, esponenzialmente (quadraticamente, cubicamente...) per p > 1, e metodi non convergenti se p < 1.

                             2
|ei+1| ≈ c|ei|, |ei+2| ≈ c|ei+1| ≈ c |ei|

ovvero, per induzione:

|e   | ≈ ck|e|
  i+k       i

Abbiamo quindi che per i molto grandi, ovvero dopo molte iterazioni, anche k molto grande, e si deve avere che 0 c 1 altrimenti la funzione sarebbe divergente.

2.2 Il metodo di Newton

Il metodo di Newton si basa sull’idea di sfruttare, come successiva approssimazione della radice, il punto di intersezione della retta tangente alla funzione nel punto (x0,f(x0)) con l’asse delle ascisse, considerando x0 l’approssimazione iniziale. Consideriamo quindi y = f(xi) + f(xi)(x xi), ponendo y a zero si ha che:

         f(x )
x1 = x0 −-′-0--
         f (x0)

che definita per f(x0)0. Possiamo estendere il ragionamento a tutti i passi successivi:

           f(xi)
xi+1 = xi − f′(x-), i = 0,1,2...
              i

Teorema

Se xi la successione sopre descritta, che converge a x, zero semplice di f(x) funzione di classe C2, allora il metodo di Newton converge almeno quadraticamente.
Dim.
Per ipotesi abbiamo che xi x per i →∞. Quindi:

                 ′            1  ′′          2
0 = f (x) = f(xi)+ f (xi)(x− xi)+ 2 f (ψi)(x− xi) =

       (               )
= f′(xi)  f(′xi)-+ (x− xi) + 1 f′′(ψi)(x− xi)2 = f′(xi)(− xi+1 + x)+ 1f ′′(ψi)e2i
         f(xi)            2                      ◟-−◝e◜i+1-◞   2

e     1f ′′(ψ )
i+21-= ---′-i-
ei    2 f(xi)

che ben definita in quanto f(xi) sicuramente diverso da zero perch stiamo trattando il caso di una radice semplice e siamo in un opportuno intorno di x, e se il metodo converge si ha:

    1 f′′(ψi)-  1f-′′(x)
ili→m∞ 2 f′(xi) = 2 f′(x)

Osservazione:

Se xi x con i →∞ radice di molteplicit m > 1, allora l’ordine di convergenza del metodo di Newton solo lineare (p = 1) con c = (m 1)∕m

2.2.1 Convergenza locale

Rimane ancora da stabilire quali siano le condizione sufficienti a garantire la convergenza del metodo verso la soluzione. Per prima cosa necessario formalizzare un metodo iterativo per approssimare una radice x dell’equazione f(x) = 0:

xi+1 = Φ (xi), i = 0,1,2,...,

ed ovviamente la radice deve essere un punto fisso per questa funzione di iterazione. Possiamo quindi definire la funzione di iterazione per il metodo di Newton:

Φ(x) = x − f(′x).
          f(x)

Teorema del punto fisso

Sia Φ(x) la funzione di iterazione che definisce un metodo numerico. Supponiamo che δ > 0,0 L < 1 tali che x,y (x δ,x + δ) I; |Φ(x) Φ(y)|≤ L|x y|
allora:

1.
x l’unico punto fisso di Φ in I
2.
se x0 Ixi Ii 0
3.
xi x per i →∞

Dim(1).
Supponiamo per assurdo che ˜x | ˜x Φ(˜x) I ulteriore punto fisso per Φ, ovviamente xx˜. Ma allora si ha che: 0 < |x ˜x| = |Φ(x) Φ(˜x)|≤ L|x ˜x| ma essendo L < 1 si ha che |x ˜x| < |x ˜x|, assurdo.
Dim(2).
Sia x0 I⇐⇒|x x0| < δ.
Allora: |x xi| < δ e |x xi| = |Φ(x) Φ(x0)|≤ L|x x0| < Lδ < δ
e xi I.
Dim(3).
|xi x| = |Φ(xi1) Φ(x)|≤ L|xi1 x|, ovvero:
|xi x| = L|Φ(xi2) Φ(x)|≤ L2|xi2 x|
Ragionando per induzione possiamo quindi concludere che:
|xi x|≤ Li|x0 x|, che tende a zero per i che tende a infinito.

Se si assume f(x) C(2), tramite questo teorema possibile dimostrare la convergenza locale del metodo di Newton. In questo caso infatti la funzione di iterazione di classe C(1) e Φ(x) = 0. Fissato quindi un arbitrario L [0,1), esister un δ > 0 tale che |Φ(x)|≤ L,x I(xδ,x + δ). Sviluppando in serie di Taylor al primo ordine si ottiene che:

∀x,y ∈ I : |Φ (x) − Φy| = |Φ(x)− Φ(x)− Φ′(ε)(y − x )| ≤ L |x − y|,

con ε compreso tra x e y ovvero apppartenente ad I. La funzione di iterazione soddisfa le ipotesi del teorema, e quindi il metodo di Newton localmente convergente se x0 I

2.2.2 Criteri di arresto

Differentemente dal metodo di bisezione, a causa della convergenza solo locale del metodo di Newton, non possibile stabilire a priori il numero massimo di iterazioni necessarie per raggiungere la soluzione con l’accuratezza richiesta. Come nel caso precedente vorremmo stabilire un criterio sulla dimensione dell’errore, ovvero vogliamo fermarci quando |ei|≤ tolx dove tolx rappresenta appunto la tolleranza che siamo disposti ad accettare intorno alla soluzione. Notiamo che nel caso di convergenza verso radici semplici il metodo di Newton converge quadraticamente ed in prossimit della radice si ha:

|xi+1 − xi| = |xi+1 −-x+ x◟−◝x◜i ◞| = |ei − ei+1| ≈ |ei| ≤ tolx
            ◟−◝e◜i+1◞    ei

e quindi un controllo del tipo |xi+1 xi| < tolx appropriato ed equivalente al controllo |f(xi)|≤|f(xi)|tolx in quanto xi+1 = xi -f(′xi)
f (xi). Bisogna notare per che questo controllo essenzialmente non altro che una verifica sull’errore assoluto su x, ovvero tolx O(x). Bench questo possa risultare appropriato per valori piccoli di x, non lo necessariamente anche per valori grandi, dove risulta molto pi efficiente un controllo sull’errore relativo. Possiamo quindi definire la grandezza “errore relativo su x” come |ei|≤ rtol|x|≈ rtol|xi+1|, ed un criterio di arresto del tipo:

-------|xi+1-−-xi|--------
   tolx    +  rtolx|xi+1| ≤ 1
toller◟an◝z◜a◞assoluta   ◟--◝◜---◞
             tolleranzarelativa

che un criterio di arrsto ibrido che si arresta quando soddisfatto il criterio pi stringente. In particolare si ha che per valori piccoli di x, rtolx|xi+1| un valore piccolo rispetto ad x e quindi il criterio assoluto predomina, viceversa invece per valori di x grandi; ponendo rtolx uguale a zero questo criterio ritorna equilvalente a |xi+1 xi| < tolx.

Convergenza lineare

Quando la convergenza del metodo solo lineare, il criterio di arresto pu essere modificato. Si ha infatti che:

                                      m − 1
|xi+1 − xi| = |ei − ei+1| ≈ |ei|(1− c) con c =-m

Noi vogliamo |ei| < tolx e quindi:

     |x   − x |
|ei| =--i+1---i ≤ tolx
       (1− c)

Ovvero abbiamo ottenuto un criterio tanto pi efficiente quanto pi la molteplicit della radice alta (c 1). Noi conosciamo xi+1, e vogliamo che |ei+1|≤ tolx ma sappiamo che |ei+1| = c|ei| e possiamo quindi controllare -c-
1−c|xi+1 xi|≤ tolx. Il problema che non conosciamo c. Possiamo per calcolare le prime tre iterazioni (x0,x1,x2) in questo modo: |x1 x0|≈ (1 c)|e0|; |x2 x1|≈ (1 c)|e1|≈ (1 c)c|e0|; e quindi c |x  − x |
-2----1-
|x1 − x0|.

Implementazione

Visto quanto detto fino ad ora, possiamo pensare alla seguente implementazione del metodo di Newton:


PIC

Figura 2.3: Esempio di funzionamento del metodo di Newton che ne mostra il comportamento nella determinazione dello zero della funzione (x2)3+2 con punto di innesco x0 = 3.5.

PIC

Figura 2.4: Ecco come cambia il comportamento del metodo di Newton nella determinazione dello zero della funzione (x2)3 + 2 se il punto di innesco invece x0 = 2.5.



PIC

Figura 2.5: In questo caso il metodo di Newton non converge, stiamo cercando la radice di x3 + 5x a partire dal punto x0 = 1 ma il metodo genera una successione infinita del tipo {1,1,1,1,..},


Costo Computazionale

Rispetto al metodo di bisezione il metodo di Newton ha un costo leggermente pi alto: questo richiede infatti 7 operazioni floating-point e due valutazioni di funzione ad ogni iterazione. In questo caso il costo massimo della determinazione dello zero non determinabile a priori, e dipende dalla scelta di itmax effettuata dall’utente. Il costo in termini di memoria, come nel caso del metodo di bisezione anch’esso molto contenuto e si utilizzano soltanto poche variabili.

2.2.3 Ripristino della convergenza quadratica in caso di radici multiple

Molteplicit di una radice:

Una radice x ha moltiplicit esatta m se:

f(x) = f′(x) = ...= f m−1(x) = 0, ,f(m )(x) ⁄= 0.

Osservazione:

Se f(x) sviluppabile in serie di Taylor in x, sua radice di molteplicit esatta m, allora scrivibile nella forma f(x) = (xx)mg(x) con g(x) ancora sviluppabile in serie di Taylor in x e tale che g(x)0.

Molteplicit nota

Supponiamo di avere una radice x di molteplicit nota m. Allora f(x) sar del tipo (xx)m. Avendo x0 approssimazione di x:

x →  x = x −  --(x0 −-x)m--= x  − 1-(x − x) = x − m-f(x0)-
 0    1   0   m(x0 − x)(m−1)   0  m   0       0    f ′(x0)

ovvero

xi+1 = xi − m f(xi), i = 0,1,2...
            f′(xi)

In questo caso particolare si riesce ad ottenere la convergenza alla radice in un solo passo iterativo, utilizzando il fattore di correzione m.
In generale, f(x) una funzione del tipo (x x)mg(x) (con g(x) qualsiasi) e non si riesce quindi ad ottenere il risultato in un solo passo iterativo come nel caso particolare detto sopra; dimostrabile che questa modifica riesce per a ripristinare la convergenza quadratica del metodo di Newton verso radici di molteplicit nota.

Implementazione

Molteplicit non nota

Supponiamo invece di avere una radice x di molteplicit ignota. In questo caso al generico passo i-esimo si ha che ei cei1 ed ei+1 cei, con c costante asintotica dell’errore, incognita. Dalla combinazione delle due si ottiene che ei+1ei1 ei2 e quindi (xi+1 x)(xi1 x) (xi x)2. Se si suppone esatta l’uguaglianza, si ha che:

    -xi+1xi−1 −-x2i-   (i)
x ≈ xi+1 − 2xi + xi−1 ≡ x

In questo modo si pu pensare ad una macchinetta iterativa che agisce su due livelli:

1.
Livello interno in cui vengono eseguiti due passi del metodo di Newton, ovvero a partire da una approssimazione iniziale x0(0) si calcolano x1(0),x2(0) approssimazioni di Newton alla prima iterazione,
2.
Livello esterno in cui, con la formula:
       x(0)x(0)− (x(0))2
x(01)= -(20)-0--(0)-1-(0)
      x2  − 2x1 + x0

si ottiene una successiva approssimazione pi accurata x0(1) della radice da utilizzare per il successivo ciclo interno. Questo sistema, chiamato “metodo di accelerazione di Aitken”, permette di ripristinare la convergenza quadratica del metodo di newton nel caso di radici multiple di molteplicit ignota, al prezzo di aumentare il numero di operazioni da eseguire ad ogni iterazione (ciclo a due livelli).

Implementazione
Costo computazionale

Come descritto nella teoria, questo algoritmo effettua, ad ogni iterazione, due passi del metodo di Newton ed un passo di accelerazione, per un costo totale quindi di 16 operazioni floating point e 4 valutazioni di funzione. Per quanto riguarda l’occupazione di memoria, anche qui usiamo solo poche variabili.

2.3 Metodi quasi-Newton

I metodi detti quasi-Newtoniani sono delle varianti del metodo di Newton che permettono di risparmiare il calcodo della derivata prima della funzione ad ogni passo iterativo. Questi metodi sono descrivibili tramite uno schema generale che li accomuna:

x    = x − f(xi),  i = 1,2,..., φ ≈ f′(x)
 i+1    i   φi                 i

2.3.1 Metodo delle corde

Questo metodo si basa sull’osservazione che, in prossimit della radice, la derivata prima della funzione in esame varia di poco, e si pu quindi pensare di calcolarla una sola volta e di utilizzare poi un fascio di rette parallele a questa per approssimarla nelle iterazioni successive. Supponendo di essere in presenza su una funzione sufficientemente regolare ed abbastanza vicini alla radice, si pu definire il metodo come:

x   = x − f-(xi),  i = 0,1,2,...,
 i+1    i  f′(x0)

Il vantaggio di questo metodo risiede nel basso costo computazionale per iterazione (si valuta solo una volta f(x) e una sola volta per iterazione f(x)), ma si pu dimostrare che l’ordine di convergenza solo lineare.

Implementazione

Costo computazionale

Questo metodo ad ogni ciclo iterativo svolge soltanto 7 operazioni floating point ed una valutazione di funzione, il costo per iterata quindi molto basso. Bisogna inoltre notare come la convergenza del metodo sia molto legata al tipo di funzione e alla scelta del punto di partenza; se questo non abbastanza vicino alla radice oppure la funzione non abbastanza regolare, il mtodo pu convergere molto lentamente.


PIC

Figura 2.6: Esempio di funzionamento del metodo delle corde che ne mostra il comportamento nella determinazione dello zero della funzione (x2)3+2 con punto di innesco x0 = 3.5.

PIC

Figura 2.7: Le prime venti iterazioni del metodo delle corde per trovare una radice di x3 5x con punto di innesco x0 = 0.99: in questo caso il metodo converge molto lentamente e pu girare per pi di 100000 iterazioni senza avvicinarsi ad una delle radici con una tolleranza di 106


2.3.2 Metodo delle secanti

Tenendo conto della definizione della derivata prima mediante rapporto incrementale, si pu pensare di sfruttare due approssimazioni successive (ovvero due precedenti valutazioni di f(x)) per approssimarla. Considerando quindi di essere all’iterazione i-esima e di aver gi calcolato xi1,xi,f(xi1),f(xi), si approssima f(x) in questo modo:

             f(x + δ)− f(x )  f(x )− f(x  )
φ(x) ≈ f′(x) =------------ = ---i------i−-1
                   δ            xi − xi−1

Il metodo iterativo risultante:

xi+1 = f(xi)xi−1-−-f(xi−1)xi,  i = 1,2,...,
         f(xi)− f(xi−1)

con x0,x1 approssimazioni iniziali assegnate.
In questo modo si ha la possibilit di calcolare una sola volta per itearazione il valore di f(x), ottenendo cos un metodo con costo per iterazione molto basso; in pi, quando la funzione st convergendo si ha che xn molto vicino a xn1, e il rapporto incrementale risulta una buona approssimazione per la derivata. Questo metodo caratterizzato da un’ordine di convergenza 1 p 2; in particolare si pu dimostrare che p = √5+1
--2- 1,618 nel caso di radici semplici, ma solo lineare nel caso di radici multiple.

Implementazione

Costo computazionale

Questo metodo ad ogni ciclo iterativo svolge soltanto 5 operazioni floating point ed una valutazione di funzione, il costo per iterata quindi veramente molto basso.


PIC

Figura 2.8: Esempio di funzionamento del metodo delle secanti che ne mostra il comportamento nella determinazione dello zero della funzione x2 5x con punto di innesco x0 = 10.


Confronto tra il metodo di Newton, metodo di Newton nel caso di radici di molteplicit nota e metodo di accelerazione di Aitken

Queste tabelle mostrano il comportamento di questi tre metodi nella valutazione delle funzioni:

             10               10x
f1(x) = (x− 1) , f(x) = (x− 1) e

per valori decrescenti della tolleranza tolx e punto iniziale x0 = 10.








f1(x) f2(x)







tolx i x tolx i x







NTN







102 43 1.08727963568088 102 52 1.08957265092411







104 87 1.00084641497829 104 96 1.00087725269306







106 1311.00000820831010 106 1401.00000850818739







108 1741.00000008844671 108 1841.00000008251024







10102181.0000000008577310102271.00000000088907







NTM







102 1 1 102 4 1.00000041794715







104 1 1 104 5 1.00000000000002







106 1 1 106 5 1.00000000000002







108 1 1 108 6 1







1010 1 1 1010 6 1







ATK







102 2 1 102 6 0.99999998332890







104 2 1 104 7 1







106 2 1 106 7 1







108 2 1 108 7 1







1010 2 1 1010 7 1







Confronto tra il metodo di Newton, metodo delle corde, metodo delle secanti e metodo di bisezione

Questa tabella mostra il comportamento di questi metodi nella valutazione della funzione:

f(x) = x − cos(x)

per valori decrescenti della tolleranza tolx e punto iniziale x0 = 0, mentre per il metodo di bisezione consideriamo l’intervallo [0,1].




f(x)



tolx i x



NTN



102 3 0.73908513338528



104 3 0.73908513338528



106 4 0.73908513321516



108 4 0.73908513321516



1010 5 0.73908513321516



SCT



102 3.0.73911936191163



104 4 0.73908511212746



106 5 0.73908513321500



108 6 0.73908513321516



1010 6 0.73908513321516



CRD



102 120.73560474043635



104 240.73905479074692



106 350.73908552636192



108 470.73908513664657



1010590.73908513324511



BIS



102 7 0.734375



104 13 0.739013671875



106 200.73908424377441



108 250.73908513784409



1010310.73908513318747



Capitolo 3
Risoluzione di sistemi lineari

Un sistema lineare pu essere scritto nella forma:
Ax = b

dove A = (aij) m×n matrice dei coefficienti, b = (bi) m vettore dei termini noti e x = (xi) n vettore delle incognite. Assumeremo che la matrice abbia un numero di righe maggiore o uguale a quello delle colonne (m n) e che abbia rango massimo (rank(A) = n); queste assunzioni riescono comunque a coprire una buona parte dei problemi derivati dalle applicazioni.

3.1 Matrici quadrate

Iniziamo esaminando il caso in cui m = n, ovvero A una matrice quadrata. In quanto rank(A) = n, A una matrice nonsisngolare e la soluzione del sistame esiste ed unica:

     −1          − 1   − 1
x = A  b  con  AA   = A  A = I

3.1.1 Matrici diagonali

Analizziamo adesso il caso in cui A sia una matrice diagonale, ovvero:

(              )
  a11       0
|(      ..      |)
        .
   0       ann

in cui valgono le propriet:

In questo caso il sistema lineare Ax = b assume la forma:

a11x1  =.  b1
        ..
annxn  =  bn

e quindi la soluzione ottenibile:

     b
xi =-i-, i = 1,...,n.
    aii

ben posta in quanto come detto prima aii0 per i = 1,...,n. Il costo di questo metodo sia in termini di occupazione di memeoria che di quantit di operazioni floating-point lineare con la dimensione del problema; sono sufficiendi n operazioni floating-point ed un vettore di dimensione n per memorizzare gli elementi della diagonale della matrice A.

3.1.2 Matrici triangolari

Le matrici triangolari, che contengono informazione soltanto in una porzione triangolare, possono essere di due tipi:

Nel caso in cui la matrice A sia triangolare inferiore (il caso in cui sia triangolare superiore analogo) il sistema lineare della forma:

a11x1                             =   b1,
a21x1 + a22x2                     =   b2,
a31x1 + a32x2 + a33x3             =   b3,
..                         ..       ..  ..
.                           .      .  .
an1x1 + an2x2 + ⋅⋅⋅+ annxn          =   bn,

e quindi le soluzioni sono ottenibili tramie sostituzioni successive in avanti:

x   =   b∕a
x1  =   1(b −11a  x )∕a
x2  =   (b2− a21x1− a22 x )∕a
 3   .   3   31 1   32 2   33
     ..       ∑
xn  =   (bn −   nj=−11anjxj)∕ann

In quanto A nonsingolare, si ha che aii0, i = 1,...,n e quindi tutte le operazioni sopra descritte sono ben definite. In quanto la matrice contiene informazione solo in una sua porzione triangolare, sufficiente memorizzare tale porzione che contiene:

 n
∑     n(n-+-1)  n2
i=1 i =   2    ≈  2

posizioni di memoria. Per quanto riguarda il costo computazionale invece, da quanto detto sopra si ha che sono necessari 1flop per calcolare x1,3flop per calcolare x2,...,2n1flop per calcolare xn, e quindi:

∑n           2
   (2i− 1) = n flop.
i=1

3.1.3 Matrici ortogonali

Una matrice ortogonale se e solo se vale la propriet: A1 = AT. In questo caso quindi la soluzione del sistema lineare ottenibile facilmente calcolando x = ATb, al costo di un prodotto matrice-vettore; saranno richiesti pertanto circa 2n2flops ed n2 posizioni di memoria.

3.1.4 Propriet

Mostriamo adesso alcune propriet interessanti di questi tipi di matrici che andremo ad utilizzare con i metodi di fattorizzazione.

Teorema

La somma ed il prodotto di matrici triangolari inferiori (superiori) ancora una matrice triangolare inferiore (superiore).

Dimostrazione:

Dimostriamo il caso con matrici triangolari inferiori, l’altro si domostra in maniera analoga.
Siano L1,L2 matrici triangolari inferiori, ed indichiamo con (aij) e (bij i corrispettivi elementi.
Per ipotesi abbiamo che aij = bij = 0, j > i.
Allora gli elementi della matrice somma saranno composti:
L1 + L2 = (aij + bij) ma se j > i si ha (aij + bij) = 0 + 0 = 0, ovvero L1 + L2 triangolare inferiore.
Vediamo adesso il caso del prodotto.
Chiamaiamo cij gli elementi della matrice L1L2. Vogliamo dimostrare che cij = 0, se j > i.
L’elemento cij, quando j > i, scrivibile nella forma: eiT(L1L2)ej. Svolgendo le operazioni in quest’ordine: (eiTL1)(L2ej) si ha:

                  (   0 )
                  |   . |
(          n− i)  ||   .. ||
(a   a ...a ◜0◞.◟..◝0)  ||   0 ||
   i1 i2   ii       ||  bij ||
                  |(   ... |)
                     b
                      nj

ed il vettore b contiene 0 nelle prime j 1 posizioni. Ma j > i e quindi i primi j zero del vettore b annullano tutti i termini diversi da zero del vettore a, e cij = 0.

Teorema

Se L1 = (aij)eL2 = (bij) sono matrici triangolari inferiori, allora gli elementi diagonali della matrice L1L2 sono dati da (aiibii).
Dimostrazione:
Sia cij elemento della matrice risultante. Per definizione di prodotto tra matrici, si ha che

     ∑n
cij =   (aikbkj)
     k=1

Ricordando che aij = bij = 0 se j > i, si nota facilmente che la somma non interessata elementi che coinvolgono valori di k minori di j o maggiori di i in quanto il prodotto annullato dall’elemento rispettivamente della matrice L2 o L1; l’unico prodotto che non si annulla quando i = j (stiamo calcolando la diagonale) si ha con k = i = j e coinvolge appunto solamente i fattori akk e bkk. Cvd.

Teorema

La somma ed il prodotto di matrici triangolari inferiori (o superiori) ancora una matrice dello stesso tipo.
Dimostrazione:
Siano L ed S due matrici tiangolari inferiori. Se i < j quindi L(i,j) = S(i,j) = 0, e sommandole, si ha quindi che L + S(i,j) = L(i,j) + S(i,j) = 0 + 0 = 0, ovvero anche la loro somma una matrice triangolare inferiore. Nel caso del prodotto si ha che

            n
L ⋅S(i,j) = ∑  L(i,k) ⋅S (k,j)
           k=1

ma se i < j, L(i,j) = S(i,j) = 0 e quindi S(k,j) = 0 per k < j e L(i,k) = 0 per ,k > i e tutti i prodotti si annullano. Lo stesso ragionamento vale per matrici triangolari superiori, con i > j.

Lemma

Il prodotto di due matrici triangolari inferiori (o superiori) a diagonale unitaria ancora una matrice dello stesso tipo.
Dimostrazione:
Abbiamo appena dimostrato che il prodotto conserva la triangolarit, vediamo cosa succede alla diagonale:

          ∑n
L⋅S (i,i) =    L(i,k) ⋅S (k,i)
          k=1

ma in questo caso l’unico elemento non nullo nella sommatoria L(1,1) S(1,1,) = 1.

Lemma

Se L nonsingolare e triangolare allora L1 triangolare dello stesso tipo.
Dimostrazione:
Supponiamo L triangolare inferiore e L1 triangolare superiore. Allora si avrebbe che (LL1)(1,2) = L(1,1) L1(1,1) + 0 + 00, in quanto per ipotesi L(1,1),L1(1,1)0. Ne consegue risultato non pu essere la matrice identit, e l’inversa di una matrice triangolare deve pertanto essere anch’essa triangolare dello stesso tipo.

Implementazione

Questo codice implementa gli algoritmi di risoluzione dei sistemi lineari di tipo semplice, quali i triangolari inferiori e superiori, diagonale, ortogonale e bidiagonale superiore e inferiore.

3.2 Metodi di fattorizzazione

Da quanto descritto fino ad ora si pu concludere che siamo in grado di trattare in maniera efficiente matrici particolari quali quelle diagonali e quelle triangolari. Ovviamente non tutte le matrici sono di questi tipi, ma siamo in grado di definire alcuni metodi, detti di fattorizzazione, che permettodo di descrivere una matrice generica come prodotto di fattori che siano matrici particolari, ovvero del tipo a noi pi congeniale. Vogliamo quindi poter scrivere A = F1F2...Fk,dove Fk n×n sono matrici diagonali, triangolari oppure ortogonali. In questo modo, tenendo conto della fattorizzazione di A, possibile risolvere il sistema Ax = b nel seguente modo:

F1x1 = b, ,F2x2 = x1, ... ,Fkxk = xk−1, x ≡ xk

importante notare che tutti i sistemi appena citati sono di tipo semplice e quindi immediatamente risolvibili, e che non necessario memorizzare tutte le soluzioni intermedie, perch queste perdono di utilit una volta calcolata quella successiva.

3.3 Fattorizzazione LU

In base a quanto detto fin’ora, una situazione interessante si ha quando stiamo lavorando su di una matrice A che pu essere scritta come il prodotto di due matrici particolari che noi chiameremo Lower ed Upper, tali che L sia triangolare inferiore a diagonale unitaria e U triangolare superiore; ovvero la matrice A fattorizzabile LU.

Teorema

Se A una matrice nonsingolare e la fattorizzazione A = LU esiste , allora questa unica.
Dimostrazione:
Siano A = L1U1 = L2U2 due fattorizzazioni di A; si ha che:

0 ⁄= det(A) = det(L2U2) = det(L2)det(U2) = det(U2).

U2 quindi nonsingolare e si ha che L11L2 = U1U21 D. Ma L11L2 e U1U21 sono rispettivamente triangolare inferiore e triangolare superiore, pertanto D deve essere diagonale. Ma L11L2 una matrice a diagonale unitaria e questo implica che D sia la matrice Identit; ne consegue quindi che L1 = L2 e U1 = U2. cvd.
Rimane adesso da stabilire l’esistenza di tale fattorizzazione, e procederemo illustrando un metodo di fattorizzazione.
Sia v = (v1...vn)T n vettore tale che vk0, si supponga di volerne azzerare tutte le componenti successive alla k-esima esclusa, lasciando invariate tutte le altre, moltiplicandolo a sinistra con una matrice triangolare inferiore a diagonale unitaria L n×n. possibile definire il vettore elementare di Gauss come:

    1  ◜k◞◟◝
g = v-(0...0,vk+1,...,vn)T
     k

grazie al fatto di avere vk0, e quindi la matrice elementare di Gauss:

L ≡ I − geTk.

Ovvero:

        (  1                        )
        |     .                     |
        ||     ..                    ||
        ||           1               ||
    L = ||          vk+1  ...        || ,
        ||         − vk              ||
        ||           ...         ...    ||
        (           vn              )
    (     )        −vk            1
       v1
    ||   .. ||
    ||   . ||                    ◜-vk◞◟ ◝
Lv = || vk || = (I − geTk)v = Iv− g(eTkv) = v − gvk
    ||   0. ||
    (   .. )
        0

Questa matrice opera come ci eravamo preposti, infatti le prime k componenti del vettore Lv coincidono con le rispettive di v.
L’inversa di L facilmente ottenibile come: L1 = I + gekT, infatti:

L− 1L = (I + geTk )(I − geTk ) = I − geTk + geTk − g(eTkg)eTk = I.

in quando ekTg,k-esima componente del vettore g, nulla.
Procediamo adesso alla trasformazione della matrice nonsingolare A in una matrice triangolare superiore mediante la procedura iterativa denominata metodo di eliminazione di Gauss, che dura n 1 passi.
Definiamo:

          (   (1)       (1))
             a11  ⋅⋅⋅ a1n
A ≡ A (1) = |(  ...         ... |)
              (1)       (1)
             an1  ⋅⋅⋅ ann

la matrice al primo passo. L’indice superiore di ciascun elemento tiene traccia del passo pi recente in cui il valore dell’elemento stato modificato dalla procedura. Vogliamo quindi come primo passo definire una matrice triangolare inferiore a diagonale unitaria che sia in grado di trasformare A rendendo la prima colonna strutturalmente uguale alla relativa di una matrice triangolare superiore.
Se a11(1)0, possiamo definire il primo vettore elementare di Gauss:

g1 ≡ -1-(0,a(121),...,a(n1)1)T
     a(11)1

e la prima matrice elementare di Gauss:

              (    1             )
              |    a(1)           |
              ||  − -2(11)- 1        ||
L1 ≡ I − g eT= ||   a.11           ||
         1 1  ||    ..       ...    ||
              |(    a(1n)1           |)
                 − -(1)-        1
                   a11

per avere:

      (                    )
        a(1) ⋅⋅⋅  ⋅⋅⋅  a(1)
      ||  101  a(2) ⋅⋅⋅  a1n(2) ||
L1A = ||   .   2.2        2n.  || ≡ A(2).
      (   ..   ..         ..  )
         0   a(n22) ⋅⋅⋅  a(2n)n

Possiamo osservare come, per la struttura di L, la prima riga del prodotto L1A non venga modificata.
Generalizziamo quindi lo stesso procedimento al passo i-esimo, tenendo conto che ajj(j)0 per ogni j < i:

              (  (1)                         (1) )
              | a11  ⋅⋅⋅    ⋅⋅⋅    ⋅⋅⋅  ⋅⋅⋅   a1n  |
              ||  0   ...                     ...   ||
              ||  .   .     (i−1)             (i− 1) ||
Li−1...L2L1A  = ||  ..    ..  ai−1,i−1  ⋅⋅⋅  ⋅⋅⋅  ai−1,n || ≡ A (i)
              ||  ...          0     a(i)  ⋅⋅⋅   a(i)  ||
              ||  .           .     i.i        i.n  ||
              (  ..           ..     ..         ..   )
                 0   ⋅⋅⋅    0     a(ni)i  ⋅⋅⋅   a(ni)n

A questo punto solo se aii(i)0 sar possibile definire l’i-esimo vettore di Gauss:

     -1-      (i)      (i)T
gi ≡ a(i)(0.◟◝.◜.0◞,ai+1,i,...,ani),
      ii   i

e quindi la i-esima matrice elementare di Gauss:

               (                          )
               | 1  .                     |
               ||     ..                   ||
               ||           1              ||
Li ≡ I − g eT = ||         a(ii+)1,i- ..        ||
         i i   ||        −  a(iii)    .       ||
               ||           ...        ...   ||
               (           a(i)            )
                         − ani(iii)          1

e quindi:

                   (                                   )
                      a(11)1  ⋅⋅⋅  ⋅⋅⋅    ⋅⋅⋅    ⋅⋅⋅   a(11)n
                   ||       ..                      ..   ||
                   ||   0     .                     .   ||
   (i)              ||   ...   ...  a(iii)   ⋅⋅⋅    ⋅⋅⋅   a(ii)n  ||     (i+1)
LiA   = Li...L2L1A = ||   ..             (i+1)         (i+1) || ≡ A
                   ||   .        0   ai+1,i+1  ⋅⋅⋅  ai+1,n ||
                   |(   ...         ...     ...           ...   |)
                       0   ⋅⋅⋅  0    a(i+1)   ⋅⋅⋅  a(i+1)
                                      n,i+1        nn

Se si osserva che nella i-esima matrice di Gauss (Li I gieiT)) le prime i righe corrispondono con quelle della matrice identit, si pu notare che le prime i righe delle matrici A(i) e A(i+1) coincidono e gli elementi azzerati al passo i-esimo (che stanno nelle prime i 1 colonne) non vengono alterati dalla moltiplicazione per Li.
Se risulta possibile iterare i passi fin qui descritti i = n 1 volte, si ottiene che:

                       (                          )
                         a(111) ⋅⋅⋅    ⋅⋅⋅     a(11)n
                       ||       ..              ..   ||
A (n) ≡ U ≡ Ln−1...L1A = ||        .   (n−1)     (.n− 1) ||
                       (           an−1,n−1  an−1,n )
                                             a(nnn)

Adesso, osservando che la matrice Ln1...L1 il prodotto di matrici triangolari inferiori a diagonale unitaria, anch’essa di questo tipo, e cos sar anche la sua inversa. Tenendo conto che:

 T
ek gi ≡ gki = 0, per k ≤ i,

e ponendo formalmente la matrice Ln1...L1 uguale a L1:

                −1    −1   −1
L  =   (Ln−1...L1T )  = L1 ...LTn−1         T           T
   =   ((I + g1e1)...(I + gn−1en−1)) = I + g1e1 + ...+ gn−1en−1
          1
       || g12   1            ||
   =   |(  ...   ...   ...     |)
         gn1  ⋅⋅⋅ gn,n−1  1

Siamo quindi riusciti ad ottenere due matrici L ed U tali che A = LU.
Consideriamo adesso le condizioni sufficienti a garantire l’esistenza della fattorizzazione LU: il metodo di eliminazione di Gauss definito se e solo se possibile costruire la matrice di Gauss al generico passo i-esimo.

Lemma

Se A una matrice nonsingolare, la fattorizzazione A = LU definita se e soltanto se aii(i)0, i = 1,2,...,n, ovvero se e solo se la matrice U nonsingolare.

Definizione

Si dice sottomatrice principale di ordine k di una generica matrice A n×n l’intersezione tra le sue prime k righe e k colonne, ovvero:

     (               )
        a11 ⋅⋅⋅  a1k
Ak = |(   ...        ...  |)
        a   ⋅⋅⋅  a
         k1       kk

In particolare si ha che il determinante det(Ak) il minore principale di ordine k di A.

Lemma

Una matrice triangolare nonsingolare se e solo se tutti i suoi minori principali sono non nulli.

Lemma

Il minore di ordinde k di A in A = LU coincide con il minore di ordine k di U in A(n) U

Teorema

Se A una matrice nonsingolare, la fattorizzazione A = LU esiste se e solo se tutti i minori principali di A sono non nulli.

La fattorizzazione quindi esiste se e solo se U nonsisingolare, e questo equivale a richiedere che tutti i minori principali di U siano non nulli, e quindi anche tutti quelli di A devono essere non nulli.

3.3.1 Costo computazionale

Esaminiamo adesso il costo dell’algoritmo di eliminazione di Gauss in termini di occupazione di memoria.
Strutturalmente, al passo i-esimo l’algoritmo azzera, nella matrice A(1), le componenti i + 1,...,n della colonna i; al contempo si ottiene l’i-esimo vettore di Gauss che contiene informazioni soltanto nelle ultime n i componenti. Possiamo quindi pensare di sovrascrivere A durante il processo, ottenendo una matrice che contiene la parte significativa della matrice U nella porzione triangolare superiore e le componenti significative della matrice L (ovvero gli elementi significativi dei vettori di Gauss) nella parte strettamente inferiore, tenendo conto che la diagonale della matrice L ha tutti elementi pari a 1 e non necessita quindi di essere memorizzata esplicitamente. Possiamo quindi concludere che il metodo di eliminazione di Gauss non richiede memoria addizionale alla matrice di partenza A, che quanto di meglio si possa richiedere.
Per quanto riguarda invece il costo computazionale, ad ogni passo i dell’algoritmo necessario calcolare:

Considerando che le prime i componenti del vettore gi sono nulle, che il vettore eiTA(i) corrisponde alla i-esima riga di A(i) che ha le prime i componenti nulle e che sappiamo a priori che le componenti da (i + 1)an della colonna i della matrice A(i+1) sono nulle, si nota che solo la sottomatrice delimitata da (i + 1,i + 1),(n,n) verr modificata al passo i. Il costo computazionale di una iterazione consiste in n i divisioni per calcolare il vettore di Gauss, (n i)2 sottrazioni ed altrettante moltiplicazioni per aggiornare la matrice; in totale la fattorizzazione avr un costo approssimativo di:

 n
∑  (2(n − i)2 + (n − i)) ≈ 2n3.
i=1                    3

Implementazione

Questo codice implementa l’algoritmo di fattorizzazione LU:

3.4 Matrici a diagonale dominante

Esistono delle classi di matrici per cui vero che se A nonsingolare allora lo anche Ak, k = 1,...n, e sono le matrici a diagonale dominante.

Definizione

La matrice A n×n si dice a diagonale dominante per righe se:

|a | > ∑ |a |,  i = 1,...,n;
  ii   j⁄=i  ij

e a diagonale dominante per colonne se:

|aii| > ∑ |aji|,  i = 1,...,n;
      j⁄=i

Per questa classe di matrici valgono le seguenti propriet:

Lemma

Se A una matrice a diagonale dominante per righe (o per colonne), allora anche tutte le sue sottomatrici principali sono dello stesso tipo

Lemma

Una matrice A a diagonale dominante per righe (o per colonne) se e solo se AT a diagonale dominante per colonne (per righe).

Teorema

Se una matrice A diagonale dominante (per righe o per colonne), allora A non singolare.
Dimostrazione:
Sia A una matrice a diagonale dominante per righe, supponiamo per assurdo che A sia singolare. Allora deve esistere x0 tale che Ax = 0. Assumiamo

                         x
xk = max (|xi|) = 1, e quindi |-i| ≤ 1 ∀i = 1,...,n.
                         xk

Ma se Ax = 0 allora ekTAx = ekT0 = 0, ne consegue che:

n∑                             ∑n
   akjxj = 0 ⇒ akkxk = akk = −      akjxj =
j=i                          j=1,j⁄=k

         |       |
         ||∑       ||  ∑          ∑
= |akk| = || akjxj||≤    |akjxj| ≤   |akj|
         |j⁄=k     |  j⁄=k        j⁄=k

ma allora A non diagonale dominante. Assurdo.

Lemma

Se A una matrice a diagonale dominante, allora A fattorizzabile LU.

3.5 Fattorizzazione LDLT

Un’altra classe di matrici fattorizzabili LU quella delle matrici simmetriche e definite positive (sdp).

Definizione

Una matrice A n×n sdp se simmetrica (ovvero A = AT) e se vale:

      n          T
∀x ∈ ℝ ,x ⁄= 0 : x Ax > 0.

Lemma

Se A sdp allora aii > 0, i = 1,...,n.
Dimostrazione:
Infati, se A sdp si dee avere che aii = eiTAei > 0. Cvd.

Lemma

Se A sdp allora A nonsingolare.
Dimostrazione:
Supponiamo per assurdo che A sia contemporaneamente sdp e singolare. Essendo A singolare deve esistere x0 tale che Ax = 0, ovvero xTAx = 0 ovvero A non sdp. Cvd.

Teorema

A sdp se e solo se k = 1,...,n, Ak sdp.
Dimostrazione:
Procediamo per assurdo, supponendo Ak non sdp. Allora deve esistere y k,y0, tale che yTAky 0. Allora deve esistere anche x n tale che x = (y,0..0)T0, che soddisfi 0 < xTAx in quanto A sdp per ipotesi. Ma allora avremmo:

         n−k         |
 T       ◜◞◟◝ T (-Ak---BT--)(  y )    T
x Ax = (y0...0)     B  |D       0   = y Aky ≤ 0

che va contro l’ipotesi.

Corollario

Se A sdp allora A fattorizzabile LU
Dimostrazione:
Da quanto dimostrato fino ad ora sappiamo che tutte le sottomatrici principali di A sono sdp, e quindi anche i corrispondenti minori principali sono tutti non nulli. Cvd.

Teorema

Gli elementi diagonali di una matrice simmetrica e definita positiva sono positivi.
Dimostrazione:
Se A sdp, per definizione si ha che aii = eiTAei > 0.

Teorema

Una matrice A sdp se e solo se scrivibile come:

A = LDLT

dove L triangolare inferiore a diagonale unitaria e D diagonale con elementi diagonali positivi.
Dimostrazione:
Per prima cosa dimostriamo che se A scrivibile come LDLT allora sdp. Dobbiamo quindi dimostrare che x0,xTLDLTx > 0. Possiamo considerare il vettore y = LTx e quindi yT = xTL, sicuramente non nullo in quanto L una matrice non singolare e x0 per ipotesi. Si pu quindi scrivere:

        n
 T     ∑     2
y Dy = i=1diiy > 0.

Adesso possiamo dimostrare che se A sdp allora A = LDLT. Abbiamo dimostrato precedentemente che se A sdp allora fattorizzabile LU. Possiamo pensare U come D dove:

                   (               )
                      U11       0
D = diag(U11...Unn) = |(      ...      |)
                       0       Unn

ed

            (     )                    (  1      ∗)
 ˆ    −1     --Uij-                 ˆ   |     .    |
U = D   U =    Uii   ,  ∀ij = i,...,n U = (     ..   )
                                          0      1

ovvero D una matrice diagonale e triangolare superiore a diagonale unitaria. Osserviamo che

AT = (LDLT )T = (LT)TDT LT = LDLT  = A

e che:

LU = LD ˆU = A = AT = (LD ˆU )T = ˆUT (DT LT )

ma quindi:

L(DUˆ) = UˆTDLT ⇒ L = ˆUT ⇒  ˆU = LT.

Per l’unicit della fattorizzazione LU, in quanto T triangolare inferiore a diagonale unitaria e DLT triangolare superiore, si ha che T = L e quindi la fattorizzazione A = LDLT ben definita. Dimostriamo adesso che gli elementi diagonali di D sono positivi, facendo vedere che questa sdp. Banalmente D simmetrica in quanto diagonale, e, qualunque si fissi x0 esiste un unico vettore y0 che soddisfi LTy = x. Si ha quindi che:

xT Dx = (LT y)TD(LT y) = yTLDLT y = yTAy > 0

in quanto A sdp.

3.5.1 Costo computazionale

Trattiamo innanzitutto gli aspetti che riguardano la memorizzazione delle informazioni. Per quanto riguarda la memorizzazione di matrici strettamente triangolari, dobbiamo pensare di memorizzare solo la parte che contiene informazioni, evitando di sprecare memoria per la parte strettamente triangolare nulla. Questo pu essere fatto memorizzando la parte triangolare in un vettore e di utilizzare un piccolo programma che permetta di recuperare l’elemento ij della matrice triangolare dal relativo vettore. Nel caso quindi della fattorizzazione LDLT, in quanto stiamo trattando matrici sdp (e quindi triangolari) possiamo pensare di memorizzarne solo la parte triangolare attraverso il codice appena visto, ed in pi possiamo pensare di riscrivere la matrice di partenza con la porzione strettamente triangolare di L (o LT) e la diagonale D, mantenendo quindi il costo di memorizzazione inferiore ad n2.
Analizziamo adesso il costo in termidi di operazioni floating-point.
Stiamo lavorando su una matrice A = (aij) = LDLT, una matrice L = (lij che ha le propriet di avere lii = 1 ed lij = 0 se j > i e D diagonale. Si ha quindi che:

                         T       T     T
∀i ≥ j, aij  =           ei AejT = ei LDTL T ej
            =              (ei L )D (ej L)   (     )
                                           |  lj.1 |
                           (  d1         ) ||  ..  ||
                           |      .      | ||  ljj ||
            =   (e◟11..◝◜.eii◞0...0)(       ..     ) ||  0  ||
                   i                  dn   |(  ..  |)
                                              .
                           ∑j                 0
            =                k=1likljkdk

e quindi aij = k=ijlikljkdk, i = j,...,n;j = 1,...,n.
Si posso verificare due casi, i = j e i > j:

              ∑                 ∑
i = j    ajj =  jk=1l2jkdk = l2jjdj +  jk−=11l2jkdk
        ma ljj = 1 quindi dj = ajj − ∑jk−=11 ljkljkdk,
                     j = 1,...,n.

i > j aij = ∑jk=1likljkdk = lijljjdj + ∑jk−=11likljkdk
            aij − ∑j−1 likljkdk
       lij = -------k=d1-------,  i = j + 1,...,n.
                    j

Si ha quindi che per effettuare il calcolo della colonna k-esima necessaria soltanto l’informazione relativa a quelle precedenti. e che l’elemento aij sia necessario solo per il calcolo della colonna j-esima; questo ci lascia liberi di sovrascriverlo con le nuove informazioni derivanti dal calcolo appena effettuato. In particolare possiamo memorizzare le informazioni relative agli elementi dk sulla diagonale e i fattori lij nella parte strettamente triangolare inferiore della matrice A. Si deve inoltre tenere conto che ljkdk viene calcolato nel ciclo esterno e sar poi necessario nel ciclo interno, quindi opportuno appoggiarsi ad un array temporaneo per memorizzare il risultato di questa operazione; a questo scopo sufficiente un array di j 1 locazioni di memoria in quanto l’ultimo elemento diagonale ad essere calcolato potr essere memorizzato direttamente sulla matrice A.

Implementazione

Questo codice implementa l’algoritmo di fattorizzazione LDLT descritto fino ad ora:

3.6 Pivoting

Prendiamo adesso in esame il caso in cui la miatrice A che vogliamo fattorizzare sia nonsingolare ma che non abbia la propriet di avere tutti i minori principali non nulli. Abbiamo prima dimostrato che in questo caso la fattorizzazione A = LU non esiste; possiamo per definire una fattorizzazione che esista sotto la sola ipotesi della nonsingolarit di A.
Supponiamo di star procedendo al primo passo del metodo di eliminazione di Gauss: per poter procedere necessario che a110. Supponiamo per che a11 = 0: in quanto A nonsingolare deve esistere, sulla prima colonna, un elemento non nullo. Possiamo quindi scegliere:

|a(1k1)1| ≡ max (|a(k11)|) > 0, k = 1,...,n.

Definiamo una matrice elementare di permutazione del tipo:

                                 (   eTk  )
                     |           ||   eT2  ||
     (  0   0T    1  |       )   ||    ..  ||
     ||  0  Ik1−2  0  |       ||   ||   T.  ||
P1 ≡ ||  1   0T    0  |       || ≡ ||  ek−T1 ||
     (---------------|-------)   ||   eT1  ||
                     |In − k1    ||  ek+1 ||
                                 (    ...  )
                                     eTn

e non altro che la matrice identit di ordine n con le righe 1 e k1 scambiate. Notiamo che P1 una matrice simmetrica ed ortogonale, ovvero vale: P1 = P1T = P11.
Segue quindi che P1A(1) non altro che la matrice A(1) con le righe 1 e k1 scambiate. Adesso abbiamo reso possibile l’effettuazione del primo passo dell’algoritmo di eliminazione di Gauss, abbiamo infatti reso non nullo a11:

g = --1-(0,a(1),...,a(1),...,a(1))T,
 1  a(k1)1    21     11    n1
      1

e la prima matrice elementare di Gauss:

L1 = I − g1eT1

che consente quindi di ottenere:

          (   (1)            (1) )
             ak11  ⋅⋅⋅  ⋅⋅⋅  ak1n
      (1)  ||   0    a(22)2  ⋅⋅⋅  a(22n) ||    (2)
L1P1A   = ||    ..    ..         ..  || ≡ A
          (    .    .(2)       (.2) )
              0    an2  ⋅⋅⋅  ann

impotante notare che:

det(A(2)) = det(L1P1A ) = d◟et◝(L◜1◞)d◟et(◝P◜1)◞ d◟et◝(◜A)◞ ⁄= 0
                        =1     1  ⁄=0peripotesi

Ma A(2) della forma:

(      |   )
  a(k11) |rT2
 --0---|A2-

e quindi:

det(A (2)) = det(a(1))det(A2) ⁄= 0 ⇒ det(A2 ) ⁄= 0.
              k1

Possiamo quindi procedere in modo analogo al metodo di fattorizzazione LU, avendo al passo i-esimo la matrice:

                  (   (1)                           (1)  )
                     ak11  ⋅⋅⋅    ⋅⋅⋅    ⋅⋅⋅ ⋅⋅⋅   ak1n
                  ||   0   ...                      ...   ||
                  ||   .   .                            ||
                  ||   ..    ..  a(iki−−11),i−1  ⋅⋅⋅ ⋅⋅⋅  a(ik−i−11),n ||     (i)
Li− 1Pi−1 ⋅⋅⋅L1P1A = ||   ..                  (i)        (i)  ||  ≡ A
                  ||   ..           0.     aii.  ⋅⋅⋅   ai.n  ||
                  |(   ..           ..      ..         ..   |)
                      0   ⋅⋅⋅     0     a(i) ⋅⋅⋅   a(ni)n
                                         n1

Avremo quindi akii(i)0 se A nonsingolare. Definiamo quindi la matrice Pi che permuta le righe i e ki (con ki i ):

            |               |
    ( -Ii−1------------------------)
    |       |               |      |
    ||       | 0    0T    1  |      ||
Pi ≡ ||      | 0  Iki−i−1  0  |      ||
    ||       | 1    0T    0  |      ||
    |(       |               |      |)
      ------|---------------|I-----
                              n−ki

e si avr quindi che l’elemento (i,i) della matrice PiA(i) diventa akij(i). Definiamo quindi l’i-esimo vettore di Gauss:

     1         (i)     (i)    (i)T
gi =-(i)(0◟,..◝◜,0◞,ai+1,i,...,aii ,...,ani) ,
    akii   i

e quindi l’i-esima matrice elementare di Gauss:

Li = I − gieTi

tali che:

                        (  a(1)  ⋅⋅⋅  ⋅⋅⋅     ⋅⋅⋅    ⋅⋅⋅  a(1) )
                        |   k11  .                        k1.n |
                        ||   0    ..                       ..  ||
                        ||   ...   ...  a(i)     ⋅⋅⋅    ⋅⋅⋅  a(i)  ||
LiPiA(i) = LiPi⋅⋅⋅L1P1A = ||            kKii                ki,n ||  ≡ A(i+1)
                        ||   ...         0    a(ii++11,)i+1  ⋅⋅⋅  a(ii++11,n)||
                        ||   ..         ..      ..            ..  ||
                        (   .         .      .            .  )
                            0   ⋅⋅⋅   0    a(ni,+i+11)  ⋅⋅⋅  a(in+n1)

Con questo procedimento, se A nonsingolare sar sempre possibile iterare i passi appena descritti fino a i = n 1, ottenendo la fattorizzazione:

Ln−1Pn−1Ln −1Pn−1...L1P1A  = A(n) ≡ U.

Se si tiene conto che le matrici di permutazione elementari Pi sono sia simmetriche che ortogonali e se moltiplicate per un vettore sortiscono l’effetto di permutare le componenti i e ki con ki i, si pu pensare di riscrivere:

  (n)       ˆ   ˆ     ˆ
A    ≡ U = Ln−1Ln−2...L1PA

dove

Lˆn −1 ≡   Ln−1,
   ˆLi ≡   Pn−1...Pi+1LiPi+1...Pn−1,  i = 1,...,n − 1,
   P  ≡   Pn−1...P1

P una matrice di permutazione, ovvero una matrice ortogonale che se moltiplicata per un vettore ne permuta le componenti. Per come abbiamo definito Pi:

eTPj = eT, con j > i,
 i      i

e quindi:

ˆLi = Pn −1...Pi+1(I − gieTi )Pi+1...Pn−1
  =  I − (Pn −1...Pi+1gi)(eTi Pi+1...Pn−1) ≡ I − ˆgieTi .

Il vettore ˆgi ha la stessa struttura del vettore elementare di Gauss prima definito ma con le ultime ni componenti permutate tra loro. Ne consegue che la struttura della matrice ˆLi analoga a quella della matrice elementare di Gauss prima definita; si pu quindi concludere che la matrice ˆLn1...ˆL1 L1 triangolare inferiore a matrice unitaria. Abbiamo quindi ottenuto una fattorizzazione con pivotin parziale:

P A = LU

Teorema

Se A una matrice nonsingolare, allora esiste una matrice di permutazione P tale che PA fattorizzabile LU.

Osservazione:

Questo metodo permette di estendere l’utilizzo della fattorizzazione LU, infatti si ha che:

Ax = b ⇔ P Ax = P b ⇔ (Ly = Pb ∧   U x = y).

e la matrice di permutazione P serve soltanto per permutare il vettore dei termini noti b, operazione che precede direttamente la risoluzione dei sistemi triangiolari sopra descritti nell’ordine specificato.

Implementazione

Questo algoritmo implementa la fattorizzazione LU con pivoting appena descritta:

3.7 Condizionamento del problema

Studieremo adesso in che modo delle perturbazioni sui dati di un sistema lineare del tipo Ax = b si ripercuotono sulla sua soluzione.

Dobbiamo prima dare delle definizioni:

Norma

Sia x n, la norma di un vettore una funzione definita su uno spazio vettoriale V e codominio in tale che:

1.
x V : x∥≥ 0, x= 0 x = 0.
2.
x V e α : αx= |α|⋅∥x.
3.
x,y V, x + y∥≤∥x+ y

Norma “p”

Sia V n. Definiamo:

      (∑n     ) 1p
∥x∥p =     |xi|p   ,  p > 0.
        i=1

Piu in dettaglio:

   (                             ∑
   {  = 1      Norma 1    ∥x∥1 =   ni=1|xi|
p =   = 2   Norma euclidea ∥x∥2 = (∑n   |xi|2)12
   (  = ∞    Norma infinito  limp →∞ ∥xi∥=p1 ⇒  max (|xi|),  i = 1,...,n.

Verifichiamo adesso che le norme appena descritte verificano le tre propriet sopra elencate.
Per quanto riguarda la prima propriet si ha che:

Per la seconda si ha che:

E infine per la terza:

Notiamo inoltre che le norme sono tra loro equivalenti, vale infatti che:

∀p,q,  ∃  γ1,γ2  |  γ1,γ2 > 0, γ1 < γ2, γ1∥x∥p ≤ ∥x∥q ≤ γ2∥x∥p

Norma indotta su matrici

Sia V = m×n e A m×n. Definiamo la norma “p” su una matrice indotta dalla corrispondente norma “p” su vettore come:

∥A ∥p = m∥axx∥ ∥Ax∥p
          p

dove A e xpossono essere su spazi diversi se A rettangolare.

Andiamo a studiare il sistema lineare:

(A + ΔA )(x+ Δx ) = b + Δb

dove le perturbazioni ΔA e Δb determinano la perturbazione Δx sulla soluzione. Consideriamo, per semplicit, il caso in cui i dati dipendano da un parametro scalare di perturbazione di dimensioni prossime a zero:

                     n×n
A (ε) = A+ εF,  F ∈ ℝ    =⇒  ΔA = εF,

b(ε) = b+ εf,f ∈ ℝn =⇒  Δb = εf.

Da cui consegue:

A (ε)x(ε) = b(ε),

Si osservi che:

A (0) = A, b(0) = b, =⇒  x(0) = x.

Sviluppando in serie di Taylor in ε = 0, si ottiene:

x(ε) = x + εx′(0)+ O(ε2) ≈ x+ εx′(0)

cio

Δx ≡ x(ε) − x ≈ εx′(0)

Si pu inoltre ricavare che:

  ′       ′      ′
A (0)x+ Ax (0) = b(0).

Quindi:

d-A(ε)x(ε) = d-b(ε) ≡ f
dε          dε

    ′            ′     ||            ′   ||                 ′
= A (ε)x(ε) +A (ε)x (ε) =  F x(ε)+ A (ε)x (ε)  ε=0 = Fx + A(0)+ x(0) = f

che permette di ottenere:

x′(0) = A−1(f− Fx).⇒ εx′(0) = A−1(εf− εFx) ⇒ Δx ≈ A− 1(Δb − ΔAx ).

Passando alle norme:

         −1
∥Δx∥ = ∥A  (Δb − ΔAx )∥

Ma per la disugualianza triangolare si pu affermare che:

∥Δx ∥ ≤ ∥A −1∥⋅∥Δb − ΔAx ∥ ≤ ∥A −1∥(∥Δb ∥+ ∥ΔA ∥ ⋅∥x∥)

            (             )              (                )
∥Δx∥-    −1   ∥Δb-∥                  −1   --∥Δb∥--  ∥ΔA-∥
∥x∥  ≤ ∥A  ∥   ∥x∥ + ∥ΔA ∥  = ∥A∥ ⋅∥A  ∥  ∥x∥⋅∥A ∥ + ∥A ∥

ma tenendo condo che:

b = Ax,  ∥b ≤ ∥A∥ ⋅∥x∥;

si pu concludere che:

           (             )
       − 1  ∥Δb-∥   ∥ΔA-∥
∥A ∥⋅∥A  ∥   ∥b∥ +  ∥A ∥

e quindi:

                          (                )
    ∥Δx ∥                 |  ∥Δb ∥   ∥ΔA ∥ |
    -----     ≤ ∥A∥ ⋅∥A−1∥||  ----- + ----- ||
    ◟∥x◝∥◜ ◞                (  ◟∥b∥--◝◜-∥A-∥◞ )
erroresuidatiinuscita              erroresuidatiiningresso

Abbiamo ottenuto una equazione che misura una sorta di errore relativo sul risultato, ed in particolare abbiamo che k = A∥⋅∥A1 il numero di condizionamento del problema. importante notare che k(A) = A∥⋅∥A1∥≥∥A A1= I= 1.

3.8 Sistemi lineari sovradeterminati

Poniamo adesso il caso di voler risolvere un sistema di equazioni lineari in cui ci siano pi equazioni che incognite ma in cui la matrice dei coefficienti abbia rango massimo. Formalmente siamo quindi nel caso:

A ∈ ℝm ×n,
x ∈ ℝn,  b ∈ ℝm,
n = rank(A).

In particolare possiamo considerare la matrice A in questo modo:

  m×n                         m
ℝ     ∋ A = (C1C2...Cn),  Cj ∈ R

si ha che:

         (       ||    n      )
span(A ) = { y ∈ ℝm ||y = ∑ α C }
         (       ||    j=1  j j)

dove

    (     )
       α1
x = |(  ..  |) ⇒ y = Ax
       .
      αn

ed

n = rank (A ) = dim(span(A))

Possiamo quindi affermare che:

Ax = b ha soluzione ⇐⇒ b ∈ span(A).

Definizione: spazio nullo
null(A) = {x ∈ ℝn |Ax = 0}

Se la matrice ha rango massimo null(A) contiene un solo vettore; generalmente si ha:

dim(null(A )) = n− rank(A )

se il rango della matrice non massimo.
In particolare avremo che:

             {
n− rank(A) =   = 0,rank(A) = n  ⇐⇒  null(A) = {0}
               > 0,< n          =⇒ ∃ infk vettori in null(A).

Dobbiamo notare che:

∀v ∈ null(A ), A (x + v) = Ax + A◟◝v◜◞ = b
                            =0

ovvero se rank(A) non massimo ci sono infinite soluzioni.

Problema dell’esistenza della soluzione

Sia b m un vettore e span(A) = n < m il rango della matrice. Allora si ha che span(A) m. Questo significa che l’insieme differenza (m span(A)) molto grande e la situazione b span(A) risulta essere un evento piu “fortunato” che probabile, ovvero in generale non si pu assumere che b sia nello spazio di A.

Definizione: Vettore residuo

Per risolvere il sistema sopra descritto proveremo dunque a trovare il vettore x che minimizza il vettore:

   (     )
      r1
r = |( ..  |) ≡ Ax − b
      .
      rm

Ovvero andremo a cercare la suluzione del sistema nel senso dei minimi quadrati:

∑m
    |ri|2 = ∥r∥22 = ∥Ax − b∥22
 i=1

Teorema

Se A una matrice in m×n,m > n,rank(A) = n allora esistono Q m×m,QTQ = Im e Rˆ n×n triangolare superiore e nonsingolare tali che:

            ( Rˆ )        m ×n
A = QR  ≡ Q   0   ,  R ∈ ℝ

importante ricordare la seguente propriet delle matrici ortogonali:

∥Qv∥2= (Qv )T Qv = vT QTQv = vT v = ∥v∥2
    2                               2

ovvero la norma euclidea invariante per moltiplicazioni per matrici ortogonali a sinistra o a destra.
Vediamo su quali fattori lavorare per minimizzare il vettore residuo:

∥Ax − b∥22 = ∥QRx − b∥22 = ∥Q(Rx − QTb)∥22 = ∥Rx − g∥22
  ∥∥( Rˆ )     ∥∥2  ∥∥(  ˆR )    ( g  ) ∥∥2  ∥∥(  ˆRx− g  ) ∥∥2
= ∥∥   0  x − g∥∥ = ∥∥   0   x−   g1   ∥∥ = ∥∥    − g 1   ∥∥
  ∥∥       ∥∥2   2                2    2          2     2
= ∥ˆRx − g1∥2 + ∥g2∥22

in cui abbiamo definito e quindi utilizzato:

          (    )
g = QT b =   g1   tale che g1 ∈ ℝn, g2 ∈ ℝm −n.
             g2

Generalmente quindi b non appartiene allo span(A). Ne deriva quindi che se ˆRx = g1 allora Ax b22 = g222 min. importante osservare che:

3.8.1 Esistenza della fattorizzazione QR

Matrice elementare di Householder

Dato il vettore:

z = (z1,...,zm )T ∈ ℝm,  z ⁄= 0

vogliamo determinare una matrice ortogonale H tale che:

Hz = αe ,  α ∈ ℝ
        1

Si ha che:

∥z∥2 = zTz = zTHT Hz = α2eTe = α2,
  2                      1 1

ovvero possiamo definire α come:

α =  ∥z∥2.

Consideriamo quindi una matrice della seguente forma:

         2
H = I − -T--vvT,  v ⁄= 0,
        v◟◝v◜◞
        =∥v∥22

dove v m sar diverso dal vettore nullo se anche z lo . La matrice H cos costruita simmetrica per definizione, ed anche ortogonale:

  T        2  (    -2-- T) (    -2-- T )      -4-- T   --4---  T 2
H  H   = H  =  I − vTvvv     I − vTvvv  = I − vTvvv  + (vTv)2(vv )
            --4-  T  -----4----(   T   T)      -4-- T   -4-- T
       = I − vT vvv + (vTv)(vTv) v(v v)v = I − vTvvv  + vTvvv  = I

Scegliamo adesso un vettore

v = z − αe1

e verifichiamo se Hz = αe1, α .

      (     2     )         2
Hz  =   I −-T--vvT  z = z−-T--vvTz
    (tenendvo covnto che vTz = v(z−vαe )z = zT z− αe Tz = ∥z∥2 − αz )
           2  (        )        12  (        )1       2     1
    = z( − vTv-zT z− αz1 v =) z − vTv-zT z− αz1 (z− αe1)
             2 (        )        2  (        )
    =   1− vT-v zTz− αz1   z+ α vTv- zT z− αz1 e1
           ◟-----◝◜-----◞       ◟-----◝◜-----◞

Dall’equazione sopra scritta notiamo che se gli scalari sopra evidenziati fossero di valore pari ad 1 allora otterremmo Hz = αei, dimostriamolo:

 2                 2∥z∥2− 2αz1
vTv-(zTz − αz1) = (z−-αe-2)T(z−-αe-) = 1  ⇒   2(zTz− αz1) = vTv
                      1        1

dove svolgendo vTv otteniamo

vTv = (z− αe1)T(z− αe1) = zT z+ α2− 2α eT1 z considerando  α2 = ∥z∥2 = zTz
                                    ◟◝z◜◞
                                      1

       T    T     T           T
allora  v v = z z+ z z− 2αz1 = 2z z− 2αz1

per cui abbiamo verificato

   T            T
2(z z − αz1) = 2z z− 2αz1

abbiamo quindi dimostrato che Hz = αe1.

Discussione del segno di α
            (        )
               z1 − α
            ||   z2   ||
v = z − αe1 = |(  ...   |)
                z
                 n

tenendo conto del condizionamento della somma algebrica imponiamo:

α = − sign (z1)⋅∥z∥2

Il calcolo di v avr pertanto numero di condizionamento pari a uno.
Vediamo adesso come si possa risparmiare una locazione di memoria nella memorizzazione di v:

|v1| > 0, se∥z∥ > 0

possiamo quindi riscrivere v come:

      (     )
         1v2
      ||  v1 ||      ˆ   ˆ   1-
v = v1|(   ... |) = v1 ⋅v ⇒ v = v1v
         vm
         v1

Notiamo adesso che:

       (           )
H (v) =  I −-2-vvT
            vTv

e quindi

                    (             )
     (    -2-- T)   |    --2--1  T|    (    -2-- T)
H (ˆv) I − ˆvTˆvˆvˆv   = |(I −  vT vv2vv |)  =  I − vTvvv   = H (v)
                          v2--1
                           1

ovvero si pu sostituire il vettore v con un qualunque altro vettore multiplo scalare di esso e la matrice di Householder non cambia. Possiamo quindi usare il vettore normalizzato e risparmiare una locazione di memoria.

Esistenza della fattorizzazione QR

Dimostriamo adesso in maniera costruttiva l’esistenza della fattorizzazione QR:
Sia:

    (   (0)       (0))
       a11   ⋅⋅⋅  a1n
A = |(   ...        ...  |)  ≡ A(0).
      a(0)  ⋅⋅⋅  a(0)
       m1        mn

dove l’indice superiore indica il passo pi recente in cui stato modoficato l’elemento corrispondente.
Considerando la prima colonna della matrice A(0) possiamo definire la matrice di Householder H1 m×m tale che:

              (      )
   (  a(0) )      a(111)
   |  ⋅11⋅⋅ |   ||   0  ||
H1 (   (0) ) ≡ (  ⋅⋅⋅ )
      am1         0

e se A ha rango massimo si deve necessariamente avere che a11(1)0, altrimenti la prima colonne della matrice risulterebbe nulla. Otterremo quindi:

         (  (1)  (1)       (1))
         | a11  a1(21)  ⋅⋅⋅  a1(n1)|
    (0)  ||   0  a21   ⋅⋅⋅  a2n ||     (1)
H1A    = |(   ...    ...        ...  |)  ≡ A
                 (1)       (1)
             0  am2   ⋅⋅⋅  amn

Possiamo adesso iterare il passo fatto in precedenza sulla seconda colonna della matrice A(1) a partire dalla seconda componente e definire quindi una matrice di Householder in questo modo:

                               (   (1) )   (  a(1) )
     ( 1 |     )                  a22     |   202 |
H2 ≡  ---|--(2)- ∈ ℝm ×m,  H (2)|(   ...  |) ≡ ||   ..  ||
         |H                       a(1)     (   .  )
                                   m2         0

La seconda matrice di Householder cos ottenuta quindi la composizione della matrice identit con una sottomatrice di Householder di dimensione m 1. Notiamo che essa ancora ortogonale e vale:

                  (  (1)  (1)  (1)       (1) )
                    a11  a12  a13  ⋅⋅⋅  a1n
                  ||  0   a(222) a(223) ⋅⋅⋅  a(22)n ||
H2A (1) = H2H1A = ||  0    0   a(323) ⋅⋅⋅  a(32)n || ≡ A (2)
                  ||   ..    ..    ..        ..  ||
                  (   .    .   (.2)       .(2) )
                     0    0   am3  ⋅⋅⋅  amn

Anche a questo passo si dovr avere che a22(2)0 se la matrice A ha rango massimo.
Generalizzando al passo i-esimo si avr quindi:

                     (   (1)                          (1)  )
                     |  a11  ⋅⋅⋅ ⋅⋅⋅    ⋅⋅⋅   ⋅⋅⋅   a1n  |
                     ||   0   ...                     ...   ||
                     ||   ..   ..   (i)                 (i)  ||
A(i) ≡ HiHi−1⋅⋅⋅H1A = ||   .     . aii    ⋅⋅⋅   ⋅⋅⋅   ain  ||
                     ||   ...        0   a(i)     ⋅⋅⋅  a(i)   ||
                     ||   .        .    i+1.,i+1       i+1.,n ||
                     (   ..        ..      ..           ..   )
                         0   ⋅⋅⋅  0    a(im),i+1  ⋅⋅⋅   a(im)n

con ajj(j)0, j = 1,...,i.
Possiamo quindi definire la matrice di Householder al passo i + 1 come:

                                                    (         )
       (    |       )                 (   (i)    )      a(ii)+1,i+1
         -Ii|--------                 |  ai+1,.i+1 |   ||     0   ||
Hi+1 = (    |       ) ∈ ℝm ×m,  H (i+1)|(     ..   |) ≡ |(     ..   |)
            |H (i+1)                      a(i)              .
                                          m,i+1           0

ottenendo cos ancora una matrice ortogonale che soddisfa:

Hi+1A(i) = Hi+1Hi ⋅⋅⋅H1A ≡ A (i+1).

con ai+1,i+1(i+1)0. Procedento in questo modo dopo n passi si arriva alla matrice:

                  (   (1)       (1) )
                     a11  ⋅⋅⋅ a1n
                  ||   0   ...   ...  ||
                  ||   .   .    (n) ||
A (n) ≡ Hn ⋅⋅⋅H1A = ||  ..    .. ann  || ≡ R
                  ||   ...        0   ||
                  ||   .         .  ||
                  (   ..         ..  )
                      0   ⋅⋅⋅  0

con aii(i)0, i = 1,...,n. Se poniamo quindi QT = Hn⋅⋅⋅H1, abbiamo ottenuto la fattorizzazione A = QR.

Implementazione





Flops Memoria



LU O(2
3n3) n2



LDL O(13n3)  2
n2 + O(n)



PLU 2
3n3 n2 + O(n)



QR O(23n2(3m n)) O(m n)




Tabella 3.1: Tabella comparativa dei costi di fattorizzazione

Capitolo 4
Approssimazione di funzioni

Il problema che ci poniamo in questo capitolo quello di determinare una conveniente approssimazione di una generica funzione
f : [a,b] ⊂ ℝ → ℝ

Questo pu rendersi necessario per diversi motivi:

4.1 Interpolazione polinomiale

Avendo a disposizione i seguenti dati:

(xi,fi),  i = 0,...,n, fi ≡ f(xi),

ricerchamo un polinomio interpolante la f(x) sulle ascisse xi [a,b], xi < xj se i < j,xixj se ij, tale che:

p(x) ∈ Πn, p(xi) = fi,  i = 0,..,n.

Teorema
∃! pn ∈ Πn tale che pn(xi) = fi, i = 0,..,n.


Dimostrazione:
Un generico polinomio di grado n sar nella forma:

      ∑n    k
p(x) =    akx
      k=0

che descrive il sistema di equazioni lineari:

              (                 ) (    )    (    )
                x00  x10  ⋅⋅⋅  xn0     a0        f0
              || x01  x11  ⋅⋅⋅  xn1 || || a1 ||    || f1 ||
Va = f,  −→    |(  ...   ...       ...  |) |(  ... |)  = |(  ... |)
                x0  x1  ⋅⋅⋅  xn     a         f
                 n   n        n      n         n

in cui V una matrice di Vandermonde trasposta univocamente definita dalle ascisse xi. Il determinante di tale matrice vale:

det(V) = ∏ (x − x)
        i>j  i   j

In quanto per ipotesi le ascisse sono tra loro distinte, abbiamo che V nonsingolare e quindi esiste ed unica la soluzione del sistema lineare, ovvero esiste ed unico il polinomio suddetto. Non utilizzeremo per questa matrice per determinare in polinomio, in quanto tra le altre cose sappiamo anche che essa diviene molto mal condizionata al crescere di n.

4.2 Forma di Lagrange e forma di Newton

Per ottenere un problema discreto con propriet piu favorevoli rispetto a quello V a = f sopra descritto, si pu pensare di utilizzare basi diverse per Πn.
Ad esempio la base di Lagrange:

          n∏    x− xj
Lkn(x) =      x--−-x-, k = 0,1,...,n.
        j=0,j⁄=k k    j

I polinomi appena descritti sono tutti ben definiti in quanto, per ipotesi, xixj con ij.

Lemma

Dati i polinomi di Lagrange appena descritti definiti sulle ascisse a x0 < x1 < ... < xn b:

         {
           1, se k = i,
Lkn(xi) =   0, se k ⁄= i,

inoltre:

Teorema (forma di Lagrange)

Il polinomio

      ∑n
p(x) =   fkLkn(x)
      k=0

appartiene a Πn e soddisfa i vincoli di interpolazione p(xi) = fi, i = 0,1,...,n.
Dimostrazione:
p(x) Πn in quanto per definizione produttoria di n termini di grado 1; inoltre si ha che:

      ∑n             ∑n
p(x) =   fkLkn(xi) =      fkLkn(xi)+fiLin (xi) = fi.
      k=0           k=0,k⁄=i
                    ◟-----◝◜-----◞
                        =0perdef.

La base di Lagrange genera una forma del polinomio interpolante in cui i coefficienti sono facili da calcolare; questa rappresentazione per si presta male alla computazione in quanto non consente di generare il polinomio in maniera incrementale ma, dovendo calcolare e poi sommare k polinomi di grado n, ha un costo computazionale molto elevato.

Base di Newton

Indichiamo con pr(x) il polinomio in Πr che interseca la funzione da approssimare sulle ascisse x0,x1,...,xr che supponiamo noto, possiamo definire il successivo pr+1(x) a partire da questo, e questo pu essere fatto in modo semplice utilizzando la base di Newton cos definita:

ω0(x) ≡ 1
ωk+1 = (x− xk)ωk(x), k = 0,1,2,...

che gode delle seguenti propriet:

            ′
1)  ωk(x) ∈ Πk,∏kovvero `e un polinomio monico di grado k,
2)  ωk+1(x) =  j=0(x − xj),
3)  ωk+1(xj) = 0, per j ≤ k,
4)  ω0(x),...,ωk(x) costituiscono una base per Πk.

Teorema

La famiglia di polinomi interpolanti pr(x)r=0n tali che:

si genera ricorsivamente in questo modo:

p0(x) = f0ω0(x)
pr(x) = pr−1(x)+ f[x0,x1,...,xr]ωr(x), r = 1,..,n.
(4.1)

dove:

                r
f[x ,x ,...,x ] = ∑ ∏------fk-------.
   0  1    r   k=0  rj=0,j⁄=k(xk − xj)
(4.2)

e la quantit f[x0,x1,...,xr] detta differenza divisa di ordine r della funzione f(x) sulle ascisse x0,x1,...,xr.
Dimostrazione:
Ragionando per induzione, la tesi banalmente vera per r = 0. Supponendola quindi vera per r 1 vogliamo dimostrare che vera anche per r:

pr(x) ∈ Πr, infatti pr−1(xi) ∈ Πr−1 e ωr(x) ∈ Π′r.

Inoltre per i < r si ha:

pr(xi) = pr−1(xi)+ f[x0,x1,...,xr]ωr(xi) = pr−1(xi) = fi

in quanto ωr(xi) = 0 e per l’ipotesi di induzione. Se i = r:

p(x ) = p  (x )+ f[x,x ,...,x]ω (x ) = f
 r r    r−1  r      0 1     r r  r    r

e quindi

f[x0,x1,...,xr] = fr −-pr−-1(xr).
                  ωr(xr)

L’espressione ben definita in quanto, essendo le ascisse distinte, ωr(xr)0. Verifichiamo adesso che l’espressione appena ottenuta formalmente uguale a quella (4.2):
esprimendo il polinomio pr(x) secondo la base di Lagrange:

       ∑r
pr(x) =   fkLkr(x)
       k=0

Per l’unicit del polinomio interpolante questi polinomi devono conicidere con quelli in (4.1). I loro coefficienti principali devono quindi a loro volta conicidere.
Considerando la natura iterativa della forma di Newton si ottiene l’espressione:

              n
p(x) ≡ p (x) = ∑ f[x ,...,x ]ω (x).
       n     r=0   0    r  r

Osservazione:

Il denominatore di Lkn(x) dato da ωn+1(xk), mentre quello di f[x0,x1,...,xr] ωr+1(xk)
Le differenze divise di f(x) sulle ascisse a x0 < x1 < ... < xn b definisco i coefficienti della rappresentazione del polinomio interpolante p(xi) rispetto alla base di Newton.

Teorema

Le differenze divise di f(x) sulle ascisse a x0 < x1 < ... < xn b godono delle seguenti propriet:

1.
se α,β e f(x),g(x) sono funzioni di una variabile reale,
(α ⋅f +β ⋅g)[x0,...,xr])α⋅f[x0,...,xr]+ β ⋅g[x0,...,xr];

2.
per ogni i0,i1,...,ir permutazione di 0,1,...,r,
f[xi0,...,xir] = f [x0,...,xr];

3.
sia f(x) = k=0kaixi Πk, allora
              {
                 ak,  se r = k,
f[x0,x1,...,xr] =   0,   se r > k;

4.
se f(x) C(r), allora:
                (r)
f[x0,x1,...,xr] = f-(ε), ε ∈ [min xi,max xi];
                 r!          i     i

5.
                   f[x1,...,xr]− f[x0,...,xr−1]
f[x0,x1,...,xr−1,xr] =---------xr −-x0-------

Osservazione

Dal punto quattro consegue che se le ascisse coincidono (x0 = ... = xn), il polinomio pn(x) = r=0nf[x0,...,xr]ωr(x) diventa il polinomio di Taylor di f(x) di punto iniziale x0.
La propriet descritta nel punto cinque invece molto interessante dal punto di vista algoritmico: spiega infatti come sia possibile costruire una differenza divisa di ordine r a partire da due differenze divise di ordine r 1 che differiscono solo per una delle ascisse. Tenendo conto che f[xi] = fi, i = 0,1,...,n, siamo in grado di generare la seguente tabella triangolare, in cui la diagonale principale contiene i coefficienti del polinomio interpolante nella forma di Newton:

0 1 2 ⋅⋅⋅ n 1 n







x0 f[x0]
x1 f[x1] f[x0,x1]
.
.. .
.. f[x1,x2] f[x0,x1,x2]
..
. ..
. ..
. ..
. ...
xn1f[xn1]f[xn2,xn1] ..
. f[x0,...,xn1]
xn f[xn] f[xn1,xn] f[xn2,xn1,xn]⋅⋅⋅ f[x1,...,xn] f[x0,...,xn]

Notiamo che le colonne della tabella vanno calcolate da sinistra verso destra; possiamo per calcolare ciascuna colonna dal basso verso l’alto: in questo modo gli elementi adiacenti a sinistra non sono pi necessari nel seguito e si posso quindi riscrivere con il nuovo elemento appena calcolato.

Implementazione

Questo codice implementa il calcolo delle differenze divise appena descritto:


PIC

Figura 4.1: Questo grafico mostra l’andamento del polinomio di Newton (in rosso) interpolante la funzione xsin(3x) + x sulle ascisse equidistanti [0 : 0.5 : 3].


4.3 Interpolazione di Hermite

Consideriamo le seguenti ascisse di interpolazione:

a ≤ x0 < x12 < ...< xn < xn+12 ≤ b = 2n+ 2 ascisse,

                                   1        1
p(x) ∈ Π2n+1 tali che p(xi) = fi, i = 0, 2,1,...,n+ 2 .

Se facciamo in modo che:

xi+ 1→  xi,  i = 0,1,...,m, e  p(xi) = f (xi),p(xi+1) = f(xi+ 1),
   2                                       2         2

si ha che:

p(xi+ 1)− p(xi) = f(xi+1)− f(xi)
     2              2

e dividendo membro a membro per xi+1
2 xi si ottiene:

                        ′      ′
p[xi,xi+ 12] = f[xi,xi+12] ⇒ p(xi) = f(xi)

Come conseguenza abbiamo qundi che il polinomio p(x) ancora definito e soddisfa i vincoli di interpolazione:

p(xi) = f(xi), p′(xi) = f ′(xi), i = 0,1,...,n.

Abbiamo appena definito il polinomio interpolante di Hermite.
Passando al polinomio di Taylor, considerando di avere una sola ascissa:

                                      ′
p(x)) = f[x0]+ f[x0,x0](x− x0) ≡ f(x0) +f (x0)(x − x0)

Per generalizzare, se si considerano le ascisse in questo modo:

a ≤ x = x  < x = x  < ...< x = x  ≤ b,
     0   0    1   1        n   n

il polinomio di Hermite sar dato da:

p(x)  =  f [x0]+ f[x0,x0](x− x0)+ f[x0,x0,x1](x − x0)2+
        f [x0,x0,x1,x1](x− x0)2(x − x1)+
        f [x0,x0,x1,x1,x2](x − x0)2(x − x1)2 + ⋅⋅⋅+
        f [x0,x0,x1,x1,...,xm,xm ](x− x0)2⋅⋅⋅(x − xm−1)2(x− xm ).

Per calcolare i coefficienti ci serviremo dell’algoritmo implementato qui di seguito, che non altro che una modifica all’algoritmo per il calcolo delle differenze divise. In questo algoritmo avr in ingresso un vettore f contenente i valori:

f(x ),f′(x ),f(x ),f′(x ),...,f(x ),f′(x  ),
  0     0     1     1       m     m

ed un vettore x contenente le ascisse.

Implementazione

Questo codice implementa il calcolo delle differenze divise appena descritto:


PIC

Figura 4.2: Questo grafico mostra l’andamento del polinomio di Hermite (in rosso) interpolante la funzione xsin(3x)+x sulle ascisse {0,1.5,2,3}; notiamo come gi con 4 punti invece dei 7 usati per il polinomio di Newton si abbia un’approssimazione della funzione molto buona.


4.4 Errore nell’interpolazione

Sia p(x) Πn il polinomio interpolante tale che, note le coppie (xi,fi), p(xi) = fi con i = 0,1,...,n. Possiamo definire l’errore di interpolazione come:

e(x) = f (x) − p(x)

dove e(x) l’errore che si commette nell’approssimare f(x) mediante il suo polinomio interpolante p(x). Segue naturalmente che, per definizione, e(xi) = 0, i = 0,1,...,n, ma di nostro interesse misurare questa funzione quando x non appartiene a x0,...,xn.

Teorema

Sia

       n
p(x) = ∑  f[x ,...,x ]ω (x),
      r=0   0     r r

tale da soddisfare i vincoli di interpolazione p(xi) = fi con i = 0,1,...,n. Il suo errore di interpolazione :

e(x) = f [x0,x1,...,xn,x]ωn+1(x).


Dimostrazione:
Si consideri un generico punto x per semplicit distinto dalle ascisse di interpolazione x0,...,xn, ed il polinomio p(x) Πn+1 che soddisfi sia gli stessi vincoli di interpolazione di p(x) che:

p(x) = f (x)

Dal teorema sulla forma di Newton segue:

p(x) = p(x)+ f[x ,x ,...,x ,x]ω  (x).
               0  1    m    n+1

Tenendo conto del vincolo aggiunto e della genericit del punto x:

f(x)− p(x) = f[x ,x ,...,x ,x]ω  (x),
               0  1    m    n+1

e quindi abbiamo appena ottenuto la tesi. Cvd.

Corollario

Se la funzione f(x) C(n+1), allora:

       f(n+1)(ξx)
e(x) = -(n-+-1)! ωn+1(x)

con ξ [min{x0,x},max{xn,x}].
Dalla struttura dell’errore di interpolazione si capisce che questo fondamentalmente guidato da due componenti:

Si ha pertanto che ωn+1(x) oscilla quando x [x0,xn] e si annulla nelle ascisse di interpolazione; ma |ωn+1(x)| cresce quanto |xn+1| se x < x0 o x > xn. Si pu concludere quindi che generalmente p(x) pu essere utilizzato in maniera utile per approssimare la funzione di partenza f(x) solo all’interno del range [x0,xn].


PIC

Figura 4.3: Questo grafico mostra l’andamento del polinomio di Newton interpolante la funzione 1(1 + x2) al variare del numero delle ascisse; abbiamo in nero la funzione originaria, in ciano verde e rosso rispettivamente i polinomi interpolanti in: [5 : 2 : 5],[5 : 1 : 5],[5 : 0.5 : 5]. Notiamo come l’errore aumenti al crescere del numero dei punti di interpolazione.

PIC

Figura 4.4: Questo grafico mostra l’andamento dell’errore rispetto al numero di punti di interpolazione scelti in modo equidistante nell’intervallo [5,5] per il polinomio di Newton interpolante la funzione 1(1 + x2). Si nota facilmente che l’errore diverge al crescere del numero dei punti di interpolazione.



PIC

Figura 4.5: Questo grafico mostra l’andamento del polinomio di Newton interpolante la funzione |x| al variare del numero delle ascisse; abbiamo in nero la funzione originaria, in ciano verde e rosso rispettivamente i polinomi interpolanti in 5, 10 e 15 punti equidistanti. Come nel caso della funzione di Runge, l’errore aumenta al crescere del numero dei punti di interpolazione.

PIC

Figura 4.6: Questo grafico mostra l’andamento dell’errore rispetto al numero di punti di interpolazione scelti in modo equidistante nell’intervallo [5,5] per il polinomio di Newton interpolante la funzione di Bernstein. Come nel caso della funzione di Runge, l’errore diverge in maniera molto pronunciata al crescere del numero di punti di interpolazione.


4.5 Condizionamento del problema

Per effettuare lo studio del condizionamento del problema occorre per prima cosa identificare quali siano i dati in ingresso per i quali una perturbazione di questi si rifletta sul risultato finale. Per definire il polinomio p(x) ci siamo serviti di un insieme di coppie di numeri (fi,xi), delle quali, in generale, possiamo dire che:

Possiamo quindi considerare l’insieme dei valori delle ascisse come parametri fissati del problema, e quindi considerare le ordinate gli unici dati in ingresso soggetti ad errore. Consideriamo la funzione f˜ (x) una perturbazione di f(x), e indicheremo con ˜f i le valutazioni di ˜f (xi) come gi fatto in precedenza con la f(x). Definiamo adesso, espressi nella forma di Lagrange, i due polinomi:

      ∑n                  ∑n
p(x) =   fkLkn(x),  ˜p(x) =    ˜fkLkn(x)
      k=0                 k=0

che sono rispettivamente il polinomio costruito a partire dai dati esatti e quello invece costruito sui dati perturbati. Confrontandoli si ottiene:

             |                        |  |               |
             ||∑n           ∑n         ||  ||∑n              ||
|p(x )− ˜p(x)| = ||  fkLkn(x)−    ˜fkLkn(x)|| = || (fk − f˜k)Lkn(x)||
              k=0          k=0           k=0

e per disugualianza triangolare possiamo affermare:

|               |                       (           )
||∑n      ˜       ||   n∑       ˜             ∑n                   ˜
||  (fk − fk)Lkn(x)|| ≤   |fk − fk|⋅|Lkn(x)| ≤    |Lkn(x)|  makx |fk − fk|
k=0                 k=0                    k=0

                ˜
≡ λn(x)maxk |fk − fk|

dove λn(x) detta funzione di Lebesgue. Per come stata definita, λn(x) dipende solo dalla scelta delle ascisse di interpolazione. Definiamo adesso:

∥f∥ = maa≤xx≤b|f (x)|.

Abbiamo quindi che:

∥p− ˜p∥ ≤ ∥λn∥⋅∥f −f˜∥ ≡ Λn∥f − ˜f∥,

dove Λn detta costante di Lebesgue, e dipende soltanto dalla scelta dell’intervallo [a,b] considerato e dalla scelta delle ascisse di interpolazione all’interno di questo.
Abbiamo ottenuto che f ˜
f una misura dell’errore sui dati in ingresso, mentre p ˜p misura l’errore sul risultato. La costante di Lebesgue appena definita misura pertanto l’amplificazione massima sul risultato dell’errore dei dati in ingresso, e definisce quindi il numero di condizionamento del problema. utile ricordare che Λn O(log n) per n →∞, ovvero il problema diviene sempre pi malcondizionato al crescere di n. Inoltre la scelta di ascisse equidistanti genera una successione {Λn} che diverge in modo esponenziale al crescere di n; ne consegue che questo metodo non una scelta molto arguta se si vogliono usare valori molto elevati per n.
Vediamo adesso come l’errore di interpolazione ed il condizionamento del problema siano strettamente interconnessi:

Teorema

Sia f(x) una funzione continua in [a,b], allora esiste p(x) Πn tale che:

     ∗
∥f − p ∥ = mp∈iΠnn ∥f − p∥.

Questo polinomio detto “polinomio di miglior approssimazione di grado n di f(x) sull’intervallo [a,b]”.

Teorema

Sia p(x) il polinomio di migliore approssimazione di grado n di f(x). Allora si ha che:

∥e∥ ≤ (1 +Λn )∥f − p∗∥.


Dimostrazione:
considerando che p(x) Πn coincide con il proprio polinomio interpolante sulle ascisse a x0 < ... < xn b, si ha:

                   ∗    ∗            ∗     ∗                   ∗
∥e∥ = ∥f − p∥ = ∥f − p + p − p∥ ≤ ∥f − p∥ + ∥p − p∥ ≤ (1+ Λn )∥f − p ∥

come volevasi dimostrare.
Come conseguenza, non detto che l’errore diminuisca al crescere di n se la costante di Lebesgue cresce esponenzialmente. Per precisare, introduciamo il “modulo di continuit di una funzione” relativo all’intervallo [a,b]:

ω(f;h) ≡ sup{|f(x)− f(y)| : x,y ∈ [a,b],|x− y| ≤ h}

dove h > 0 un parametro assegnato. Possiamo osservare che se f(x) C(1), allora ω(f;h) 0 se anche h 0,e che se f(x) Lipschitziana con costante L, allora ω(f;h) Lh.

Teorema di Jackson

Il polinomio di approssimazione di una funzione p(x) = minpΠnf p,f(x) C(0), tale che:

             (       )
∥f − p∗∥ ≤ α ⋅ω f; b−-a
                   n

in cui la costante α indipendente da n. Possiamo quindi concluedere che, per una funzione generica vale:

                (       )
∥e∥ ≤ α(1+ Λ )ω  f; b-− a .
            n        n

4.6 Ascisse di Chebyshev

Trattiamo adesso la risoluzione del seguente problema:

   min     ∥ωn + 1∥ ≡    min       max |ωn+1(x)|
a≤x0<...<xn≤b          a≤x0<...<xn≤b  a≤x≤b

che deriva dallo studio del condizionamento del problema, in cui ci si accorge che una scelta oculata delle ascisse di interpolazione deve essere fatta tenedo conto che:

Per semplicit, assumiamo che:

[a,b] ≡ [− 1,1]

In questo modo non si perde in generalit, infatti sempre possibile costruire una corrispondenza biunivoca tra gli insiemi Z = {z : z [1,1]} e Y = {y : y [a,b]}. Si ha che:

   ( a+ b   b− a )
y ≡  --2--+ -2--z  ∈ [a,b]

infatti se z = 1, y = a e se z = 1, y = b; per tutti i valori di z compresi tra 1 ed 1 si ha che alla quantit a+b
 2 (che indica la met esatta dell’intervallo [a,b]) viene sommata una quantit a+b
 2 b−a
 2z a+b
 2 e quindi tutti i valori risultanti devono appartenere all’intervallo [a,b].
Viceversa, si ha che:

    2y−-a-−-b
z ≡   b− a   ∈ [− 1,1]

infatti se y = a, z = 1 e se y = b, z = 1; per tutti i valori di y compresi tra a e b si ha che |2y ab|≤|ba| e quindi tutti i valori risultanti appartengono all’intervallo [1,1].
Possiamo adesso definirela famiglia di polinomi di Chebyshev di prima specie come:

T0(x) ≡ 1
T1(x) = x
Tk+1(x) = 2xTk(x)− Tk−1(x), k = 1,2,...

Valgono le seguenti propriet:

1.
Tk(x) un polinomio di grado esatto k
2.
Il coefficiente principale di Tk(x) 2k1, k = 1,2,...
3.
La famiglia di polinomi ˆTk,, in cui
 ˆ             ˆ       1− k
T0(x) = T0(x), Tk(x) = 2   Tk(x),  k = 1,2,...

una famiglia di polinomi monici di grado k, k = 1,2,....

4.
Se si pone x = cosθ, con θ [0], per parametrizzare i punti dell’intervallo [-1,1] rispetto a θ si ottiene:
Tk(x) ≡ Tk(cosθ) = cos(kθ), k = 0,1,....

Teorema

Gli zeri di Tk(x), tra loro distinti, sono dati da:

         (        )
x(ik)= cos (2i+-1)π  ,  i = 0,1,...,k− 1
             2k

Per x [1,1], i valori estremi di Tk(x) sono assunti nei punti:

        ( i )
ξ(ik)= cos  -π  ,  i = 0,1,...,k
          k

e corrispondono a:

Tk(ξ(k)) = (− 1)i,  i = 0,1,...,k.
   i

Si ha quindi che Tk= 1, e ˆTk= 21kTk= 1. Inoltre per k = 1,2,...:

21− k = ∥ˆTk∥ = min ∥φ∥.
             φ∈Π′k

Segue quindi che, sulle ascisse di interpolazione a x0 < ...axn b per l’intervallo [1,1,] scelte in questo modo:

          ((2i+ 1)π)
xn− i = cos--------  ,  i = 0,1,...,n
            2(n + 1)

si ottiene:

          n∏
ωn+1 (x) =   (x− xi) ≡ Tˆn+1(x)
          i=0

che la soluzione del problema che ci eravamo posti. Con questa scelta delle ascisse, se la funzione sufficientemente regolare si ottiene che:

      ∥f(n+1)∥
∥e∥ ≤ (n-+-1)!2n-

Si ha quindi che la corrispondente costante di Lebesgue vale:

Λ  ≈ 2-logn
 n   π

che risulta quindi avere una cresctia ottimale anche per un numero di ascisse tendente all’infinito.


PIC

Figura 4.7: Questo grafico mostra l’andamento del polinomio di Newton interpolante la funzione 1(1 + x2) al variare del numero delle ascisse; abbiamo in nero la funzione originaria, in ciano verde e rosso rispettivamente i polinomi interpolanti nelle ascisse di Chebyshev colcolate su 5,11 e 21 punti nell’intervallo [5,5]. Notiamo come, al contrario di come succedeva con le ascisse equidistanti, l’errore si mantenga molto pi basso.

PIC

Figura 4.8: Questo grafico mostra l’andamento dell’errore rispetto al numero di punti di interpolazione scelti con il criterio di Chebyshev per il polinomio di Newton interpolante la funzione 1(1 + x2). Si nota facilmente che l’errore diminuisce all’aumentare del numero di punti di interpolazione, contrariamente a quanto accadeva con le ascisse equidistanti.


4.7 Funzioni spline

Abbiamo appena visto come al crescere del grado del polinomio interpolante sia necessario effettuare una scelta delle ascisse di interpolazione che non faccia crescere in maniera importante la costante di Lebesgue; d’altro canto se si mantiene basso il numero n delle ascisse di interpolazione il modulo di continuit di f, ω = (f; b−-a)
    n non pu tenedere a zero una volta fissato l’intervallo [a,b]. Possiamo per fissare una partizione dell’intervallo originario Δ = {a = x0 < x1 < ... < xn = b}, per la quale h = maxi=1,...,n(xi xi1) tende a zero se n tende all’infinito, e considerare, per ogni sottointervallo [xi1,xi] della partizione, un polinomio di grado m fissato interpolante f negli estremi del sottointervallo. In questo modo si ha che il problema del condizionamento perde importanza in quanto m rimane costante (e quindi anche l’errore), mentre a questo punto vale che ω(f;h) tende a zero per n che tende all’infinito. Quella che abbiamo appena descritto una funzione polinomiale a tratti, definiamola pi esattamente:

Definizione

sm(x) una spline di grado m sulla partizione Δ se:

Se si aggiunge la condizione:

sm(xi) = fi, i = 0,1,...,n

allora la spline quella interpolante la funzione f(x) nei nodi della partizione Δ.

Teorema

Se sm(x) una spline di grado m sulla partizione Δ = {a = x0 < x1 < ... < xn = b}, allora sm(x) una spline di grado m 1 sulla stessa partizione.
Dimostrazione:
Se sm(x) (m1), allora sm(x) (m2) e se sm|[xi1,xi](x) pi(x) Πm, allora p(x) Πm1.

Teorema

L’insieme delle funzioni spline di grado m definite sulla partizione Δ uno spazio vettoriale di dimensione m + n.

Ne consegue quindi che sono necessarie m + n condizioni indipendenti per individuare univocamente la spline interpolante una funzione su una partizione assegnata (m + n1 se si considera n l’esatto numero di punti di interpolazione); in quanto tipicamente abbiamo in totale n + 1 condizioni di interpolazione (le coppie (xi,fi)), queste permettono l’individuazione unica della spline lineare, che coincide con la spezzata congiungente i punti {(xi,fi)}i=0,...,n in questo modo:

            (x − x   )f + (x − x)f
s1|[xi−1](x) = -----i−-1-i----i-----i−1-, i = 1,...,n.
                   xi − xi−1

Per identificare una spline di ordine pi elevato necessaria l’imposizione di oppurtune condizioni aggiuntive.

4.7.1 Spline cubiche

Come detto prima, per ottenere una spline cubica occorre imporre, oltre alle n + 1 condizioni di interpolazione, due condizioni aggiuntive, dalla cui scelta si ottengono spline cubiche differenti.

Spline naturale

Questa scelta consiste nell’imporre:

 ′′          ′′
s3(a) = 0, s3(b) = 0

Spline completa

In questo caso supponiamo di conoscere i valori di f(x) negli estremi della partizione, e quindi imponiamo:

s′(a) = f′(x), s′(b) = f ′(b)
 3             3

Spline periodica

Se la funzione che stiamo interpolando periodica e l’intervallo [a,b] contiene un numero finito e intero di periodi, si possono imporre le condizioni:

s′(a) = s′(b), s′′(a) = s′′(b).
 3     3      3      3

Condizioni not-a-knot

In questo caso imponiamo che uno stesso polinomio in Π3 costituisca la restrizione della spline sull’intervallo [x0,x1] [x1,x2], ed un’altro polinomio analogo faccia lo stesso su [xn1,xn1] [xn1,xn]. Le condizioni da imporre saranno quindi:

s′3′′|[x0,x1](x1) = s′′3′ |[x1,x2](x1), s′3′′|[xn−1,xn−1](xn −1) = s′′3′ |[xn−1,xn](xn−1)

Se si osserva che s3′′′|[xi1,xi](x) Π0, si pu esprimere le condizioni suddette come:

s′′3(x1)− s′3′(x0)   s′′3(x2)− s′3′(x1)
---x1 −-x0---=  ---x2 −-x1---

e

s′3′(xn−1)−-s′′3(xn−2)= s′3′(xn)−-s′3′(xn−1).
    xn−1 − xn−2         xn − xn−1

4.7.2 Calcolo di una spline cubica

Andiamo adesso ad approssimare una funzione f(x) definita su un’intervallo [a,b] con una spline cubica, prendendo in esame una partizione dell’intervallo Δ = {a = x0 < ... < xn = b}. Osserviamo come se s3(x) una spline cubica, allora s3(x) una spline quadratica e s3′′(x) una lineare. Vale quindi che:

 ′′           (x−-xi−1)mi +-(xi −-x)mi−1
s3|[xi−1,xi](x) =           hi           ,  x ∈ [xi−1,xi]

dove:

      ′′
mi ≡ s3(xi),  i = 0,1,...,n,
hi = xi − xi− 1, i = 1,2,...,n.

Integrando questa equazione si ottiene :

             (x − xi−1)2mi +(xi − x)2mi−1
s′3|[xi−1,xi](x) =------------2h------------ + qi,  x ∈ [xi−1,xi]
                           i

dove qi una costante di integrazione, e integrando ancora una volta si ottiene:

             (x-−-xi−1)3mi-+(xi −-x)3mi−1
s3|[xi−1,xi](x) =            6hi            + qi(x− xi−1)+ ri,  x ∈ [xi− 1,xi]

dove ri una costante di integrazione. Imponiamo adesso i vincoli di interpolazione, ovvero che s3(xi) = fi, i = 0,1,...,n, per ricavare ri e qi:

s (x) = h2im  + qh + r = f
 3  i   6  2i   ii   i   i
s3(xi− 1) = himi−1 + ri = fi−1
          6

da cui:

r = f   − h2im
 i   i−1   6   i−1
qi = fi −-fi−-1− hi(mi − mi−1)
       hi      6

Abbiamo quindi che :

                      2           2
s′3|[xi−1,xi](x) = (x−-xi−1)-mi +-(xi −-x)-mi−1 + fi −-fi−1-− hi(mi − mi −1),
                         2hi                 hi     6

sempre con x [xi1,xi]. Imponiamo adesso la continuit di s3(x), tenendo conto che:

quindi:

          s′3|[xi−1,xi](xi)            =              s′3|[xi,xi+1](xi)
 hi     hi                             hi+1-               hi+1
 2 mi − 6 (mi − mi−1)+ f[xi−1,xi]   =  −  2  mi + f[xi,xi+1]−  6 (mi+1 − mi )
himi −1 + mi(2hi + 2hi+1)+ hi+1mi+1 =         6(f[xi,xi+1]− f[xi−1,xi])

e dividendo membro a membro per hi + hi+1(= xi+1 xi1):

        hi           hi+1
φi = hi +-hi+1, ξi = hi-+-hi+1, i = 1,...,n− 1.

e quindi:

φ m   + 2m  + ξmi + 1 = 6f [x   ,x,x   ],  i = 1,...,n− 1.
 i i−1     i   i           i−1  i i+1

Si osserva che φii > 0 e φi + ξi = 1 per i = 1,...,n 1.
Poniamoci adesso nel caso delle spline naturali, e imponiamo quidni la condizione m0 = mn = 0. Si ottiene immediatamente il sistema:

(                         ) (       )    (                 )
   2   ξ1                       m1             f[x0,x1,x2]
|| φ2   2   ξ2             || ||  m2   ||    ||    f[x1,x2,x3]   ||
||     ...  ...   ...        || ||   ...   || = 6||         ...       ||
||          .    .         || ||   .   ||    ||         .       ||
(          ..   ..   ξn−2 ) (   ..   )    (         ..       )
               φn−1   2       mn −1         f[xn−2,xn− 1,xn]

Notiamo come questa matrice goda di diverse propriet:

Condizioni not-a-knot

Andiamo adesso ad imporre le condizioni not-a-knot:

s′3′(x1)−-s′′3(x0)  s′3′(x2)−-s′′3(x1)   s′3′(xn)−-s′3′(xn−1)-  s′′3(xn−-1)-− s′3′(xn−2)
    x1 − x0   =     x2 − x1   ;     xn − xn− 1  =     xn−1 − xn −2

che significa porre:

m2-−-m1-= m1-−-m0-,  mn-−-mn-−1=  nn−1 −-nn−2
   h2        h1          hn          hn−1

da cui si ricava:

                 m2h1 − m1h2  =   h2m1 − h2m0
-m0h2--− m1(h1 +-h2)-+-h1m2-- =   0
h1 + h2 h  h1 + h2   hh1 + h2
   m0 ---2--− m1 + ---1---m2  =   0
     ◟h1 +◝◜h2◞      h◟1 +◝◜h2◞
        ξ1           φ1

in modo analogo si ricava che:

 ---hn----mn−2 − mn−1 + --hn−1---= 0
◟hn +◝h◜n−1◞              h◟n +◝h◜n−-1◞
   ξn−1                   φn−1

che da origine alla matrice:

(                          ) (       )     (        0       )
   ξ1  − 1  φ1                   m0        |   f [x0,x1,x2]  |
||  φ1   2    ξ1             || ||   m1  ||     ||   f [x1,x2,x3]  ||
||      ...  ...   ...       || ||    ...  || = 6 ||        .       ||
|(          φ      2   ξ    |) |(  m    |)     ||        ..       ||
            ξn− 1 − 1  φn−1       nm−1       ( f[xn−2,xn−1,xn])
            n−1        n−1        n                 0

Sommando la seconda riga alla prima e la penultima all’ultima si ottiene:

                                               (                )
(  1   1     1            ) (   m0 + m1  )     |    f[x0,x1,x2]  |
||  φ1  2    ξ1            || ||     m1     ||     ||    f[x0,x1,x2]  ||
||      ..   ..   ..       || ||      ..     || = 6 ||    f[x1,x2,x3]  ||
|       .     .    .      | |      .     |     ||        ...       ||
(          φn−1   2  ξn−1 ) (    mn −1   )     ( f [xn−2,xn−1,xn])
             1    1    1       mn + mn−1         f [xn−2,xn−1,xn]

che ancora non una matrice tridiagonale a causa degli elementi sottolineati. Possiamo per sottrarre la prima colonna dalla seconda e dalla terza e l’ultima dalla penultima e terzultima per ottenere finalmente una matrice tridiagonale. Per fare questo ci serviremo di una matrice siffatta:

    (                            )         (                      )
      1  − 1  − 1                             1  1  1
    || 0   1   0                  ||         ||  0  1  0             ||
    || 0   0   1                  ||         ||  0  0  1             ||
T = ||             ...             || ; T −1 = ||          ...         || ;
    ||                  1   0   0 ||         ||               1 0  0 ||
    |(                  0   1   0 |)         |(               0 1  0 |)
                      − 1 − 1  1                           1 1  1

Quello che andremo a svolgere invece di Am = b diviene ATT1m = b, e la matrice AT quindi:

(  1     0                                            )
|  φ1  2− φ1  ξ1 − φ1                                 |
||        φ2      2      ξ2                             ||
||               .      .         .                    ||
||                ..     ..       ..                   ||
||                     φn− 2      2        ξn−2        ||
(                           φn −1 − ξn−1 2− ξn−1  ξn− 1)
                                            0      1

(                   )    (                 )
    m0 + m1 + m2              f[x0,x1,x2]
||        m1         ||    ||    f[x0,x1,x2]   ||
||         ...         || = 6||         ...       ||
|(       m           |)    |(  f[x   ,x   ,x  ]|)
  m  + m  n−1+ m             f[xn−2,xn−1,x n]
    n    n−1    n−2            n−2  n−1  n

che tridiagonale e fattorizzabile LU in quanto ha tutti i minori principali non nulli.

Spline completa

Per avere una spline completa invece dobbiamo imporre:

s′3(a) = f′(a), s′3(b) = f′(b)

ovvero (considerando per comodit gli indici che partono da 1 anzich da 0 ed arrivano ad n, in accordo con quanto utilizzato nelle implementazioni):

s′(a) ≡ s′(x1) = (x1-−-x1)2m2-−-(x2 −-x1)2m1-+ f2 −-f1+ (x2 −-x1)(m2− m1) = f′(a);
 3     3              2(x2 − x1)         x2 − x1    6

 ′     ′      (xn-−-xn−1)2mn-−-(xn-−-xn)2mn-−1  fn −-fn−1 (xn-−-xn−1)              ′
s3(b) ≡ s3(xn) =         2(xn − xn−1)         + xn − xn− 1+    6     (mn − mn−1) = f (b);

da cui si ottiene:

                       ′                                          ′
− 3h1m1+h1 (m2− m1) = 6f(a)− 6f[x1,x2]; 3hn−1mn+hn −1(mn − mn −1) = 6f (b)− 6f [xn,xn−1];

ovvero:

− 4h1m1+h1m2 = 6(f′(a)− f[x1,x2]);  4hn−1mn − hn−1mn −1 = 6(f′(b)− f[xn−1,xn]).

Si ricava quindi la matrice:

                                                     (                  )
(  − 4h1 h1                            ) (  m1 )         f′(a)− f[x1,x2]
||   φ1    2   ξ1                       || ||   .. ||     ||      f[x1,x2]     ||
||        φ2   2    ξ2                  || ||   . ||     ||         ...        ||
||             ..    ..     ..           || ||   ... ||  = 6||         .        ||
|              .     .      .          | |   .. |     ||         ..        ||
(                 φn− 1    2     ξn−1  ) (   . )     (      f[x1,x2]     )
                         − hn−1  4hn −1     mn           f′(b)− f[xn−1,xn]

che ancora una volta tridiagonale con tutti i minori principali non nulli, e pertanto fattorizzabile LU.

Spline periodica

Per calcolare la spline cubica periodica invece dobbiamo imporre:

s′(a) = s′(b); s′′(a) = s′′(b).
 3     3      3      3

Sostituendo nelle formule (come nel caso precedente per semplicit sono stati usati gli indici Matlab) andiamo a ricavare:

Periodicit della derivata prima (s3(a) = s3(b)):
 (x2 − x1)m1           h1           (xn − xn−1)mn            hn− 1
−-----2-----+f [x1,x2]− -6 (m2 − m1 ) =----2------+f [xn,xn−1]−--6--(mn − mn−1)

− 3h1m1 − 3hn−1mn − h1(m2 − m1 )+ hn−1(mn − mn −1) = 6(f[xn,xn− 1]− f[x1,x2])

da cui otteniamo quindi la prima condizione:

− 2h1m1 − 2hn−1mn − h1m2 − hn−1mn −1 = 6(f[xn,xn −1]− f [x1,x2]).

Periodicit della derivata seconda (s3′′(a) = s3′′(b)):
(x1-− x1)m2-+-(x2 −-x1)m1 = (xn −-xn−-1)mn-+-(xn-− xn-)mn-−1
           h1                         hn−1

h m    h   m
-1-1-= -n−1--n   da cui la condizione m1 = mn
 h1      hn−1

Andiamo quindi a costruire la matrice relativa al sistema lineare risultante:

(                                                   )
     1     0    0    0    0     ...    0       -1
||   φ1     2    ξ1   0    0     ...    ...      0    || (  m  )
||    0    φ2    2   ξ2    0     ...    ...      0    || |   .1|
||    ...          ...  ...   ...                    ...   || ||   .. ||
||    .              .     .     .               .   || ||   .. ||
||    ..               ..   ..     ..             ..   || ||   .. ||
||    ..                                              || (   .. )
|(    .                  φn −2    2    ξn−2     0    |)    mn
     0                         φn−1    2      ξn−1
   − 2h1  − h1  0   ...  ...     0    -hn−1   − 2hn−1

   (                       )
               0
   ||        f [x1,x2]       ||
   ||            ...          ||
   ||            .          ||
= 6||            ..          ||
   ||            ...          ||
   |(       f[x ,x   ]      |)
      f[x ,x   n]−n f−1[x − x ]
         n  n−1     1   2

4.8 Approssimazione polinomiale ai minimi quadrati

Il problema che ci poniamo in questa sezione quello di determinare il polinomio che meglio approssima una serie di dati sperimentali, i quali tipicamente sono in sovrabbondanza ma soggetti ad errori. Vogliamo quindi determinare un polinomio di grado m:

   ∑m     k
y =    akx
   k=0

che meglio approssimi i dati:

(xi,yi),  i = 0,1,...,n, n ≥ m

dove le ascisse non sono necessariamente tutte distinte, anche se dobbiamo imporre che nella determinazione di un polnomio in Πm, almeno m + 1 lo siano. Definiamo quindi i vettori f e y rispettivamente dei valori osservati e previsti:

   (  f0 )              ( y0 )
   |   .. |    n+1       |  .. |     n+1               m∑     k
f = (  . ) ∈ ℝ   ;  y = (  . ) ∈ ℝ    , yi = pm (xi) =   akxi;
      fn                  yn                         k=0

rimane da determinare il vettore a = (a0,a1,...,am)T che minimizza la quantit:

      2   ∑n       2
∥f − y∥2 =    |fi − yi|
          i=0

ovvero:

                  (     )
                     a0
pm (xi) = (x0ix1i...xmi )|(  ... |)
                     am

Possiamo esprimere quindi il problema in questa forma:

             (   0   1       m )
                x0  x0  ⋅⋅⋅  x0
y = V a; V = |(  ...            ... |)  ∈ ℝn+1 ×m+1
               x0  x1   ⋅⋅⋅  xm
                 n   n       n

Si ha che:

pmi∈nΠ ∥f− y∥22 =  mimn+1 ∥Va − f∥22
m   m          a∈ℝ

Teorema

La matrice V ha rango m + 1.
Dimostrazione:
Se ci si limita a considerare m + 1 righe corrispondenti ad ascisse distinte, esse formano una matrice di Vandermonde nonsingolare, e quindi il rango appunto m + 1.
Il problema diventa quindi quello di risolvere, nel senso dei minimi quadrati, il sistema sovradeterminato:

V a = y

Possiamo quindi considerare:

           (    )
              ˆR      ˆ    m+1×m+1     T
V = QR = Q    0   , R ∈ ℝ         , Q  Q = Ir+1

e quindi:

mpi∈nΠ ∥f− y∥22 =  mimn+1 ∥Va− f∥22 =  mimn+1∥QRa  − f∥22 = ∥Q(Ra − Q◟T◝◜f◞)∥22
   m          a∈ℝ              a∈ℝ                         g

           ∥∥( ˆ)    (  ) ∥∥2  ∥∥( ˆ     ) ∥∥
= ∥Ra− g∥ = ∥∥ R  a−  g1  ∥∥ = ∥∥  Ra− g1  ∥∥= ∥Rˆa − g1∥22+ ∥g2∥22 ⇒ ˆRa = g1
           ∥  0      g2  ∥2  ∥   − g2   ∥

Ne deriva quindi che la soluzione di questa equazione, a = Rˆ1g1 esiste ed unica, ed quindi tale anche il polinomio di approssimazione ai mini quadrati. Notiamo che nel caso particolare in cui m = n, il polinomio di approssimazione ai minini quadrati coincide esattamente con quello interpolante sulle stesse ascisse; in tal caso infatti il vettore residuo g2 il vettore vuoto e il polinomio quindi interpola tutti i valori assegnati.

Osservazione

Pu succedere di aver bisogno di minimizzare, invece della norma del vettore dei residui, un vettore pesato D(V a f22 che permetta appunto di dare un’importanza maggiore o minore a certe misurazione rispetto alle altre. Definiamo quindi:

    (             )
    | w0          |
D = (      ...     ) ,  wi > 0,i = 0,...,m.
               wn

abbiamo quindi che:

          (             )
             p0(x0) − f0
r = y − f = |(    ..      |) .
                 .
            pm (xn)− fn

Quindi il problema:

      ∑n
∥r∥22 =   |ri|2
      i=0

diventa:

        ∑n        ∑n
∥Dr ∥22 =   |ri|w2i =   w2i(pm (xi)− fi)2.
        i=0        i=0


PIC

Figura 4.9: Questo grafico mostra l’andamento della spline cubica naturale (verde) e di quella not-a-knot (rossa) interpolanti la funzione 1(1 + x2) in sette punti equidistanti nell’intervallo [5,5]

PIC

Figura 4.10: Questo grafico mostra invece l’andamento delle stesse spline di cui sopra ma interpolanti la funzione in otto punti equidistanti nello stesso intervallo.



PIC

Figura 4.11: Questo grafico mostra l’andamento della spline cubica naturale (verde) e di quella not-a-knot (rossa) interpolanti la funzione |x|) in sette punti equidistanti nell’intervallo [5,5]

PIC

Figura 4.12: Questo grafico mostra invece l’andamento delle stesse spline di cui sopra ma interpolanti la funzione in otto punti equidistanti nello stesso intervallo.



PIC

Figura 4.13: Questo grafico mostra l’andamento dell’errore di interpolazione per le funzioni spline interpolanti la funzione 1(1+x2) all’aumentare del numero di punti di interpolazione. Notiamo come l’errore vari molto quando il numero di punti di interpolazione passa da dispari a pari (per valori piccoli): questo causato dal fatto che nel caso di numeri pari viene saltato il punto centrale dell’intervallo in cui la funzione assume valore massmo.

PIC

Figura 4.14: Questo grafico mostra l’andamento dell’errore di interpolazione per le funzioni spline interpolanti la funzione |x| all’aumentare del numero di punti di interpolazione.



PIC

Figura 4.15: Questo grafico mostra l’andamento di una spline completa (rosso) confrontata con quello di una spline not-a-knot (verde) nell’interpolare la funzione xsin(3x) + x.

PIC

Figura 4.16: Questo grafico mostra l’andamento di una spline periodica interpolante la funzione sin(x) in 4 punti equidistanti.



PIC

Figura 4.17: Esempio di approssimazione ai minimi quadrati di una serie di dati.

PIC

Figura 4.18: Esempio di come cambia l’approssimazione all’aumentare progressivo del perso del punto evidenziato. In ordine crescente di peso le funzioni risultanti sono di colore nero (1), blu(2), ciano(3), verde(10).


Capitolo 5
Formule di quadratura

Consideriamo adesso il problema di voler calcolare, in modo approssimato, il valore di un integrale definito:
      ∫ b
I(f) =   f(x)dx,  a < b
       a

con [a,b] intervallo limitato, ed f(x) continua su tale intervallo. Per fare ci andremo a calcolare l’integrale di una funzione polinomiale (o polinomiale a tratti) che approssimi fx), cosa che siamo in grado di fare in modo semplice e senza approssimazioni. Quello che adremo a fare sar quindi calcolare:

      ∫ b        ∫ b
I(f) =  a f(x)dx ≈  a φ(x)dx,  φ(x) ≈ f(x)

Studiamo adesso il condizionamento del problema:

             |∫               |
             ||  b             ||
|I(f)− I(φ)| = || a (f(x)− φ(x))dx||

questa quantit pu essere maggiorata con:

             ∫                                  ∫
|I(f)− I(φ)| ≤  b|f (x) − φ (x )|dx ≤ max |f(x)− φ(x)| bdx
              a                x∈[a,b]            a

         ∫ b
≤ ∥f − φ ∥ a 1dx = (b − a)∥f − φ∥.

Il fattore (b a) rappresenta quindi il fattore di amplificazione dell’errore sul risultato, e quindi:

k = b− a

definisce il numero di condizinamento del problema.

5.1 Formule di Newton-Cotes

Consideriamo l’approssimazione di f(x) fornita dal polinomio interpolante su n + 1 ascisse equidistanti:

p(xi)  =  f(xi) ≡ fi, i = 0,1,...,n;
  xi  =  a+ ih,     h = b−-a-
                         n

Se si considera la forma di Lagrange del polinomio interpolante si ottiene:

      ∫   n               n   ∫             n    ∫    n
I(f) ≈  b∑  f L  (x)dx = ∑  f   bL  (x)dx = ∑ f   b  ∏    x-−-xjdx
       a k=0 k  kn        k=0 k a  kn       k=0 k  a j=0,j⁄=k xi − xj

ponendo x = a + th, t = 0,...,n, si ottiene:

          (                             )        (                  )
  b− a ∑n     ∫ n  n∏   a + th − (a+ jh)       ∑n     ∫ n  ∏n   t− j
= --n--   (fk          a-+-kh−-(a+-jh)dt) = h    (fk          k−-jdt)
       k=0     0 j=0,j⁄=k                      k=0     0  j=0,j⁄=k

Abbiamo quindi ottenuto la formula che definisce l’approssimazione di I(f) che stavamo cercando:

         ∑n
In(f) ≡ h   cknfk
         k=0

in cui

     ∫ n  n∏    t− j
ckn =          k-−-jdt, k = 0,...,n.
      0 j=0,j⁄=k

detta formula di quadrtura di Newton-Cotes.

Teorema

Per i coefficienti ckn vale che:

  n
1-∑
n    ckn = 1
  k=0


Dimostrazione:
considerando f(x) 1 si ha che:

∫ b                  n            (   n    )
   1dx = b− a ≡ b−-a-∑ ckn = (b− a) 1∑   ckn
 a               n  k=0             nk=0

da cui deriva direttamente la tesi.

5.1.1 Formula dei trapezi

Poniamo n = 1. Sappiamo che c01 + c11 = 1, e quindi:

     ∫
c  =   1t−-0dt = 1 =⇒ c  = 1− c  =  1
11    0 1− 0     2     01      11   2

e quindi:

       b− a
I1(f ) =----(f(a)+ f(b))
         2

che approssima l’aria sottesa dal grafico di f(x) con quella sottesa dal trapezio di vertici (a,0),(a,f(a)),(b,f(b))(b,0).

5.1.2 Formula di Simpson

Se invece poniamo n = 2, sappiamo che c02 + c12 + c22 = 2, pertanto:

     ∫ 2t(t− 1)    1 ∫ 2          1 (t3   t2) ||2   1
c22 =   -------dt =-   (t2 − t)dt =  -- − -- || = -
      0   2⋅1      2  0           2   3   2  0   3

e

     ∫ 2(t−-1)(t-−-2    1 ∫ 2                   1   ∫ 2         1
c02 = 0     2⋅1    dt = 2 0 t(t− 1)− 2(t− 1)dt = 3 − 0 (t− 1)dt= 3
                                                   ◟---◝◜----◞
                                                       =0

e si pu ricavare c12 = 4
3. Quindi:

       b− a (         (a + b)      )
I2(f) =--6-- f(a)+ 4f  --2--  +f (b)

Implementazione

Questo il codice che implementa quanto appena descritto:

5.1.3 Condizionamento del problema

Indichiamo con φ(x) la perturbazione della funzione f(x). Si ottiene:

                    | n           |        n
|I (f )− I (φ )| = b−-a||∑  (f  − φ )c ||≤  b−-a∑  |f − φ |⋅|c  |
  n      n       n  ||k=0 k    k kn||    n  k=0  k   k    kn

  (      ∑n     )
≤   b−-a-   |ckn| ∥f − φ∥
      n  k=0

e quindi otteniamo che il numero di condizionamento del problema :

     b−-a-n∑
kn =  n      |ckn|
          k=0

Osserviamo che :

kn = k,  n = 1,...,6;
  kn > k,  n > 6.

Ne consegue quindi che i problemi di calcolo I(f) e In(f) hanno lo stesso numero di condizionamento solo per n 6, e quindi le formule di Newton-Cotes sono sconvenienti al crescere di n.

5.2 Errore e formule composite

Esaminiamo adesso l’errore di quadratura, definito in questo modo:

E  (f ) = I(f)− I (f)
  n            n

Dall’espressione dell’errore di interpolazione discende che:

        ∫ b        ∫ b                 ∫ b
En(f) =  a e(x)dx =  a (f(x)− pn(x))dx =  a f[x0,...,xn,x ]ωn+1(x)dx.

Teorema

Se f(x) C(n+k), con

    { 1,  se n `e dispari,
k =   2,  se n `e pari

allora l’errore di quadratura dato da:

           (n+k)(ξ)(     )n+k+1
En (f) = νnf------- b-−-a
          (n + k)!    n

per un opportuno ξ [a,b], ed:

     { ∫ n∏nj=0(t− j)dt,  se n `e dispari
νn =   ∫0nt∏n  (t− j)dt,  se n `e pari
        0    j=0

Per il metodo dei trapezi si ottiene quindi :

         1
E1(f) = −--f(2)(ξ)(b− a)3,  ξ ∈ [a,b]
         12

e per quello di simpson:

          1  (4)   (b− a )5
E2 (f ) = − 90f (ξ) --2--  ,  ξ ∈ [a,b]

5.2.1 Formule composite

Per ridurre il rapporto (b a)∕n, possibile pensare di suddividere l’intervallo [a,b] in pi sottointervalli della medesima ampiezza, ed utilizzare la stessa dormula di Newton-Cotes per ciascun sottointervallo.

Formula dei trapezi composita

Consideriamo l’applicazione della formula dei trapezi su ciascun sottointervallo [xi1,xi],i = 1,...,n di una partizione uniforme di [a,b]:

I(f)  ≈  I(n)(f) ≡ b−-a(f + f + f + f + ...+ f   + f    + f )
          1   (    2n∑   0   1   1)   2       n−1   n−1    n
      =  b−-a- f0 + 2 ni−=11fi + fn
         b2−na ( ∑n        f0 − fn)
      =  --n-- (  i=0fi)− --2----

Il cui errore di quadratura vale:

  (n)      n  (2)  ( b− a)3
E 1 (f) = − 12f (ξ)  -n---  , ξ ∈ [a,b].

Formula di Simpson composita

Considerando solo valori pari di n:

        (n)         b − a
I(f) ≈ I2  (f) =   −-3n--(f0 + 4f1 + f2 + f2 + ...+ fn−2 + 4fn−1 + fn)
                  b−-a-( ∑n ∕2        ∑n ∕2          )
              =    3n   4  i=1f2i−1 + 2  i=0 f2i − f0 − fn

ed il corrispondente errore di quadratura:

 (n)       n       ( b− a)
E2 (f) = −180f(4)(ξ)  -n--- ,ξ ∈ [a,b]

Notiamo come, per n →∞, Ek(n) 0 con k = 1,2.

Implementazione

Questo codice implementa la formula dei trapezi e quella di Simpson composite:

5.3 Formule adattative

Non infrequente che la funzione integranda f(x) abbia delle variazioni rapide solo in uno o pi piccoli sottointervalli dell’intervallo [a.b], e al contempo sia molto regolare negli altri tratti. Quello che si vuole ottenere in questi casi un metodo che si sappia adattare automaticamente a questo tipo di situazioni, ovvero che permetta di scegliere in modo opportuno i nodi della partizione dell’intervallo [a,b].

5.3.1 Formula dei trapezi

applicando il metodo dei trapezi sull’intervallo [a,b] diviso in 1 e 2 sottointervalli si ha:

      (
      |{  I(1)(f) = b-−-a(f(a) +f (b))
I(f) ≈     1      b 2− a
      |(  I(12)(f) =-----(f(a) +2f (a+2b) +f (b))
                   4

ovvero:

                      (     )
       (1)      -1  ′′    b−-a- 3
I(f )− I1  (f ) ≈ −12 f (ξ)   1

       (2)      -1  ′′  (b−-a)3
I(f )− I1  (f ) ≈ −12 f (ξ)  4

facendo la differenza membro a membro si ottiene:

                               (     )
I(2)(f)− I(1)(f) ≈ − 1-f′′(ξ)(b− a)3 1− 1
 1       1        12             4

e ancora dividendo per 3:

− I1(2)(f) + I1(f)     1  (2)        31         (2)
------3------- ≈ −12f   (ξ)(b− a) 4 ≈ I(f)− I1 (f).

Abbiamo appena trovato una stima dell’errore:

       (2)     1  (2)      (1)
I(f)− I1 (f) ≈ 3(I1 (f)− I1 (f))

che possiamo utilizzare come soglia in un algoritmo ricorsivo per decidere se il caso di suddividere ulteriormente il sottointervallo considerato.

5.3.2 Formula di Simpson

applicando il metodo di Simpson sull’intervallo [a,b] diviso in 2 e 4 sottointervalli si ha:

      (
      ||| I(2)(f) = b−-a(f(a)+ 4f(a+b)+ f(b))
      {  2        6             2
I(f) ≈ ||
      |( I(24)(f) = b−-a(f(x0)+ 4f(x1)+ 2f(x2) +4f (x3)+ f(x4))
                  12

dove x0 = a,x2 = a+b-
 2,x4 = b,x1 = x0+x2
  2,x3 = x2+x4
  2.
Si ottiene quindi:

       (2)       1  (4)  ( b− a)5
I1(f)− I2 (f) ≈ − 180f (ξ)  -2---

       (4)       1       ( (b − a)5)
I1(f)− I2 (f) ≈ − 180f(4)(ξ) 25-⋅16-

da cui:

                   1       (b − a)5            1       (b− a)5
I(42)(f)− I(22)(f) ≈ −---f(4)(ξ)---5--(1− 16) = −---f(4)(ξ)--5----(16− 1)
                  180        2              180       2 ⋅16

Abbiamo appena trovato una stima dell’errore:

              I(4) − I(2)
I(f )− I2(4)(f ) ≈-2----2--
                 15

che come nel caso della formula dei trapezi possiamo utilizzare come soglia per decidere se il caso di suddividere ulteriormente il sottointervallo considerato.

Implementazione

Questi codici implementano in maniera ricorsiva e molto semplice le formule adattive dei trapezi e di Simpson:


PIC

Figura 5.1: Questo grafico mostra il metodo di simpson adattativo sulla funzione sin(10(1 + x2)) + x con tolleranza 103.

PIC

Figura 5.2: Ecco il comportamento sulla stessa funzione del metodo adattivo dei trapezi, ma con tolleranza 5 102. A parit di tolleranza il metodo dei trapezi lavora su 697 punti contro solo 97 del metodo di Simpson del grafico sopra.


Capitolo 6
Altre implementazioni

6.1 Modifica al metodo delle secanti

Questo metodo, invece di richiedere la formula della derivata prima come fa il metodo originale, calcola un’approssimazione del valore di questa nel punto iniziale x0:

In questa tabella confrontiamo i valori ottenuti con i due metodi:

Cercando uno zero della funzione f(x) = x3 + sin(10x):








x0 tolx iter.S ris. Secanti iter.SM ris. Sec. Mod. errore







0.2 102 2 0.31674403271112 2 0.316744038045635e-09







0.2 106 4 0.31735605718908 4 0.317356057189081e-15







0.2 1010 5 0.31735605713204 5 0.317356057132045e-17














1.5 102 15 -0.31053322118703 17 -0.31735615577029 6e-3







1.5 104 18 -0.31735605713204 18 -0.317356057132454e-13







1.5 1010 19 -0.31735605713204 19 -0.317356057132045e-17














1.06 102 47 -8.549970108032228e-06 25 0.31666506455617 -







1.06 104 49 -8.844156346925046e-21 27 0.31735605728114 -







1.061010 49 -8.844156346925046e-21 29 0.31735605713204 -







Vediamo quindi che la valutazione approssimata della derivata porta in un caso a convergere verso una radice diversa rispetto a quella cui converge il metodo originale (ma comunque un risultato corretto), e negli altri casi non peggiora il risultato.


PIC
Figura 6.1: Questo grafico mostra la funzione f(x) = x3 + sin(10x), i cerchi rossi indicano gli x0 utilizzati come punti di partenza dei metodi.


6.2 Fattorizzazione LU

6.2.1 Algoritmo ottimizzato

Questa la versione ottimizzata dell’algoritmo presentato nella parte teorica. Con queste semplici modifiche si riesce ad ottenere un aumento di prestazioni notevole: fattorizzando una matrice di elementi random di 2000x2000 elementi, si passa dai 195.73 secondi dell’algoritmo originale a 74.32.

6.2.2 Ottenere i fattori L ed U

Con questo algoritmo si recuperano dalla matrice A restituita dall’algoritmo di fattorizzazione i fattori L ed U necessari alla effettiva risoluzione del sistema:

6.2.3 Risolvere il sistema

Questo semplice codice automatizza la risoluzione del sistema Ax = b:

6.2.4 Ottimizzazione per sistemi tridiagonali

In caso si debba fattorizzare un sistema tridiagonale si pu migliorare l’algoritmo:

6.3 Fattorizzazione LDLT

6.3.1 Algoritmo ottimizzato

Questa versione ottimizzata dell’algoritmo presentato nella parte teorica utilizza un vettore contenente la porzione triangolare inferiore della matrice A anzich tutta la matrice, sfruttandone la simmetria. Il codice pu essere ulteriormente ottimizzato integrando le funzionalit dei metodi vget e vstore nel codice; viene effettuata una chiamata al metodo per ogni accesso ad un elemento della matrice, e questo appesantisce parecchio dal punto di vista computazionale (ma facilita la lettura del codie).

6.3.2 Trasformare una matrice sdp in un vettore utile all’utilizzo dell’algoritmo ottimizzato

Questo codice serve per vettorizzare per righe la porzione triangolare inferiore di una matrice:

Questa un’ottimizzazione del codice di vettorizzazione:


PIC

Figura 6.2: Questo grafico mostra il costo medio temporale per la fattorizzazione LDL di matrici di dimensioni variabili da 50×50 a 200×200 elementi. In Blu rappresentato il tempo impiegato dalla fattorizzazione con matrice, in rosso la fattorizzazione LDL vettorizzata, in verde il tempo impiegato per trasformare la matrice in vettore. Notiamo come la vettorizazione, bench permetta un grosso risparmio in termini di memoria (memorizziamo solo circa la met degli elementi della matrice) aumenta non di poco il tempo di computazione, pagato in termini di trasferimenti di porzioni di memoria e chiamate a funzioni esterne.

PIC

Figura 6.3: Questo grafico confronta il costo temporale all’aumentare della dimensione della matrice dell’algoritmo di fattorizzazione con le chiamate di funzione (blu), quello senza chiamate (rosso) e quello originale (verde). Notiamo come le chiamate di funzione appesantissero in maniera notevole il metodo. Il giallo indica il tempo di vettorizzazione (non ottimizzato).


6.3.3 Piccolo miglioramento

Questa versione del metodo elimina le chimate di funzione per il calcolo degli indici:


PIC

Figura 6.4: Questo grafico confronta il costo temporale all’aumentare della dimensione della matrice dell’algoritmo di vettorizzazione non ottimizzato (blu) con quello ottimizzato (verde). Questo confronto mostra come l’efficienza aumenti lavorando su blocchi di dati anzich dati singoli.


6.3.4 Ottenere i fattori L e D

Con questo algoritmo si recuperano dalla matrice A restituita dall’algoritmo di fattorizzazione non ottimizzato i fattori L e D necessari alla effettiva risoluzione del sistema:

E con questo si fa la stessa cosa a partire dal vettore restituito dalla versione vettorizzata:

6.3.5 Risolvere il sistema

Questo semplice codice automatizza la risoluzione del sistema Ax = b:

6.4 Fattorizzazione LU con pivoting parziale

6.4.1 Algoritmo ottimizzato

Come nel caso della fattorizzazione LU, svolgendo a mano i cicli che Matlab svolge in maniera inefficiente guadagnamo del tempo:

6.4.2 Risolvere il sistema

Questo codice automatizza la risoluzione del sistema Ax = b:

in cui utilizziamo il metodo di permutazione:

6.5 Fattorizzazione QR

6.5.1 Ottenere i fattori Qed R

Con questo algoritmo recuperiamo i fattori Qed R direttamente utili alla risoluzione del sistema; l’algoritmo di fattorizzazione restituisce nella porzione triangolare inferiore della matrice i vettori di Householder ed pertanto necessario calcolare la matrice Q, R facilmente costruibile a partire da ˆR.

6.5.2 Risolvere il sistema

Anche in questo caso abbiamo un semplice codice per automatizzare la risoluzione del sistema Ax = b:

6.6 Algoritmo di Horner

6.6.1 Algoritmo originale

Questo algoritmo valuta il polinomio p(x) = k=0nakxk nei punti contenuti nel vettore x a partire dal vettore a contenente i coefficienti del polinomio. Quello che si fa infatti svolgere il calcolo: a1 + x(a2 + x(a3 + x(...anx))) che svolgendo i calcoli si vede essere equivalente a: a1 + xa2 + x2a3 + ... + xnan.

6.6.2 Algoritmo greneralizzato

Questo algoritmo valuta il polinomio p(x) = a1+(xx1)a2+(xx1)(xx2)a3+...+(xx1)...(xxn1)an nei punti contenuti nel vettore x0 a partire dal vettore a contenente i coefficienti del polinomio e il vettore x contenente le radici del polinomio.

6.6.3 Polinomio interpolante di Newton

6.6.4 Polinomio interpolante di Hermite

6.7 Spline

6.7.1 Spline lineare

Questo semplice codice calcola e valuta una spline lineare:

6.7.2 Spline cubica naturale

Questi codici effettuano il calcolo e la valutazione in un array di punti di una spline cubica naturale:

6.7.3 Spline cubica not-a-knot

Questi codici effettuano il calcolo e la valutazione in un array di punti di una spline cubica con condizioni not-a-knot:

6.7.4 Spline cubica completa

Questi codici effettuano il calcolo e la valutazione in un array di punti di una spline cubica Completa:

6.7.5 Spline cubica periodica

Questi codici effettuano il calcolo e la valutazione in un array di punti di una spline cubica Periodica:

6.7.6 Approssimazione ai mini quadrati

Qusto codice calcola e valuta il polinomio che risolve il problema dei minimi quadrati:


PIC

Figura 6.5: Ecco un esempio di approssimazione ai minimi quadrati di una serie di dati sperimentali: in questo grafico stiamo approssimando i dati relativi ai tempi di fattorizzazione LDLT relativi al metodo fattLDLlib (blu) e vectorizeLDL (verde).

PIC

Figura 6.6: Approssimazione ai minimi quadrati dei dati relativi al metodo fattLDL (blu) e fattLDLlib (verde).


6.8 Formule di quadratura

6.8.1 Formula dei trapezi adattiva non ricorsiva

6.8.2 Formula di Simpson adattiva non ricorsiva