Indice

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 ∈ℝ approssimazione del risultato esatto x ∈ℝ. É possibile definire la grandezza errore assoluto come:

o in altra forma:

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(10⁻⁵) 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(10⁻⁵).

1.1.2 L’errore relativo

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

da cui:

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 = . L’errore relativo è quindi: ε_x = ≃ = −0.000029492. Notiamo che la formula: −log₁₀(ε_x) restituisce all’incirca il numero di cifre significative corrette all’interno di : 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 x₀:

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

sostituendo si ha:

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

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

ovvero:

In questo caso l’errore commesso è O(h²), 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 {x_n}, ricavabili mediante la funzione di iterazione del metodo:

Ovviamente, affinché il calcolo iterativo possa avere inizio, è necessario definire un x₀ 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:

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:

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

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 x_n = x, x_n+1≠x e non sarebbe vero che lim_n→∞x_n = x, ovvero la funzione non sarebbe convergente e non sarebbe possibile arrivare a un risultato nemmeno dopo infinite iterazioni.

Esempio

Il metodo iterativo

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

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

Dal momento che consideriamo x₀ = 2, la funzione si avvicina a per valori positivi. Quindi:

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

Esempio

Il metodo iterativo:

fornisce una successione di approssimazioni convergente a .
Ecco un esempio di successione, arrestata a un errore di convergenza ≈ 10⁻¹², considerando = 1.414213562374:

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:

per cui, assegnata una base b ∈ℕ:

A questa stringa viene fatto corrispondere il numero intero:

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

ovvero il massimo numero positivo rappresentabile è:

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

1.3 Rappresentazione dei numeri reali

Un numero reale è rappresentato come una stringa del tipo:

dove α₀ ∈{+,−},α_i,β_j ∈{0,1,...,b − 1} con i = 1...m e j = 1...s ed α₁≠0 per la normalizzazione.
La stringa suddetta rappresenta il numero reale:

in cui lo shift ν è fissato e chiameremo le quantità ρ = ±∑ _i=1ⁿα_ib¹⁻ⁱ 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:

ovvero e = (b − 1)∑ _j=1^sb^s−j − ν.

Teorema

Secondo la rappresentazione appena definita, 1 ≤ ρ < b.
Dimostrazione:
Il minimo numero rappresentabile si ha nel caso α₁ = 1,α₂ = 0...α_m = 0: 1.^m−1.
Il massimo invece si avrà nel caso α₁ = ... = α_m = (b − 1): (b − 1)._m−1
In quest’ultimo caso si ha:

che è sicuramente minore di b in quanto (1 − b^−m) < 1. Cvd.

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

1.3.1 I numeri di macchina

Definiamo adesso l’insieme dei numeri di macchina:

In particolare si ha che:

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

Teorema

ha un numero finito di elementi.
Infatti, α₂...α_mβ₁...β_s possono assumere b^s+m−1 valori, α₁ ne può assumere (b − 1) e α₀ invece soltanto 2. Ovvero i numeri rappresentabili sono: 2(b − 1)b^s+m−1 < 2b^s+m.
In particolare |ρ| = α₁α₂...α_m ≥ α₁ ≥ 1, ma anche:
|ρ| = α₁α₂...α_m ≤ (b − 1).^m−1 ≤ b. Cvd.

Funzione floating

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

Troncamento

Sia x = ±α₁α₂...α_mα_m+1...b^e ∈
Possiamo definire:

Arrotondamento

Sia x = ±α₁α₂...α_mα_m+1...b^e ∈
Possiamo definire

Teorema

∀x ∈,x≠0:

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

e quindi:

Osservazione: −log₁₀|ε_x|≃ numero di cifre corrette nella mantissa.
In particolare −log₁₀|ε_x|≤−log₁₀|u| e m = 1 − log₁₀|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 < oppure α_m+1 ≥.
Nel primo:

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

Nel secondo:

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 , e in questo caso ci troveremmo in una condizione di errore. I casi possibili sono:

• |x|≥ r₂: questa situazione è denominata overflow, e la sua gestione (recovery) dipende dal sistema di calcolo utilizzato. Nel nostro caso, lo standard IEEE754 prevede la definizione di una quantità, indicata con Inf,con eventuale segno, che rappresenta la quantità inf.
  
