Calcolo numerico

* Appunti di matematica numerica

Matteo Rorato, Politecnico di Milano

24 agosto 2020

Indice

1 I fondamenti della matematica numerica 3

1.1 Stabilità di un problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Stabilità e convergenza di un metodo numerico . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 I numeri floating point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.4 Errori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Zeri di funzioni non lineari 8

2.1 Il metodo di bisezione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Iterazioni di punto fisso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Il metodo delle corde, di Newton, delle secanti . . . . . . . . . . . . . . . . . . . . . . . . 11

3 Sistemi lineari 14

3.1 Richiami di algebra lineare . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.2 Metodi diretti: il MEG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.3 Analisi di stabilità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.4 Metodi iterativi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.4.1 Metodo di Richardson e del gradiente . . . . . . . . . . . . . . . . . . . . . . . . 26

3.5 Sistemi indeterminati . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4 Autovalori ed autovettori 28

4.1 Il metodo delle potenze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4.1.1 Localizzazione geometrica degli autovalori . . . . . . . . . . . . . . . . . . . . . . 33

4.2 Metodi basati sulle iterazioni QR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

* Nel redigere i presenti appunti è stato preso liberamente spunto dalle lezioni del corso di Calcolo Numerico tenuto dalla prof. Paola F. Antonietti nell’ambito del corso di laurea triennale in Ingegneria Aerospaziale (aa 2014/15), oltre che dai testi citati in bibliografia. 1

5 Approssimazione di funzioni e di dati 35

5.1 Interpolazione polinomiale di Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

5.2 Approssimazione nel senso dei minimi quadrati . . . . . . . . . . . . . . . . . . . . . . . 38

6 Approssimazione delle derivate 39

7 Integrazione numerica 42

7.1 Formule di quadratura interpolatorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

7.2 Formule di quadratura interpolatorie gaussiane . . . . . . . . . . . . . . . . . . . . . . . 47

7.3 Integrazione automatica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

8 Equazioni differenziali ordinarie 53

8.1 Analisi di convergenza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

8.2 Zero-stabilità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

8.3 Assoluta stabilità . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

8.4 Metodi a più passi (o multistep) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

8.5 Metodi Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

8.6 Sistemi di ODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

9 Problemi ai limiti 67

9.1 Approssimazione alle differenze finite di un problema di diffusione-trasporto a trasporto dominante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

9.2 Approssimazione alle differenze finite dell’equazione del calore monodimensionale . . . . 72

A Equazioni a derivate parziali 75

2

1 I fondamenti della matematica numerica

1.1 Stabilità di un problema

Si consideri il seguente problema astratto: trovare tale che x (1.1) F (x, d) = 0 dove è l’insieme dei dati da cui dipende la soluzione, ed esprime la relazione funzionale tra ed F x.

A seconda del tipo di problema rappresentato nella (1.1), le variabili e potranno esprimere d. x d numeri reali, vettori o funzioni. Tipicamente (1.1) viene detto problema se e sono dati e diretto F d x è incognito, se ed sono noti e è incognito, di nel caso in cui e sono inverso identificazione F x d x d dati, mentre la relazione funzionale è incognita (problemi di quest’ultimo tipo esulano dagli interessi F di questo corso).

Il problema (1.1) è o (rispetto a un certo insieme di dati) se ammette un’unica ben posto stabile soluzione la quale dipende con continuità dai dati. x Un problema che non goda della proprietà precedente si dice o e prima di affrontarne la risoluzione numerica è bene, quando mal posto instabile ha senso, regolarizzarlo ovvero trasformarlo in modo opportuno in un problema ben posto. Non è infatti appropriato pretendere che sia il metodo numerico a porre rimedio alle patologie di un problema intrinsecamente mal posto.

Sia l’insieme dei dati ammissibili, ovvero l’insieme dei dati in corrispondenza dei quali il problema D (1.1) ammette soluzione unica (ma non necessariamente dipendenza continua dai dati). La dipendenza continua dai dati significa che piccole perturbazioni sui dati debbano riflettersi in “piccole” variazioni d nella soluzione Precisamente, sia un fissato elemento di e una perturbazione ammissibile sui x. d D d dati tale che e sia la conseguente variazione nella soluzione in modo che si abbia 2 d + d D, x (1.2) F (x + x, d + d) = 0.

