Appunti completi di Metodi numerici - Parte 3

T

La H H ha la dimensione (n+1)*(n+1) e questo è fantastico perché n era la dimensione piccola, quindi la

T

matrice ha una dimensione piccola. La H H è una matrice simmetrica e definita positiva. Potrei pensare a

questo punto di applicare il Cholesky e avrei finito però purtroppo non si può fare in questo caso.

171 +

Se la matrice è rettangolare, MatLab non prende l’inversa ma prende l’inversa generalizzata (H ).

Se invece avessimo una matrice quadrata:

L’inversa generalizzata mi rida fuori l’inversa. Abbiamo visto queste cose perché quando usiamo MatLab e

usiamo \, se utilizziamo una rettangolare il comando ci fa la soluzione del sistema, ma ce la fa ai minimi

quadrati ovvero risolvendo le equazioni normali ovvero calcolando l’inversa generalizzata, quindi prende la

soluzione che mi dà l’errore più piccolo e non quello che mi fa venire zero. Quindi se il sistema è quadrato il

\ ci dà la nostra soluzione tradizionale, ma se diamo una matrice rettangolare ci risolve dandoci una inversa

molto più generale.

La soluzione a questo problema la si trova risolvendo le cosiddette equazioni normali. Attenzione al

MALCONDIZIONAMENTO:

T 2

K(H H)=K(H)

Il risultato di questo è che errori di arrotondamento nella risoluzione del sistema possono causare grandi

errori nell’approssimante dei dati. Se non ci piace questa strada perché ha un indice di condizionamento

molto brutto, non posso risolvere le equazioni normali, in generale non andrei mai a mettere un Cholesky

che mi risolve le equazioni normali perché è pericoloso, è meglio avere delle alternative.

172

Polinomio di m.a. ai minimi quadrati

In questo caso abbiamo dei passaggi, usando un esempio polinomiale, per vedere come si può ottenere

l’equazione normale. Si tratta di minimizzare:

Derivando rispetto ad a0 , a1 , …, an ed uguagliando a zero (sistema di equazioni normali):

Si ottiene un sistema lineare.

In forma matriciale:

T

La matrice H H risulta:

• Simmetrica per costruzione.

• Invertibile e definita positiva.

• Può essere malcondizionata. 173

Esempio

174

175

176

Metodo QR-LS

Si va a vedere la soluzione numerica.

Si ha l’errore indicato con S, il risultato non verrà zero quindi mi accontento di renderlo il più piccolo possibile,

quindi tra tutti gli a si va a prendere quello che lo minimizza. Si applica una fattorizzazione QR alla matrice

del sistema H=QR. La Q è una matrice quadrata con la dimensione grande ed R è una dimensione

rettangolare, fatta esattamente come H, però la R è triangolare (sopra piena e sotto vuota):

La Q ortogonale è tra le matrici più amate perché era quella che soddisfava la relazione:

T T

Q Q = QQ = I

Questo significa che la trasposta coincide con l’inversa:

T -1

Q = Q

Quindi le matrici ortogonali sono deliziose perché per fare l’inversa basta scambiare le righe con le colonne.

R1 triangolare superiore non singolare. Questa fattorizzazione me la calcolo e al posto di H andiamo a scrivere

QR: 2

‖ ‖

−

min 2 T

Ora vado a inserire all’interno della norma la matrice Q :

2

‖ ‖

−

min 2

Osservazione: Si può andare ad inserire all’interno della norma una matrice?

Se abbiamo un vettore qualsiasi, andare a fare la norma 2 al quadrato vuol dire fare il prodotto scalare del

vettore con sé stesso:

2 T

‖ ‖ = < v, v > = v v

2

Se andiamo a inserire la matrice all’interno della norma:

2 T T T

‖ ‖ = < Pv, Pv > = v P P v = v v

2

Se abbiamo usato una matrice ortogonale, la lunghezza del vettore non cambia. P, le trasformazioni lineari

rappresentate da matrici ortogonali sono le cosiddette isometrie, che non cambiano la lunghezza dei vettori.

Quindi andando a inserire la abbiamo fatto un’operazione lecita.

177

Si partiziona il vettore c:

Per minimizzare questo vado a risolvere un sistema triangolare piccolo, trovo gli a e ho finito, quindi il

problema difficilissimo con le equazioni normali l’ho trasformato nella soluzione di un piccolo triangolare.

Si risolve il sistema triangolare:

Il residuo è dato da: Scelta delle funzioni base

Scelta dei polinomi nella base delle potenze:

si vuole determinare:

che renda minimo:

Come per l’interpolazione, anche per l’approssimazione ci piace tanto il polinomio, ma può essere che nei

casi reali si abbia bisogno di funzioni molto più complicate, quindi il problema è scegliere le funzioni di base,

costruirsi la matrice dei coefficienti (ovvero la H). Però attenzione perché il ragionamento dobbiamo farlo

noi, ovvero scegliere le funzioni di base del modello, richiamare la fattorizzazione QR che ci dà MatLab

dopodiché si risolve un semplice triangolare e abbiamo risolto il problema.

178

CAPITOLO 9 - Formule di quadratura

In questo capitolo ci occuperemo del calcolo di integrali definiti da a a b, di f(x) una funzione continua nel

caso monodimensionale. Le formule di quadratura che andremo ad utilizzare si dividono in:

▪ Formule di Newton-Cotes (monodimensionali)

o Trapezi

o Simpson

o 3/8

o Punto medio

o Composite

▪ Formule di Gauss (bidimensionali), che utilizzeremo nella seconda parte di corso e sono le più efficaci

per codici ad elementi finiti

Considerazioni teoriche

Dato l’integrale esatto da [a,b] della funzione f(x) continua nell’intervallo, si può avere anche la sua

approssimazione I perché se non so risolvere la funzione f(x) posso andarla a sostituire con qualcosa di più

n

semplice f (x), che è un modellino che dipende da n parametri, e che condivide con la funzione un certo

n

numero di valori. Allora la migliore sostituzione della funzione f(x) è con un polinomio interpolante, cioè che

condivide con f un certo numero di punti e di valori, ed inoltre la risoluzione di un polinomio è abbastanza

semplice (banale). Approssimazione di f(x)

⇒

Per valutare se la sostituzione che si sta facendo è opportuna e giusta, ci viene in soccorso la teoria, che ci

dice che la cosa migliore da fare è andare a valutare l’errore. L’errore è la differenza tra l’integrale esatto e

l’integrale approssimato, e questo errore dovrà tendere allo zero.

In valore assoluto l’errore è minore uguale dell’espressione dell’integrale in valore assoluto, e poi posso

maggiorare con la grandezza dell’intervallo moltiplicato per il massimo dell’errore. La teoria dice che se riesco

a rendere piccolo questo valore massimo, ovvero calcolando la norma infinito di questo valore risulta minore

di un valore ε, allora posso dire che l’errore è ε volte (b-a).

L’analisi ci dice che se metto un polinomio al posto di fn, al cresce del grado n, posso abbattere l’errore sino

al raggiungere l’ε desiderato. 179

La formula dell’integrale approssimato con il polinomio è:

Il polinomio interpolante di Lagrange è una somma dei valori disponibili, in questo caso i campioni della

funzione (x ) che si trovano nei nodi di interpolazione, ciascuno moltiplicato per il polinomio fondamentale di

i

Lagrange (L , L , L ,…, L ) . Poi andando a risolvere l’integrale dei polinomi ricavo dei numeri perché sono degli

0 1 2 n

integrali definiti in un intervallo [a,b] che saranno indicati con w (weight) o con c (coefficienti).

i i

Quindi il calcolo dell’integrale approssimato è semplicissimo perché dato da una sommatoria di coefficienti

e campionamenti della funzione valutata nei nodi.

Inoltre, dovremo decidere quanti nodi prendere, ovvero di che grado n deve essere il polinomio da prendere.

Le formule di quadratura sono semplicissime, devo solo campionare la funzione integranda, devo dire quanti

sono i nodi opportuni e quanti ne voglio.

Grado di precisione o esattezza

Un modo di caratterizzare l'accuratezza di una formula di quadratura è tramite il suo grado di precisione, il

massimo intero m tale che tutti i polinomi grado non superiore ad m siano integrati esattamente:

La perfezione si avrebbe se l’integrale esatto è uguale all’integrale approssimato, perché come detto la

differenza è l’errore E. La precisione è legata alla scelta di n, ovvero il grado del polinomio, e il grado di

precisione è definito come il massimo numero m che fa venire zero, cioè il massimo grado del polinomio

interpolante che mi fa venire l’errore uguale a zero per ogni funzione polinomiale.

- Se n è il grado del polinomio interpolatore, E ≡ 0 per k = 0,…,n, e di conseguenza m ≥ n.

k

- E' pero possibile che anche se l'errore di interpolazione non è identicamente nullo, e di qui

la possibilità di avere m > n.

Come si fa a fare la verifica di queste funzioni polinomiali e ad inventare i polinomi?

I polinomi costituiscono uno spazio vettoriale, quindi un qualsiasi polinomio lo posso costruire basandomi su

2 m

una base. Allora provo la formula di quadratura sulla base dei monomi su 1, su x, su x , fino a x , finché

m+1

l’errore viene zero, e dico che ha grado precisione 1, 2, …, m. Poi quando per x l’errore non viene più zero

mi fermo perché oltre quel grado di precisione non posso andarci.

Se utilizzo il polinomio interpolante di Lagrange o di Newton che grado di precisione hanno?

180

+1

()) ) ( )

= − → ∫ (() − → ∫ (( − ∗ … ∗ − ∗ )

0 ( + 1)!

Partendo dalla rappresentazione dell’errore del polinomio interpolante dobbiamo determinare quando

questo errore è pari a 0. Se costruisco l’integrale approssimato su dei polinomi e la funzione è un polinomio

allora dovrebbe venirmi zero. Quindi se f è un polinomio questa formula viene zero fino a che f è un polinomio

di grado n, perché per il polinomio n+1 non è detto che l’errore viene zero.

Da cosa dipende l’errore di interpolazione?

- Regolarità della funzione

- Disposizione dei punti di interpolazione sull’asse delle ascisse

Formule di quadratura di Newton-Cotes nodi equispaziati

• Newton-Cotes Chiuse: Usano entrambi gli estremi di integrazione

o formula dei Trapezi : Lineare

o formula di Simpson 1/3: Quadratica

o formula di Simpson 3/8 : Cubica

o formula di Boole : 4°-ordine

• Newton-