• 0 ≤|x| < r₁: questa situazione è denominata underflow, e la sua gestione si può tipicamente effettuare in due modi:

  1.

  store 0: si pone fl(x) = 0;

  2.

  gradual underflow: permette la denormalizzazione della mantissa e quindi la graduale riduzione del numero di cifre significative. In questo modo si accresce l’insieme dei numeri di macchina aggiungendogli i numeri denormalizzati, tuttavia bisogna tener presente che in questo caso l’errore ε_x aumenta oltre i limiti definiti nel teorema sopra esposto.

Possiamo adesso completare la definizione di fl(x):

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:

• singola precisione: esponente: 8 bit; segno: 1 bit; mantissa: 23 bit;
• doppia precisione: esponente: 11 bit; segno: 1 bit; mantissa: 53 bit;

In quanto i numeri rappresentati sono da intendersi come normalizzati, si è stabilito che α₁ 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: α₂...α_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:

• e ∈ 1,...,254: il dato è normalizzato e ν = 127
• e = 255: rappresentazione di ±inf se f = 0 (si utilizza solo il bit del segno), Not A Number (NaN) se f≠0, tipicamente il risultato di una operazione indefinita (come potrebbero essere 0 ∗ inf,inf −inf,0∕0).
• e = 0: rappresentazione dello zero se f = 0, rappresentazione denormalizzata se f≠0: in questo caso ν = 126.

per quanto riguarda la singola precisione e:

• e ∈ 1,...,2046: il dato è normalizzato e ν = 1023
• e = 2047: rappresentazione di ±inf se f = 0 (si utilizza solo il bit del segno), Not A Number (NaN) se f≠0, tipicamente il risultato di una operazione indefinita.
• e = 0: rappresentazione dello zero se f = 0, rappresentazione denormalizzata se f≠0: in questo caso ν = 1022.

per quanto riguarda la doppia.
In questo modo si ottiene un range di numeri rappresentabili di circa 10⁻³⁸ ÷ 10³⁸ per i numeri in singola precisione e 10⁻³⁰⁸ ÷ 10³⁰⁸ 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:

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

Esempi

Il più grande numero di macchina si avrà nel caso:

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

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:

Quindi secondo lo standard:

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

La precisione di macchina è:

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 10³⁸ per i reali e circa 10¹⁰ 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 ỹ = (); 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(), ovvero considereremo il metodo numerico come eseguito in aritmetica esatta.
Consideriamo quindi le grandezze = 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:

ovvero:

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:

• k  1: gli errori sul risultato sono paragonabili a quelli sul dato iniziale. In questo caso il problema è ben condizionato;
• k >> 1: gli errori sul risultato possono essere molto più grandi di quelli sul dato iniziale. In questo caso il problema è mal condizionato;
• k  u⁻¹: se si utilizza una precisione di macchina u, i risultati saranno privi di significato in quanto affetti dall’errore di rappresentazione, ovvero |ε_x|  u;

1.5.1 Condizionamento delle operazioni elementari

Moltiplicazione: y = x₁x₂

Divisione: y =

Somma algebrica: y = x₁ + x₂

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

• se x₁ ≈ −x₂ la variabile K assume valori che rendono il problema mal condizionato; nel peggiore dei casi si può incorrere nella cancellazione numerica.
• se x₁ e x₂ hanno segno concorde abbiamo K = 1.

Radice: y =

Esempi

Matlab - Errori di rappresentazione

Il codice Matlab:

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 x≠1 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 :

Calcolare eseguendo non è il migliore dei modi: infatti, in caso x = y = r₁,fl(x) = fl(y) = 0,fl() = 0, oppure se x è vicino ad r₂, x² non può essere rappresentato. Questo problema può essere risolto in modo più efficiente eseguendolo invece in questo modo:

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:

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:

• (eps/2+1)-1*(2/eps): ans = 0
• (eps/2+(1-1))*(2/eps): ans = 1

Ovvero non vale la proprietà associativa. Nel caso invece:

si hanno i seguenti risultati:

• (1e300-1e300)*1e300: ans = 0
• (1e300*1e300)-(1e300*1e300): ans = NaN

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 ∗ 10⁻¹ − 1.234 ∗ 10⁻¹ = 1 ∗ 10⁻⁴.
L’errore relativo su y é:

mentre essendo ε_x1 = 3.5 ∗ 10⁻⁴ e ε_x2 = −1 ∗ 10⁻⁴, si ha che

Il problema risulta essere mal condizionato.

Capitolo 2
Radici di una equazione

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:

• f(x) ha un numero finito di soluzioni; ad esempio il caso f(x) = (x − 1)(x − 2).
• f(x) non ha soluzioni reali; ad esempio f(x) = x² + 1.
• f(x) ha un numero infinito di soluzioni; ad esempio f(x) = sin(x).

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

• f(x) è continua nell’intervallo [a,b]
• f(a)f(b) < 0.

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: x₁ = .
A questo punto possono verificarsi solo tre casi:

• f(x₁) = 0: abbiamo trovato la radice e quindi risolto il problema;
• f(a)f(x₁) < 0: la funzione interseca l’asse delle ascisse prima del punto x₁ e si può ripetere il procedimento sull’intervallo [a,x₁]
• f(x₁)f(b) < 0: la funzione interseca l’asse delle ascisse dopo il punto x₁ e si può ripetere il procedimento sull’intervallo [x₁,b].

2.1.1 Criteri di arresto:

Con questo metodo si ha che |x − x₁|≤. Ragionando per induzione, indicando con [a_i,b_i] l’intervallo di confidenza e con x_i = (a_i + b_i)∕2 l’approssimazione della radice all’intervallo i-esimo, si ha che:

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:

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⁽¹⁾, sviluppando in serie di Taylor:

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:

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:

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 log₂(b − a) − log₂(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.

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

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() con l’intervallo di partenza [10⁻³,10⁻¹].

2.1.2 Ordine di convergenza

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

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

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.

ovvero, per induzione:

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 (x₀,f(x₀)) con l’asse delle ascisse, considerando x₀ l’approssimazione iniziale. Consideriamo quindi y = f(x_i) + f′(x_i)(x − x_i), ponendo y a zero si ha che:

che è definita per f′(x₀)≠0. Possiamo estendere il ragionamento a tutti i passi successivi:

Teorema

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

che è ben definita in quanto f′(x_i) è 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:

Osservazione:

Se x_i →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:

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:

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 + δ) ≡; |Φ(x) − Φ(y)|≤ L|x − y|
allora:

1.

x è l’unico punto fisso di Φ in

2.

se x₀ ∈⇒ x_i ∈ ∀i ≥ 0

3.

x_i →x per i →∞

Dim(1).
Supponiamo per assurdo che ∃ | ≡ Φ() ∈ ulteriore punto fisso per Φ, ovviamente x≠. Ma allora si ha che: 0 < |x −| = |Φ(x) − Φ()|≤ L|x −| ma essendo L < 1 si ha che |x −| < |x −|, assurdo.
Dim(2).
Sia x₀ ∈⇐⇒|x − x₀| < δ.
Allora: |x − x_i| < δ e |x − x_i| = |Φ(x) − Φ(x₀)|≤ L|x − x₀| < Lδ < δ
e x_i ∈.
Dim(3).
|x_i −x| = |Φ(x_i−1) − Φ(x)|≤ L|x_i−1 −x|, ovvero:
|x_i −x| = L|Φ(x_i−2) − Φ(x)|≤ L²|x_i−2 −x|
Ragionando per induzione possiamo quindi concludere che:
|x_i −x|≤ Lⁱ|x₀ −x|, che tende a zero per i che tende a infinito.

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

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

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 |e_i|≤ 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:

e quindi un controllo del tipo |x_i+1 − x_i| < tolx è appropriato ed equivalente al controllo |f(x_i)|≤|f′(x_i)|tolx in quanto x_i+1 = x_i −. 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 |e_i|≤ rtol|x|≈ rtol|x_i+1|, ed un criterio di arresto del tipo:

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|x_i+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 |x_i+1 − x_i| < tolx.

Convergenza lineare

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

Noi vogliamo |e_i| < tolx e quindi:

Ovvero abbiamo ottenuto un criterio tanto più efficiente quanto più la molteplicità della radice è alta (c ≈ 1). Noi conosciamo x_i+1, e vogliamo che |e_i+1|≤ tolx ma sappiamo che |e_i+1| = c|e_i| e possiamo quindi controllare |x_i+1 − x_i|≤ tolx. Il problema è che non conosciamo c. Possiamo però calcolare le prime tre iterazioni (x₀,x₁,x₂) in questo modo: |x₁ −x₀|≈ (1 −c)|e₀|; |x₂ −x₁|≈ (1 −c)|e₁|≈ (1 −c)c|e₀|; e quindi c ≈.