Allora, richiediamo che valga la seguente proprietà tale che (1.3) 9K 8 2 k  k= K (d) d : d + d D, xk K dk0 0 0 (K sarà il numero di condizionamento). La norma usata per i dati e quella per le soluzioni possono 0 non coincidere qualora e denotino diversi tipi di variabile (ad esempio una matrice ed un vettore). d x

Esempio. Un semplice esempio di problema mal posto è quello della determinazione del numero di radici reali di un polinomio. Ad esempio, il polinomio presenta, al 4 2 p(x) = x x (2a 1) + a(a 1) variare con continuità di nei numeri reali, una variazione discontinua del numero di radici reali. Si a hanno infatti 4 radici reali se 2 se e nessuna se 2a 1, a [0, 1) a < 0.

Al fine di esprimere la misura della dipendenza continua dai dati, introduciamo la seguente definizione.

Definizione 1.1 di condizionamento di (1.1)). (Numero Per il problema (1.1), definiamo il come numero di condizionamento relativo → k xk/kxk (1.4) 6 2 K(d) = sup , d = 0, d + d D k dk/kdk.

Nel caso in cui o sarà più sensato introdurre il dato numero di condizionamento assoluto, d = 0 x = 0 da → k xk (1.5) 6 2 K (d) = sup , d = 0, d + d D abs k dk 3.

Il problema (1.1) si dirà mal condizionato in se è “grande”; mal condizionato in generale se lo d K(d) è in corrispondenza di ogni dato ammissibile (il significato preciso da attribuire a piccolo e grande d varierà da problema a problema).

La proprietà di buon condizionamento di un problema prescinde dal metodo numerico che verrà usato per risolverlo. In effetti, è possibile costruire metodi numerici sia stabili, sia instabili per risolvere problemi ben condizionati. Il concetto di stabilità per un algoritmo o per un metodo numerico è analogo a quello usato per il problema (1.1) e verrà introdotto nella prossima sezione.

1.2 Stabilità e convergenza di un metodo numerico

Supporremo nel seguito che il problema (1.1) sia ben posto. Un metodo numerico per la risoluzione approssimata di (1.1) consisterà, in generale, nel costruire una successione di problemi approssimati (1.6) F (x , d ) = 0, n 1 n n n dipendenti da un certo parametro (che definiremo caso per caso), con la sottintesa speranza che n per ovvero che la soluzione numerica converga alla soluzione esatta. Affinché questo! ! 1,x x n n avvenga è necessario che e che “approssimi” quando Precisamente, se il dato! ! 1.d d F F n d n n del problema (1.1) è ammissibile per , si dice che (1.6) è se consistente F n per (1.7)! ! 1 F (x, d) = F (x, d) F (x, d) 0 n n n essendo la soluzione di (1.1) corrispondente al dato Il senso di questa definizione verrà precisato x d. nei capitoli seguenti per ogni singola classe di problemi.

In analogia a quanto detto per il problema (1.1), affinché il metodo numerico sia a sua volta ben posto o richiederemo che per ogni fissato, la soluzione in corrispondenza del dato , che il stabile esista n x d n n calcolo di in funzione di sia (o riproducibile, ovvero se ripetuto più volte con lo stesso dato, unico x d n n il calcolo deve dare la stessa soluzione), che Più precisamente, sia dipenda con continuità dai dati. x n un elemento arbitrario di , essendo l’insieme di tutti i dati ammissibili per (1.6). Sia una d D D d n n n n perturbazione ammissibile, ovvero tale che , e sia la corrispondente perturbazione 2 d + d D x n n n n sulla soluzione, in modo tale che si abbia (1.8) F (x + x , d + d ) = 0 n n n n n.

Allora, richiediamo che tale che (1.9) 9K 8 2 k k  k k= K (d ) d : d + d D , x K d0 0 n n n n n n 0 n.

Definizione 1.2 di condizionamento di (1.6)). (Numero Per ognuno dei problemi (1.6) definiamo ed come numero di condizionamento relativo assoluto → k k/kx kx n n (1.10) 6 2 K (d ) = sup , d = 0, d + d D n n n n n n k k/kd kd n n → k kx n (1.11) 6 2 K (d ) = sup , d = 0, d + d D abs,n n n n n n k kd n.

Il metodo numerico è detto ben condizionato se è “piccolo” in corrispondenza di ogni dato K (n ) d n n n ammissibile, mal condizionato in caso contrario. 4

Definizione 1.3 1). (Convergenza Il metodo numerico (1.6) è se e solo se convergente tale che 8" 9 9> 0 n ("), (n , ") > 00 0 8n 8 k k  ) kx(d) > n ("), d : d x (d + d )k "0 n n n n dove è un dato ammissibile per il problema (1.1), la soluzione ad esso corrispondente ed x(d) x (d+ d ) n n la soluzione del problema numerico (1.6) con dato . d + d n.