Implementazione

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

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:

Osservazione:

Se f(x) è sviluppabile in serie di Taylor in x, sua radice di molteplicità esatta m, allora è scrivibile nella forma f(x) = (x−x)^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 (x−x)^m. Avendo x₀ approssimazione di x:

ovvero

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 e_i ≈ ce_i−1 ed e_i+1 ≈ ce_i, con c costante asintotica dell’errore, incognita. Dalla combinazione delle due si ottiene che e_i+1e_i−1 ≈ e_i² e quindi (x_i+1 −x)(x_i−1 −x) ≈ (x_i −x)². Se si suppone esatta l’uguaglianza, si ha che:

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 x₀⁽⁰⁾ si calcolano x₁⁽⁰⁾,x₂⁽⁰⁾ approssimazioni di Newton alla prima iterazione,

2.

Livello esterno in cui, con la formula:

si ottiene una successiva approssimazione più accurata x₀⁽¹⁾ 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:

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:

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.

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 x_i−1,x_i,f(x_i−1),f(x_i), si approssima f′(x) in questo modo:

Il metodo iterativo risultante:

con x₀,x₁ 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 x_n è molto vicino a x_n−1, 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 = ≈ 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.

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:

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

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:

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

Capitolo 3
Risoluzione di sistemi lineari

dove A = (a_ij) ∈ℝ^m×n matrice dei coefficienti, b = (b_i) ∈ℝ^m vettore dei termini noti e x = (x_i) ∈ℝⁿ 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:

3.1.1 Matrici diagonali

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

in cui valgono le proprietà:

• a_ij = 0 ⇔ i≠j.
• det(A) = ∏ _i=1ⁿa_ii≠0a_ii≠0 ∀i,0 ≤ i ≤ n.

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

e quindi la soluzione ottenibile:

è ben posta in quanto come detto prima a_ii≠0 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:

• Triangolari inferiori: a_ij = 0 se j > i

• Triangolari superiori: a_ij = 0 se i > j

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

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

In quanto A è nonsingolare, si ha che a_ii≠0, 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:

posizioni di memoria. Per quanto riguarda il costo computazionale invece, da quanto detto sopra si ha che sono necessari 1flop per calcolare x₁,3flop per calcolare x₂,...,2n− 1flop per calcolare x_n, e quindi:

3.1.3 Matrici ortogonali