Misure della convergenza di ad sono fornite dall’errore assoluto o da quello relativo, come vedremo x x n in seguito.

I concetti di stabilità e convergenza sono fortemente interconnessi, come provato dal seguente teorema.

Teorema 1.1 equivalenza o di Lax-Richtmyer). (di Se il problema (1.1) è ben posto, condizione necessaria affinché il problema numerico (1.6) sia convergente è che esso sia stabile, condizione sufficiente affinché il problema numerico (1.6) sia convergente è che sia consistente con il problema (1.1). consistenza convergenza stabilità! !

Si vede dunque come per un metodo numerico consistente, la stabilità è equivalente alla convergenza.

L’analisi di stabilità di un metodo numerico si può effettuare seguendo strategie alternative.

• Nell’analisi (forward si stimano le variazioni sulla soluzione dovute sia ai in avanti analysis) k kx n perturbazioni nei dati, sia ad errori intrinseci al metodo numerico.
• Nell’analisi (backward si stimano le perturbazioni che si dovrebbero “impri-all’indietro analysis) d n mere” ai dati di un determinato problema al fine di riottenere il risultato effettivamente calcolato x̂ n nell’ipotesi in cui si operi in aritmetica esatta (trovare tale che Nell’ef-d F (x̂ , d + d ) = 0). n n n n n fettuare questa stima si tiene conto in alcun modo di si sia ottenuto (ovvero di quale non come x̂ n metodo si sia impiegato per generarlo).

L’analisi in avanti e quella all’indietro sono due diverse modalità della cosiddetta analisi a priori. Essa può essere applicata per indagare non solo la stabilità di un metodo numerico, ma anche la convergenza di quest’ultimo alla soluzione del problema esatto (1.1). Si parlerà in questo caso di analisi a priori e potrà essere ancora realizzata con la tecnica in avanti o con quella all’indietro. dell’errore,

Essa si distingue dalla cosiddetta la quale mira a fornire una stima analisi a posteriori dell’errore, dell’errore sulla base di quantità effettivamente calcolate usando uno specifico metodo numerico.

1.3 I numeri floating point

Vedi quaderno.

1.4 Errori

Il processo di risoluzione di un problema fisico porta alla generazione di due principali tipologie di errore. All’inizio del processo un primo errore sarà commesso nel tentativo di far corrispondere e m alla realtà fisica un certo modello matematico (errore dato dal modello). Successivamente si avrà un errore causato dal processo di calcolo vero e proprio (errore computazionale), il cui controllo e c è quindi oggetto del calcolo scientifico. In fact, when a computed result of a given real number is not machine representable, then there are two ways it can be represented in the machine. The first method, (troncamento), is the method in which the digits after are symply chopped off. chopping d t 5 The second method is (arrotondamento), in which the digits after are not only chopped rounding d t off, but the digit is also rounded up or down with a certain criteria (da [NLA, p.46]). La scelta del d t metodo è dettata dal tipo di operazione. Di conseguenza l’errore computazionale è in genere dato dalla somma di un errore di troncamento e di un errore di arrotondamento . D’ora in poi chiameremo e et a l’errore computazionale semplicemente errore e lo indicheremo con e.

Definizione 1.4 (Errore). Indichiamo con (x) la soluzione esatta del modello matematico e con x (x̂) la soluzione ottenuta al termine del processo numerico. Si definisce errore (computazionale) la x̂ quantità data da (1.12) e = x x̂ e = x x̂ definiamo errore la quantità data da assoluto (1.13) |x |xe = x̂| e = x̂| definiamo (se errore o - (se moltiplicato esprime l’errore in percentuale relativo percentuale 6 → 100)x = 0) di quantità esatta - la quantità data da |x |xx̂| x̂| (1.14) e = e = |x| |x| definiamo norma dell’errore (se la quantità data da relativo 6 x = 0) kx x̂k (1.15) e = kxk dove il simbolo indica una qualsiasi norma vettoriale. La norma dell’errore definito in (1.14) è k.k kek diversa dalla norma dell’errore definita in (1.15). Ricordiamo che |x |x̂ kx kx̂ x̂| = x|, x̂k = xk.

1.1. Rispetto all’errore assoluto, l’errore relativo rapporta la qualità dell’approssimazione Osservazione all’ordine di grandezza della soluzione particolare, risultando quindi spesso più utile. Ad esempio l’errore assoluto tra i valori e e tra i valori e è lo stesso e pari 1, x = 1 x̂ = 2 x = 1001 x̂ = 1002 sebbene l’errore che si ha nel primo caso sia evidentemente più grave a causa del diverso ordine di grandezza, di cui tiene invece conto l’errore relativo. Ovvero nel primo caso ho un errore del 100% rispetto al valore esatto, nel secondo caso dello 0.1%.

1.2. L’errore non è in genere una quantità direttamente calcolabile, in quanto dipen-Osservazione de dall’incognita stessa del problema. Al suo posto vengono perciò utilizzati i cosiddetti stimatori dell’errore.

Definizione 1.5 2). (Convergenza Un metodo iterativo è detto convergente se, detto il valore x k calcolato al passo k e il valore esatto soluzione del nostro problema x (1.16) lim e = lim (x x) = 0 k k k!1 k!1 (notiamo che implica che lim (x x) = 0 lim (x x ) = 0) k k k!1 k!1

1.3. Da tale definizione deduciamo che un metodo iterativo che converga a un valore che Osservazione non sia soluzione del nostro problema non verrà detto convergente.

Per quanto riguarda la convergenza di quantità vettoriali, è opportuno ricordare il seguente risultato.

Teorema 1.2. (analogo a (3.5)) Sia una norma in uno spazio di dimensione finita. Allora k.k V (k) (k), kx klim x = x lim x = 0 k!1 k!1 dove e è una successione di elementi di . (condizione necessaria e sufficiente per la (k)2 {x }x V V convergenza di una successione di vettori)

Definizione 1.6 monotona). (Convergenza Diremo che un metodo iterativo converge in modo monotono se (1.17) 8e < e k 0 k k 16

Criteri d’arresto Una volta stabilito che un dato metodo iterativo è convergente, prima di avviare il calcolo dobbiamo preoccuparci di inserire all’interno dell’algoritmo un certo che ne criterio d’arresto blocchi l’esecuzione nel caso in cui la soluzione a cui siamo pervenuti alla -esima iterazione sia rispetto - k sa di un certo errore massimo da noi stabilito, detto Tale controllo dell’errore effettuato dal tolleranza. criterio d’arresto può avvenire principalmente in due modi (utilizzabili anche contemporaneamente).

• Controllando ad ogni iterazione se uno o più stimatori di errore opportunamente scelti (generalmente basati sull’andamento della funzione) rispettano la tolleranza data (controllo tramite condizione).
• Stabilendo all’inizio del ciclo il numero dell’iterata dalla quale ottenere la soluzione in base all’errore che vogliamo ottenere (controllo tramite contatore)(metodo indiretto). Per fare ciò è necessario trovare uno stimatore d’errore che sia funzione di (numero dell’iterata) e porlo minore della tolleranza; k il numero di