Una matrice è ortogonale se e solo se vale la proprietà: A⁻¹ = A^T. In questo caso quindi la soluzione del sistema lineare è ottenibile facilmente calcolando x = A^Tb, al costo di un prodotto matrice-vettore; saranno richiesti pertanto circa 2n²flops ed n² 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 L₁,L₂ matrici triangolari inferiori, ed indichiamo con (a_ij) e (b_ij i corrispettivi elementi.
Per ipotesi abbiamo che a_ij = b_ij = 0, j > i.
Allora gli elementi della matrice somma saranno composti:
L₁ + L₂ = (a_ij + b_ij) ma se j > i si ha (a_ij + b_ij) = 0 + 0 = 0, ovvero L₁ + L₂ è triangolare inferiore.
Vediamo adesso il caso del prodotto.
Chiamaiamo c_ij gli elementi della matrice L₁L₂. Vogliamo dimostrare che c_ij = 0, se j > i.
L’elemento c_ij, quando j > i, è scrivibile nella forma: e_i^T(L₁L₂)e_j. Svolgendo le operazioni in quest’ordine: (e_i^TL₁)(L₂e_j) si ha:

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 c_ij = 0.

Teorema

Se L₁ = (a_ij)eL₂ = (b_ij) sono matrici triangolari inferiori, allora gli elementi diagonali della matrice L₁L₂ sono dati da (a_iib_ii).
Dimostrazione:
Sia c_ij elemento della matrice risultante. Per definizione di prodotto tra matrici, si ha che

Ricordando che a_ij = b_ij = 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 L₂ o L₁; 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 a_kk e b_kk. 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

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:

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 L⁻¹ è triangolare dello stesso tipo.
Dimostrazione:
Supponiamo L triangolare inferiore e L⁻¹ triangolare superiore. Allora si avrebbe che (L⋅L⁻¹)(1,2) = L(1,1) ∗L⁻¹(1,1) + 0 + 0≠0, in quanto per ipotesi L(1,1),L⁻¹(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 = F₁F₂...F_k,dove F_k ∈ℝ^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:

É 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 = L₁U₁ = L₂U₂ due fattorizzazioni di A; si ha che:

U₂ è quindi nonsingolare e si ha che L₁⁻¹L₂ = U₁U₂⁻¹ ≡ D. Ma L₁⁻¹L₂ e U₁U₂⁻¹ sono rispettivamente triangolare inferiore e triangolare superiore, pertanto D deve essere diagonale. Ma L₁⁻¹L₂ è una matrice a diagonale unitaria e questo implica che D sia la matrice Identità; ne consegue quindi che L₁ = L₂ e U₁ = U₂. cvd.
Rimane adesso da stabilire l’esistenza di tale fattorizzazione, e procederemo illustrando un metodo di fattorizzazione.
Sia v = (v₁...v_n)^T ∈ℝⁿ vettore tale che v_k≠0, 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:

grazie al fatto di avere v_k≠0, e quindi la matrice elementare di Gauss:

Ovvero:

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: L⁻¹ = I + ge_k^T, infatti:

in quando e_k^Tg,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:

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 a₁₁⁽¹⁾≠0, possiamo definire il primo vettore elementare di Gauss:

e la prima matrice elementare di Gauss:

per avere:

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

A questo punto solo se a_ii⁽ⁱ⁾≠0 sarà possibile definire l’i-esimo vettore di Gauss:

e quindi la i-esima matrice elementare di Gauss:

e quindi:

Se si osserva che nella i-esima matrice di Gauss (L_i ≡ I −g_ie_i^T)) le prime i righe corrispondono con quelle della matrice identità, si può notare che le prime i righe delle matrici A⁽ⁱ⁾ e A⁽ⁱ⁺¹⁾ coincidono e gli elementi azzerati al passo i-esimo (che stanno nelle prime i − 1 colonne) non vengono alterati dalla moltiplicazione per L_i.
Se risulta possibile iterare i passi fin qui descritti i = n − 1 volte, si ottiene che:

Adesso, osservando che la matrice L_n−1...L₁ è il prodotto di matrici triangolari inferiori a diagonale unitaria, è anch’essa di questo tipo, e così sarà anche la sua inversa. Tenendo conto che:

e ponendo formalmente la matrice L_n−1...L₁ uguale a L⁻¹:

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 a_ii⁽ⁱ⁾≠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:

In particolare si ha che il determinante det(A_k) è 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⁽ⁿ⁾ ≡ 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⁽¹⁾, 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:

• l’i-esimo vettore di Gauss
• il prodotto L_iA⁽ⁱ⁾ = (I −g_ie_i^T)A⁽ⁱ⁾ = A⁽¹⁾ −g_i(e_i^TA⁽ⁱ⁾) ≡ A⁽ⁱ⁺¹⁾.

Considerando che le prime i componenti del vettore g_i sono nulle, che il vettore e_i^TA⁽ⁱ⁾ corrisponde alla i-esima riga di A⁽ⁱ⁾ 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⁽ⁱ⁺¹⁾ 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)² sottrazioni ed altrettante moltiplicazioni per aggiornare la matrice; in totale la fattorizzazione avrà un costo approssimativo di:

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 A_k, ∀k = 1,...n, e sono le matrici a diagonale dominante.

Definizione

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

e a diagonale dominante per colonne se:

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 A^T è 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 x≠0 tale che Ax = 0. Assumiamo

Ma se Ax = 0 allora e_k^TAx = e_k^T0 = 0, ne consegue che:

ma allora A non è diagonale dominante. Assurdo.

Lemma

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

3.5 Fattorizzazione LDL^T

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 = A^T) e se vale:

Lemma

Se A è sdp allora a_ii > 0, ∀i = 1,...,n.
Dimostrazione:
Infati, se A è sdp si dee avere che a_ii = e_i^TAe_i > 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 x≠0 tale che Ax = 0, ovvero x^TAx = 0 ovvero A non è sdp. Cvd.

Teorema

A è sdp se e solo se ∀k = 1,...,n, A_k è sdp.
Dimostrazione:
Procediamo per assurdo, supponendo A_k non sdp. Allora deve esistere y ∈ℝ^k,y≠0, tale che y^TA_ky ≤ 0. Allora deve esistere anche x ∈ℝⁿ tale che x = (y,0..0)^T≠0, che soddisfi 0 < x^TAx in quanto A è sdp per ipotesi. Ma allora avremmo: