Scarica Calcolo numerico e più Dispense in PDF di Calcolo Numerico solo su Docsity! Capitolo 1 L’insieme dei numeri macchina 1.1 Introduzione al Calcolo Numerico Il Calcolo Numerico è una disciplina che fa parte di un ampio settore della Matematica Applicata che prende il nome di Analisi Numerica. Si tratta di una materia che è al confine tra la Matematica e l’Informatica poichè cerca di risolvere i consueti problemi matematici utilizzando però una via algoritmica. In pratica i problemi vengono risolti indicando un processo che, in un numero finito di passi, fornisca una soluzione numerica e soprattutto che sia implementabile su un elaboratore. I problemi matematici che saranno affrontati nelle pagine seguenti sono problemi di base: risoluzione di sistemi lineari, approssimazione delle radici di funzioni non lineari, approssimazione di funzioni e dati sperimentali, calcolo di integrali definiti. Tali algoritmi di base molto spesso non sono altro se non un piccolo ingranaggio nella risoluzione di problemi ben più complessi. 1.2 Rappresentazione in base di un numero reale Dovendo considerare problemi in cui l’elaboratore effettua computazioni esclu- sivamente su dati di tipo numerico risulta decisivo iniziare la trattazione degli argomenti partendo dalla rappresentazione di numeri. Innanzitutto è oppor- tuno precisare che esistono due modi per rappresentare i numeri: la cosiddet- ta notazione posizionale, in cui il valore di una cifra dipende dalla posizione 1 CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 2 in cui si trova all’interno del numero, da quella notazione non posizionale, in cui ogni numero è rappresentato da uno, o da un insieme di simboli (si pensi come esempio alla numerazione usata dai Romani). La motivazione che spinge a considerare come primo problema quello della rappresentazione di numeri reali è che ovviamente si deve sapere il livello di affidabilità dei risul- tati forniti dall’elaboratore. Infatti bisogna osservare che i numeri reali sono infiniti mentre la memoria di un calcolatore ha una capacità finita che ne ren- de impossibile la rappresentazione esatta. Una seconda osservazione consiste nel fatto che un numero reale ammette molteplici modi di rappresentazione. Per esempio scrivere x = 123.47 è la rappresentazione, in forma convenzionale, dell’espressione x = 123.47 = 1× 102 + 2× 101 + 3× 100 + 4××10−1 + 7× 10−2, da cui, mettendo in evidenza 102: x = 102 × ( 1× 100 + 2× 10−1 + 3× 10−2 + 4××10−3 + 7× 10−4 ) mentre, mettendo in evidenza 103 lo stesso numero viene scritto come x = 103 × ( 1× 10−1 + 2× 10−2 + 3× 10−3 + 4××10−4 + 7× 10−5 ) deducendo che ogni numero, senza una necessaria rappresentazione conven- zionale, può essere scritto in infiniti modi. Il seguente teorema è fondamentale proprio per definire la rappresentazione dei numeri reali in una determinata base β. Teorema 1.2.1 Sia β ∈ N, β > 1, allora ogni numero reale x, x 6= 0, può essere rappresentato univocamente in base β nel seguente modo x = ± βp ∞∑ i=1 diβ −i dove p ∈ Z, e i valori di ∈ N (detti cifre), verificano le seguenti proprietà: 1. di ∈ {1, 2, 3, . . . , β − 1}; 2. d1 6= 0; 3. le cifre di non sono definivamente uguali a β − 1. CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 5 ± d1 d2 d3 d4 . . . dt segno mantissa esponente Figura 1.1: Locazione di memoria. • un campo di memoria che può assumere m+M + 1 differenti configu- razioni (e perciò può memorizzare i differenti valori p dell’esponente), • un campo che può assumere due differenti configurazioni (e perciò può memorizzare il segno + o −), è in grado di rappresentare tutti gli elementi dell’insieme F(β, t,m,M). In realtà poichè se β = 2 d1 = 1, allora determinati standard non memorizzano la prima cifra della mantissa. Il più piccolo numero positivo appartenente all’insieme F(β, t,m,M) si ottiene prendendo la più piccola mantissa (ovvero 0.1) ed il più piccolo esponente x = 0.1× β−m mentre il più grande ha tutte le cifre della mantissa uguali alla cifra più grande (ovvero β − 1) ed il massimo esponente x = 0. dd . . . dd ︸︷︷︸ t βM , d = β − 1. Consideriamo ora come esempio l’insieme F(2, 2, 2, 2), cioè i numeri binari con mantissa di due cifre ed esponente compreso tra -2 e 2. Enumeriamo gli elementi di questo insieme. Poichè il numero zero non appartiene all’insieme dei numeri macchina viene rappresentato solitamente con mantissa nulla ed CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 6 ••• • • • • • • • 0 1 2 3 • Figura 1.2: Elementi dell’insieme F(2, 2, 2, 2). esponente −m. p = −2 x = 0.10× 2−2 = 2−1 × 2−2 = 2−3 = 0.125; x = 0.11× 2−2 = (2−1 + 2−2)× 2−2 = 3/16 = 0.1875; p = −1 x = 0.10× 2−1 = 2−1 × 2−1 = 2−2 = 0.25; x = 0.11× 2−1 = (2−1 + 2−2)× 2−1 = 3/8 = 0.375; p = 0 x = 0.10× 20 = 2−1 × 20 = 2−1 = 0.5; x = 0.11× 20 = (2−1 + 2−2)× 20 = 3/4 = 0.75; p = 1 x = 0.10× 21 = 2−1 × 21 = 1; x = 0.11× 21 = (2−1 + 2−2)× 21 = 3/2 = 1.15; p = 2 x = 0.10× 22 = 2−1 × 22 = 2; x = 0.11× 22 = (2−1 + 2−2)× 22 = 3; Nella Figura 1.2 è rappresentato l’insieme dei numeri macchina positivi ap- partenenti a F(2, 2, 2, 2) (i numeri negativi sono esattamente simmetrici ri- spetto allo zero). Dalla rappresentazione dell’insieme dei numeri macchina si evincono le seguenti considerazioni: 1. L’insieme è discreto; 2. I numeri rappresentabili sono solo una piccola parte dell’insieme R; 3. La distanza tra due numeri reali consecutivi è βp−t, infatti, consideran- do per semplicità numeri positivi, sia x = +βp × (0.d1, d2, . . . , dt−1, dt) CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 7 il successivo numero macchina è y = +βp × (0.d1, d2, . . . , dt−1, d̃t) dove d̃t = dt + 1. La differenza è pertanto y − x = +βp(0. 00 . . . 00 ︸︷︷︸ t−1 1) = βp−t. Nello standard IEEE (Institute of Electric and Electronic Engineers) singola precisione una voce di memoria ha 32 bit, dei quali 1 riservato al segno, 8 all’esponente e 23 alla mantissa. Allora β = 2, t = 23,m = 127 e M = 128. Per la doppia precisione si utilizzano 64 bit, di cui 1 per il segno, 11 per l’e- sponente e 52 per la mantissa. Dunque β = 2, t = 52,m = 1023 e M = 1024. Dopo aver compreso la struttura dell’insieme F(β, t,m,M) resta da capire come, assegnato un numero reale x sia possibile rappresentarlo nell’insieme dei numeri macchina, ovvero quale elemento x̃ ∈ F(β, t,m,M) possa essergli associato in modo da commettere il più piccolo errore di rappresentazio- ne possibile. Supponiamo ora che la base β sia un numero pari. Possono presentarsi diversi casi: • Sia x = ± βp n∑ i=1 diβ −i con d1 6= 0, n ≤ t, e −m ≤ p ≤ M. Allora è evidente che x ∈ F(β, t,m,M) e pertanto verrà rappresentato esattamente su un qua- lunque elaboratore che utilizzi F(β, t,m,M) come insieme dei numeri di macchina. • Sia x = ± βp n∑ i=1 diβ −i con n ≤ t ma supponiamo che p /∈ [−m,M ]. Se p < −m allora x è più piccolo del più piccolo numero di macchina: in questo caso si dice che si è verificato un underflow (l’elaboratore interrompe la sequenza di calcoli e segnala con un messaggio l’underflow). Se p > M allora CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 10 • • a ≡ tr(x) b • x b− a = βp−t Figura 1.3: Stima dell’errore di rappresentazione nel caso di troncamento. maggiorazione per tali errori nel caso in cui x̃ sia il troncamento di x > 0. Nella Figura 1.3 a e b rappresentano i due numeri macchina tali che sia vera la relazione (1.1). È evidente ch risulta |tr(x)− x| < b− a = βp−t. Per maggiorare l’errore relativo osserviamo che |x| = +βp × 0.d1d2d3 · · · ≥ βp × 0.1 = βp−1. da cui 1 |x| ≤ β 1−p e quindi |tr(x)− x| |x| ≤ β p−t × β1−p = β1−t. (1.2) Passiamo ora alla valutazione degli errori quando x̃ = arr(x). Nella Figura 1.4 a e b rappresentano i due numeri macchina tali che sia vera la relazione (1.1). Se x > 0 si trova a sinistra del punto medio (a + b)/2 allora l’arrotondamento coincide con il valore a, se si trova nel punto medio oppure alla sua destra allora coincide con b. È evidente che il massimo errore si ottiene quando x coincide con il punto medio tra a e b risulta |arr(x)− x| ≤ 1 2 (b− a) = 1 2 βp−t. Per maggiorare l’errore relativo procediamo come nel caso del troncamento di x: |arr(x)− x| |x| ≤ 1 2 βp−t × β1−p = 1 2 β1−t. (1.3) CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 11 a+b 2 a b• • a ≡ arr(x) • x • x b ≡ arr(x) 1 2 βt−p1 2 βt−p Figura 1.4: Stima dell’errore di rappresentazione nel caso di arrotondamento. Le quantità che compaiono a destra delle maggiorazioni (1.2) e (1.3), ovvero u = β1−t oppure u = 1 2 β1−t sono dette precisione di macchina o zero macchina per il troncamento (o per l’arrotondamento, in base alla tecnica in uso). Posto εx = x̃− x x , |ε| ≤ u risulta x̃ = x(1 + εx) (1.4) che fornisce la relazione tra un numero x ∈ R e la sua rappresentazione macchina. 1.4.1 Operazioni Macchina Se x, y ∈ F(β, t,m,M) è chiaro che il risultato di un’operazione aritmetica tra x e y non è detto che sia un numero macchina, inoltre è chiaro che quanto detto per la rappresentazione dei numeri reali sia valido anche per tale risultato. Se · è una delle quattro operazioni aritmetiche di base allora affinchè il risultato sia un numero macchina deve accadere che x · y = fl(x · y). (1.5) CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 12 L’operazione definita dalla relazione (1.5) è detta operazione macchina. L’o- peraziona macchina associata a · viene indicata con ⊙ e deve soddisfare anch’essa la relazione (1.4), ovvero dev’essere: x⊙ y = (x · y)(1 + ε), |ε| < u (1.6) per ogni x, y ∈ F(β, t,m,M) tali che x ⊙ y non dia luogo ad overflow o underflow. Si può dimostrare che x⊙ y = tr(x · y) e x⊙ y = arr(x · y) soddisfano la (1.6) e dunque danno luogo ad operazioni di macchina. Le quattro operazioni cos̀ı definite danno luogo alla aritmetica di macchina o aritmetica finita. La somma algebrica macchina (addizione e sottrazione) tra due numeri x, y ∈ F(β, t,m,M) richiede le seguenti fasi: 1. Si scala la mantissa del numero con l’esponente minore in modo tale che i due addendi abbiano lo stesso esponente (ovvero quello dell’esponente maggiore); 2. Si esegue la somma tra le mantisse; 3. Si normalizza il risultato aggiustando l’esponente in modo tale che la mantissa sia un numero minore di 1. 4. Si arrotonda (o si tronca) la mantissa alle prime t cifre; Consideriamo per esempio i numeri x, y ∈ F(10, 5,m,M) x = 0.78546× 102, y = 0.61332× 10−1 e calcoliamo il numero macchina x⊕ y. 1. Scaliamo il numero y fino ad ottenere esponente 2 (quindi si deve spostare il punto decimale di 3 posizioni), y = 0.00061332× 102; 2. Sommiamo le mantisse 0.78546 + 0.00061332 = 0.78607332; 3. Questa fase non è necessaria perchè la mantissa è già minore di 1; 4. Si arrotonda alla quinta cifra decimale ottenendo x⊕ y = 0.78607× 102. CAPITOLO 1. L’INSIEME DEI NUMERI MACCHINA 15 Consideriamo la divisione tra i due numeri x = 0.12100× 105, y = 0.11000× 102 nell’insieme F(10, 5,m,M). 1. Scaliamo il dividendo di una cifra decimale 0.012100; 2. Dividiamo le mantisse 0.012100/0.11000 = 0.11000; 3. Il troncamento fornisce lo stesso numero 0.11000; 4. Si sottraggono gli esponenti ottenendo il risultato x⊘ y = 0.11000× 103. Si può dimostrare che valgono le seguenti proprietà: 1. L’insieme F(β, t,m,M) non è chiuso rispetto alle operazioni macchina; 2. L’elemento neutro per la somma non è unico: infatti consideriamo i due numeri macchia x = 0.15678× 103, y = 0.25441× 10−2, appartenenti all’insieme F(10, 5,m,M), innanzitutto si scala y y = 0.0000025441× 103, sommando le mantisse si ottiene 0.1567825441 mentre l’arrotondamen- to fornisce il risultato finale x⊕ y = 0.15678× 103 = x. 3. L’elemento neutro per il prodotto non è unico; 4. Non vale la proprietà associativa di somma e prodotto; 5. Non vale la proprietà distributiva della somma rispetto al prodotto. Capitolo 2 Equazioni non Lineari 2.1 Introduzione Le radici di un’equazione non lineare f(x) = 0 non possono, in generale, essere espresse esplicitamente e anche se ciò è possibile spesso l’espressione si presenta in forma talmente complicata da essere praticamente inutilizzabile. Di conseguenza per poter risolvere equazioni di questo tipo siamo obbligati ad utilizzare metodi numerici che sono, in generale, di tipo iterativo, cioè par- tendo da una (o in alcuni casi più) approssimazioni della radice, producono una successione x0, x1, x2, . . . , convergente alla radice. Per alcuni di questi metodi per ottenere la convergenza è sufficiente la conoscenza di un intervallo [a, b] che contiene la soluzione, altri metodi richiedono invece la conoscenza di una buona approssimazione iniziale. Talvolta è opportuno utilizzare in maniera combinata due metodi, uno del primo tipo e uno del secondo. Prima di analizzare alcuni metodi per l’approssimazione delle radici dell’e- quazione f(x) = 0 diamo la definizione di molteplicità di una radice. Definizione 2.1.1 Sia f ∈ Cr([a, b]) per un intero r > 0. Una radice α di f(x) si dice di molteplicità r se lim x→α f(x) (x− α)r = γ, γ 6= 0, γ 6= ±∞. (2.1) Se α è una radice della funzione f(x) di molteplicità r allora risulta f(α) = f ′(α) = · · · = f (r−1)(α) = 0, f (r)(α) = γ 6= 0. 16 CAPITOLO 2. EQUAZIONI NON LINEARI 17 2.2 Localizzazione delle radici Nei successivi paragrafi saranno descritti alcuni metodi numerici per il calcolo approssimato delle radici di un’equazione non lineare. Tali metodi numerici sono di tipo iterativo, ovvero consistono nel definire una successione (o più successioni), che, a partire da un’assegnata approssimazione iniziale (nota), converga alla radice α in un processo al limite. Infatti poichè non esistono tecniche generali che consentano di trovare l’espressione esplicita di α in un numero finito di operazioni, allora questa può essere calcolata in modo approssimato solo in modo iterativo. Questa perculiarità tuttavia richiede che sia nota appunto un’approssimazione iniziale o, almeno, un intervallo di appartenenza. Il problema preliminare è quello di localizzare la radice di una funzione, problema che viene affrontato in modo grafico. Per esempio considerando la funzione f(x) = sin ( log(x2 + 1) ) − e −x x2 + 1 risulta immediato verificare che il valore dell’ascissa in cui si annulla è quello in cui si intersecano i grafici delle funzioni g(x) = sin ( log(x2 + 1) ) h(x) = e−x x2 + 1 . Un modo semplice per stimare tale valore è quello di tracciare i grafici delle due funzioni, come riportato nella seguente figura in cui il grafico di h(x) è in rosso, mentre quello di g(x) è blu, e l’intervallo di variabilità di x è [0, 2.5]. CAPITOLO 2. EQUAZIONI NON LINEARI 20 dove ε è una prefissata tolleranza. La (2.2) determina anche un limite per il numero di iterate infatti: b0 − a0 2k ≤ ε ⇒ k > log2 ( b0 − a0 ε ) . Poichè bk − α ≤ bk − ak, il criterio (2.2) garantisce che α è approssimata da ck+1 con un errore assoluto minore di ε. Se 0 6∈ [a, b] si può usare come criterio di stop bk − ak min (|ak|, |bk|) ≤ ε (2.3) che garantisce che α è approssimata da ck+1 con un errore relativo minore di ε. Un ulteriore criterio di stop è fornito dal test: |f(ck)| ≤ ε. (2.4) È comunque buona norma utilizzare due criteri di stop insieme, per esempio (2.2) e (2.4) oppure (2.3) e (2.4). 2.3.1 Il metodo della falsa posizione Una variante del metodo delle bisezioni è appunto il metodo della falsa posi- zione. Partendo sempre da una funzione f(x) continua in un intervallo [a, b] tale che f(a)f(b) < 0, in questo caso si approssima la radice considerando l’intersezione della retta passante per i punti (a, f(a)) e (b.f(b)) con l’asse x. L’equazione della retta è y = f(a) + f(b)− f(a) b− a (x− a) pertanto il punto c1, sua intersezione con l’asse x, è: c1 = a− f(a) b− a f(b)− f(a) . Si testa a questo punto l’appartenenza della radice α ad uno dei due inter- valli [a, c1] e [c1, b] e si procede esattamente come nel caso del metodo delle bisezioni, ponendo [a1, b1] ≡    a1 = a, b1 = c1 se f(a)f(c1) < 0 a1 = c1, b1 = b se f(a)f(c1) > 0. CAPITOLO 2. EQUAZIONI NON LINEARI 21 Ad un generico passo k si calcola ck = ak−1 − f(ak−1) bk−1 − ak−1 f(bk−1)− f(ak−1) e si pone [ak, bk] ≡    ak = ak−1 bk = ck se f(ak−1)f(ck) < 0 ak = ck bk = bk−1 se f(ak−1)f(ck) > 0. Anche per questo metodo è possibile dimostrare la convergenza nella sola ipo- tesi di continuità della funzione f(x). Nella seguente figura è rappresentato graficamente il metodo della falsa posizione. a0 b0 c1 c2 c3 • •• I0 I1 I2 I3 function [alfa,k]=bisezione(f,a,b,tol) % % La funzione approssima la radice con il metodo di bisezione % % Parametri di input CAPITOLO 2. EQUAZIONI NON LINEARI 22 % f = funzione della quale calcolare la radice % a = estremo sinistro dell’intervallo % b = estremo destro dell’intervallo % tol = precisione fissata % % Parametri di output % alfa = approssimazione della radice % k = numero di iterazioni % if nargin==3 tol = 1e-8; % Tolleranza di default end fa = feval(f,a); fb = feval(f,b); if fa*fb>0 error(’Il metodo non e’’ applicabile’) end c = (a+b)/2; fc = feval(f,c); k = 0; while (b-a)>tol | abs(fc)>tol if fa*fc<0 b = c; fb = fc; else a = c; fa = fc; end c = (a+b)/2; fc = feval(f,c); if nargout==2 k = k+1; end end alfa = c; return CAPITOLO 2. EQUAZIONI NON LINEARI 25 x0x1 x2x4x3 α • •• ••• y = g(x) Figura 2.1: Interpretazione geometrica del processo xk+1 = g(xk), se −1 < g′(α) ≤ 0. con ξ ∈ [α− ρ, α + ρ]. Poichè |g′(ξ)| < 1 si ha |α− β| < |α− β| e ciò è assurdo. ✷ Nelle figure 2.2 e 2.1 è rappresentata l’interpretazione geometrica di un metodo di iterazione funzionale in ipotesi di convergenza. Definizione 2.4.1 Un metodo iterativo del tipo (2.5) si dice localmente con- vergente ad una soluzione α del problema f(x) = 0 se esiste un intervallo [a, b] contenente α tale che, per ogni x0 ∈ [a, b], la successione generata da (2.5) converge a α. Come abbiamo già visto nel caso del metodo delle bisezioni anche per metodi di iterazione funzionale è necessario definire dei criteri di arresto per il calcolo delle iterazioni. Teoricamente, una volta stabilita la precisione voluta, ε, si dovrebbe arrestare il processo iterativo quando l’errore al passo k ek = |α− xk| CAPITOLO 2. EQUAZIONI NON LINEARI 26 x0 x1 x2 x3 x4 α •• • • • • y = g(x) Figura 2.2: Interpretazione geometrica del processo xk+1 = g(xk), se 0 ≤ g′(α) < 1. risulta minore della tolleranza prefissata ε. In pratica l’errore non può essere noto quindi è necessario utilizzare qualche stima. Per esempio si potrebbe considerare la differenza tra due iterate consecutive e fermare il calcolo degli elementi della successione quando |xk+1 − xk| ≤ ε, oppure |xk+1 − xk| min(|xk+1|, |xk|) ≤ ε |xk+1|, |xk| 6= 0 se i valori hanno un ordine di grandezza particolarmente elevato. Una stima alternativa valuta il residuo della funzione rispetto al valore in α, cioè |f(xk)| ≤ ε. 2.4.1 Ordine di Convergenza Per confrontare differenti metodi iterativi che approssimano la stessa radice α di f(x) = 0, si può considerare la velocità con cui tali successioni convergono CAPITOLO 2. EQUAZIONI NON LINEARI 27 verso α. Lo studio della velocità di convergenza passa attraverso il concetto di ordine del metodo. Definizione 2.4.2 Sia {xk}∞k=0 una successione convergente ad α e tale che xk 6= α, per ogni k. Se esiste un numero reale p ≥ 1 tale che lim k→+∞ |xk+1 − α| |xk − α|p = γ con    0 < γ ≤ 1 se p = 1 γ > 0 se p > 1 (2.8) allora si dice che la successione ha ordine di convergenza p. La costante γ prende il nome di costante asintotica di convergenza. In particolare se p = 1 e 0 < γ < 1 allora la convergenza si dice lineare, mentre se p > 1 allora la convergenza si dice genericamente superlineare, per esempio se p = 2 la convergenza si dice quadratica, se p = 3 cubica e cos̀ı via. Osservazione. La relazione (2.8) implica che esiste una costante positiva β (β ≃ γ) tale che, per k sufficientemente grande: |xk+1 − α| ≤ β|xk − α|p (2.9) ed anche |xk+1 − α| |α| ≤ β|α| p−1 ∣ ∣ ∣ ∣ xk − α α ∣ ∣ ∣ ∣ p . (2.10) Le (2.9) e (2.10) indicano che la riduzione di errore (assoluto o relativo) ad ogni passo è tanto maggiore quanto più alto è l’ordine di convergenza e, a parità di ordine, quanto più piccola è la costante asintotica di convergenza. In generale l’ordine di convergenza è un numero reale maggiore o uguale a 1. Tuttavia per i metodi di iterazione funzionale di tipo (2.5) è un numero intero per il quale vale il seguente teorema. Teorema 2.4.3 Sia {xk}∞k=0 una successione generata dallo schema (2.5) convergente ad α, punto fisso di g(x), funzione sufficientemente derivabile in un intorno di α. La successione ha ordine di convergenza p ≥ 1 se e solo se g′(α) = g′′(α) = · · · = g(p−1)(α) = 0, g(p)(α) 6= 0. (2.11) CAPITOLO 2. EQUAZIONI NON LINEARI 30 Per la convergenza el’ordine del metodo di Newton-Raphson vale il seguente teorema. Teorema 2.4.4 Sia f ∈ C3([a, b]), tale che f ′(x) 6= 0, per x ∈ [a, b], do- ve [a, b] è un opportuno intervallo contenente α, allora valgono le seguenti proposizioni: 1. esiste un intervallo [α − ρ, α + ρ], tale che, scelto x0 appartenente a tale intervallo, la successione definita dal metodo di Newton-Raphson è convergente ad α; 2. la convergenza è di ordine p ≥ 2. Dimostrazione. Per valutare la convergenza del metodo calcoliamo la derivata prima della funzione iteratrice: g′(x) = 1− [f ′(x)]2 − f(x)f ′′(x) [f ′(x)]2 = f(x)f ′′(x) [f ′(x)]2 . Poichè f ′(α) 6= 0 risulta: g′(α) = f(α)f ′′(α) [f ′(α)]2 = 0 quindi, fissato un numero positivo κ < 1, esiste ρ > 0 tale che per ogni x ∈ [α − ρ, α + ρ] si ha |g′(x)| < κ e quindi vale il teorema di convergenza 2.4.2. Per dimostrare la seconda parte del teorema si deve calcolare la derivata seconda di g(x): g′′(x) = [f ′(x)f ′′(x) + f(x)f ′′′(x)][f ′(x)]2 − 2f(x)f ′(x)[f ′′(x)]2 [f ′(x)]4 . Calcolando la derivata seconda in x = α risulta g′′(α) = f ′′(α) f ′(α) (2.14) ne segue che se f ′′(α) 6= 0 allora anche g′′(α) 6= 0 e quindi, applicando il Teorema 2.4.3, l’ordine p = 2. Se invece f ′′(α) = 0 allora l’ordine è almeno CAPITOLO 2. EQUAZIONI NON LINEARI 31 pari a 3. Dalla relazione 2.14 segue inoltre che la costante asintotica di convergenza vale γ = 1 2 ∣ ∣ ∣ ∣ f ′′(α) f ′(α) ∣ ∣ ∣ ∣ . ✷ Il Teorema 2.4.4 vale nell’ipotesi in cui f ′(α) 6= 0, cioè se α è una radice semplice di f(x). Se invece la radice α ha molteplicità r > 1 l’ordine di convergenza del metodo non è più 2. In questo caso infatti si può porre f(x) = q(x)(x− α)r, q(α) 6= 0, quindi riscrivendo la funzione iteratrice del metodo di Newton-Raphson ri- sulta g(x) = x− q(x)(x− α) rq(x) + q′(x)(x− α) , da cui, dopo una serie di calcoli, risulta g′(α) = 1− 1 r . (2.15) Pertanto, poichè r > 1 risulta |g′(x)| < 1 e quindi per il Teorema 2.4.2 il metodo è ancora convergente ma, applicando il Teorema 2.4.3 l’ordine di convergenza è 1. Se si conosce la molteplicità della radice si può modificare il metodo di Newton-Raphson ottenendo uno schema numerico con ordine 2. Ponendo xk+1 = xk − r f(xk) f ′(xk) k = 0, 1, 2, . . . si definisce un metodo con la seguente funzione iteratrice g(x) = x− r f(x) f ′(x) da cui segue, tenendo conto della (2.15), che g′(α) = 0. Riportiamo nel seguito l’implementazione MatLab del metodo di Newton- Raphson. CAPITOLO 2. EQUAZIONI NON LINEARI 32 function [alfa,k]=newton(f,f1,x0,tol,Nmax) % % La funzione calcolo un’approssimazione % della radice con il metodo di Newton-Raphson % % Parametri di input % f = funzione della quale calcolare la radice % f1 = derivata prima della funzione f % x0 = approssimazione iniziale della radice % tol = precisione fissata % Nmax = numero massimo di iterazioni fissate % % Parametri di output % alfa = approssimazione della radice % k = numero di iterazioni % if nargin==3 tol=1e-8; Nmax=1000; end k=0; x1=x0-feval(f,x0)/feval(f1,x0); fx1 = feval(f,x1); while abs(x1-x0)>tol | abs(fx1)>tol x0 = x1; x1 = x0-feval(f,x0)/feval(f1,x0); fx1 = feval(f,x1); k=k+1; if k>Nmax disp(’Il metodo non converge’); alfa = inf; break end end alfa=x1; return Esempio 2.4.1 Approssimare il numero α = m √ c con m ∈ R, m ≥ 2, c > 0. CAPITOLO 2. EQUAZIONI NON LINEARI 35 2.4.4 Il Metodo della Secante Il metodo della secante è definito dalla relazione xk+1 = xk − f(xk) xk − c f(xk)− f(c) dove c ∈ [a, b]. Il significato geometrico di tale metodo è il seguente: ad un generico passo k si considera la retta congiungente i punti di coordinate (xk, f(xk)) e (c, f(c)) e si pone xk+1 pari all’ascissa del punto di intersezione di tale retta con l’asse x. Dalla formula si evince che la funzione iteratrice del metodo è g(x) = x− f(x) x− c f(x)− f(c) . Il metodo è rappresentato graficamente nella seguente figura. c (c, f(c)) • • • • x0 x1 x2 In base alla teoria vista nei paragrafi precedenti il metodo ha ordine di conver- genza 1 se g′(α) 6= 0. Può avere ordine di convergenza almeno 1 se g′(α) = 0. Tale eventualità si verifica se la tangente alla curva in α ha lo stesso coeffi- ciente angolare della retta congiungente i punti (α, 0) e (c, f(c)). Poichè il metodo delle secanti ha lo svantaggio di avere, solitamente, conver- genza lineare mentre il metodo di Newton-Raphson, pur avendo convergenza quadratica, ha lo svantaggio di richiedere, ad ogni passo, due valutazioni di funzioni: f(xk) ed f ′(xk), quindi se il costo computazionale di f ′(xk) è mol- to più elevato rispetto a quello di f(xk) può essere più conveniente l’uso di metodi che necessitano solo del calcolo del valore della funzione f(x). Capitolo 3 Metodi numerici per sistemi lineari 3.1 Introduzione Siano assegnati una matrice non singolare A ∈ Rn×n ed un vettore b ∈ Rn. Risolvere un sistema lineare avente A come matrice dei coefficienti e b come vettore dei termini noti significa trovare un vettore x ∈ Rn tale che Ax = b. (3.1) Esplicitare la relazione (3.1) significa imporre le uguaglianze tra le compo- nenti dei vettori a primo e secondo membro: a11x1+ a12x2+ · · ·+ a1nxn = b1 a21x1+ a22x2+ · · ·+ a2nxn = b2 ... ... an1x1+ an2x2+ · · ·+ annxn = bn. (3.2) Le (3.2) definiscono un sistema di n equazioni algebriche lineari nelle n inco- gnite x1, x2, . . . , xn. Il vettore x viene detto vettore soluzione. Prima di af- frontare il problema della risoluzione numerica di sistemi lineari richiamiamo alcuni importanti concetti di algebra lineare. Definizione 3.1.1 Se A ∈ Rn×n è una matrice di ordine 1, si definisce determinante di A il numero detA = a11. 36 CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 37 Se la matrice A è quadrata di ordine n allora fissata una qualsiasi riga (colon- na) di A, diciamo la i-esima (j-esima) allora applicando la cosiddetta regola di Laplace il determinante di A è: detA = n∑ j=1 aij(−1)i+j detAij dove Aij è la matrice che si ottiene da A cancellando la i-esima riga e la j-esima colonna. Il determinante è pure uguale a detA = n∑ i=1 aij(−1)i+j detAij, cioè il determinante è indipendente dall’indice di riga (o di colonna) fissato. Se A è la matrice di ordine 2 A = [ a11 a12 a21 a22 ] . allora detA = a11a22 − a21a12. Il determinante ha le seguenti proprietà: 1. Se A è una matrice triangolare o diagonale allora detA = n∏ i=1 aii; 2. det I = 1; 3. detAT = detA; 4. detAB = detA detB (Regola di Binet); 5. se α ∈ R allora detαA = αn detA; 6. detA = 0 se una riga (o una colonna) è nulla, oppure una riga (o una colonna) è proporzionale ad un’altra riga (o colonna) oppure è combinazione lineare di due (o più) righe (o colonne) di A. 7. Se A è una matrice triangolare a blocchi A = [ B C O D ] CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 40 Anche per il seguente sistema il vettore soluzione è calcolabile in modo analogo. a11x1 = b1 a21x1 +a22x2 = b2 ... ... . . . ... ai1x1 +ai2x2 . . . +aiixi = bi ... ... . . . ... an1x1 +an2x2 . . . +anixi . . . +annxn = bn (3.7) In questo caso la matrice dei coefficienti è triangolare inferiore e la soluzione viene calcolata con il metodo di sostituzione in avanti:    x1 = b1 a11 xi = bi − i−1∑ j=1 aijxj aii i = 2, . . . , n. Concludiamo questo paragrafo facendo alcune considerazioni sul costo com- putazionale dei metodi di sostituzione. Per costo computazionale di un al- goritmo si intende il numero di operazioni che esso richiede per fornire la soluzione di un determinato problema. Nel caso di algoritmi numerici le ope- razioni che si contano sono quelle aritmetiche su dati reali. Considerano per esempio il metodo di sostituzione in avanti. Per calcolare x1 è necessaria una sola operazione (una divisione), per calcolare x2 le operazioni sono tre (un prodotto, una somma algebrica e una divisione), mentre il generico xi richie- de 2i− 1 operazioni (i− 1 prodotti, i− 1 somme algebriche e una divisione), indicato con C(n) il numero totale di operazioni necessarie è: C(n) = n∑ i=1 (2i− 1) = 2 n∑ i=1 i− n∑ i=1 1 = 2 n(n+ 1) 2 − n = n2, sfruttando la proprietà che n∑ i=1 i = n(n+ 1) 2 . CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 41 Il costo computazionale viene sempre valutato in funzione di un determinato parametro (il numero assoluto in sè non avrebbe alcun significato) che, in questo caso è la dimensione del sistema. In questo modo è possibile prevedere il tempo necessario per calcolare la soluzione del problema. 3.3 Metodo di Eliminazione di Gauss L’idea di base del metodo di Gauss è appunto quella di operare delle oppor- tune trasformazioni sul sistema originale Ax = b, che non costino eccessiva- mente, in modo da ottenere un sistema equivalente1 avente come matrice dei coefficienti una matrice triangolare superiore. Supponiamo di dover risolvere il sistema: 2x1 +x2 +x3 = −1 −6x1 −4x2 −5x3 +x4 = 1 −4x1 −6x2 −3x3 −x4 = 2 2x1 −3x2 +7x3 −3x4 = 0. Il vettore soluzione di un sistema lineare non cambia se ad un’equazione viene sommata la combinazione lineare di un’altra equazione del sistema. L’idea alla base del metodo di Gauss è quella di ottenere un sistema linea- re con matrice dei coefficienti triangolare superiore effettuando opportune combinazioni lineari tra le equazioni. Poniamo A(1) =     2 1 1 0 −6 −4 −5 1 −4 −6 −3 −1 2 −3 7 −3     , b(1) =     −1 1 2 0     rispettivamente la matrice dei coefficienti e il vettore dei termini noti del si- stema di partenza. Calcoliamo un sistema lineare equivalente a quello iniziale ma che abbia gli elementi sottodiagonali della prima colonna uguali a zero. Azzeriamo ora l’elemento a21(1). Lasciamo inalterata la prima equazione. 1Due sistemi si dicono equivalenti se ammettono lo stesso insieme di soluzioni, quindi nel nostro caso la stessa soluzione. Osserviamo che se x∗ è un vettore tale che Ax∗ = b e B è una matrice non singolare allora BAx∗ = Bb; viceversa se BAx∗ = Bb e B è non singolare allora B−1BAx∗ = B−1Bb e quindi Ax∗ = b. Dunque se B è non singolare i sistemi Ax = b e BAx = Bb sono equivalenti. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 42 Poniamo l21 = − a21 a11 = −−6 2 = 3 e moltiplichiamo la prima equazione per l21 ottenendo: 6x1 + 3x2 + 3x3 = −3. La nuova seconda equazione sarà la somma tra la seconda equazione e la prima moltiplicata per l21: −6x1 −4x2 −5x3 +x4 = 1 6x1 +3x2 +3x3 = −3 −x2 −2x3 +x4 = −2 [Nuova seconda equazione]. Prcediamo nello stesso modo per azzerare gli altri elementi della prima co- lonna. Poniamo l31 = − a (1) 31 a (1) 11 = −−4 2 = 2 e moltiplichiamo la prima equazione per l31 ottenendo: 4x1 + 2x2 + 2x3 = −2. La nuova terza equazione sarà la somma tra la terza equazione e la prima moltiplicata per l31: −4x1 −6x2 −3x3 −x4 = 2 4x1 +2x2 +2x3 = −2 −4x2 −x3 −x4 = 0 [Nuova terza equazione]. Poniamo ora l41 = − a (1) 41 a (1) 11 = −2 2 = −1 e moltiplichiamo la prima equazione per l41 ottenendo: −2x1 − x2 − x3 = 1. La nuova quarta equazione sarà la somma tra la quarta equazione e la prima moltiplicata per l41: 2x1 −3x2 +7x3 −3x4 = 0 −2x1 −x2 −x3 = 1 −4x2 +6x3 −3x4 = 1 [Nuova quarta equazione]. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 45 Abbiamo ottenuto un sistema triangolare superiore: 2x1 +x2 +x3 = −1 −x2 −2x3 +x4 = 4 7x3 −5x4 = 8 3x4 = −7. La matrice dei coefficienti e il vettore dei termini noti sono diventati: A(4) =     2 1 1 0 0 −1 −2 1 0 0 7 −5 0 0 0 3     , b(4) =     −1 4 8 −7     . Vediamo come ciò sia possibile. Riconsideriamo il sistema di equazioni nella sua forma scalare (3.2): n∑ j=1 aijxj = bi, i = 1, . . . , n. (3.8) Per motivi che risulteranno chiari tra poco poniamo a (1) ij = aij e b (1) i = bi. Isoliamo in ogni equazione la componente x1. Abbiamo: a (1) 11 x1 + n∑ j=2 a (1) 1j xj = b (1) 1 (3.9) a (1) i1 x1 + n∑ j=2 a (1) ij xj = b (1) i , i = 2, . . . , n. (3.10) Dividendo l’equazione (3.9) per a (1) 11 e moltiplicandola rispettivamente per −a(1)21 ,−a (1) 31 , . . . ,−a (1) n1 si ottengono n− 1 nuove equazioni: a (1) i1 x1 + n∑ j=2 a (1) i1 a (1) 11 a (1) 1j xj = a (1) i1 a (1) 11 b (1) i , i = 2, . . . , n. (3.11) Sommando da (3.10) le equazioni (3.11) si ottiene n∑ j=2 ( a (1) ij − a (1) i1 a (1) 11 a (1) 1j ) xj = b (1) i − a (1) i1 a (1) 11 b (1) 1 , i = 2, . . . , n. (3.12) CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 46 L’equazione (3.9) insieme alle (3.12) formano un nuovo sistema di equazioni, equivalente a quello originario, che possiamo scrivere cos̀ı:    a (1) 11 x1 + n∑ j=2 a (1) 1j xj = b (1) 1 n∑ j=2 a (2) ij xj = b (2) i i = 2, . . . , n (3.13) dove    a (2) ij = a (1) ij − a (1) i1 a (1) 11 a (1) 1j i, j = 2, . . . , n b (2) i = b (1) i − a (1) i1 a (1) 11 b (1) 1 i = 2, . . . , n. (3.14) Osserviamo che la matrice dei coefficienti del sistema (3.13) è la seguente A(2) =      a (1) 11 a (1) 12 . . . a (1) 1n 0 a (2) 22 . . . a (2) 2n ... ... ... 0 a (2) n2 . . . a (2) nn      . Ora a partire dal sistema di equazioni n∑ j=2 a (2) ij xj = b (2) i i = 2, . . . , n, ripetiamo i passi fatti precedentemente: a (2) 22 x2 + n∑ j=3 a (2) 2j xj = b (2) 2 (3.15) a (2) i2 x2 + n∑ j=3 a (2) ij xj = b (2) i , i = 3, . . . , n. (3.16) CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 47 Moltiplicando l’equazione (3.15) per −a (2) i2 a (2) 22 , per i = 3, . . . , n, si ottiene a (2) i2 x2 + n∑ j=3 ( −a (2) i2 a (2) 22 a (2) 2j ) xj = a (2) i2 a (2) 22 b (2) 2 , i = 3, . . . , n. (3.17) Sommando a questo punto le equazioni (3.17) alle (3.16) si ottiene: n∑ j=3 ( a (2) ij − a (2) i2 a (2) 22 a (2) 2j ) xj = b (2) i − a (2) i2 a (2) 22 b (2) 2 , i = 3, . . . , n (3.18) ovvero scritta in forma più compatta: n∑ j=3 a (3) ij xj = b (3) i i = 3, . . . , n dove    a (3) ij = a (2) ij − a (2) i2 a (2) 22 a (2) 2j i, j = 3, . . . , n b (3) i = b (2) i − a (2) i2 a (2) 22 b (2) 2 i = 3, . . . , n. Abbiamo il nuovo sistema equivalente:    n∑ j=1 a (1) 1j xj = b (1) 1 n∑ j=2 a (2) 2j xj = b (2) 2 n∑ j=3 a (3) ij xj = b (3) i i = 3, . . . , n. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 50 Abbiamo ottenuto la seguente matrice al passo 2: A(2) =     3 3 5 0 0 −1 1 −1 0 2 0 4 0 2 −5/3 4     . Calcoliamo i due moltiplicatori l3,2 = 2, l4,2 = 2. Calcoliamo la terza riga: [3a riga di A(2) + ] 0 2 0 4 + [(2)× 2a riga di A(2)] 0 −2 2 −2 = [3a riga di A(3)] 0 0 2 2 La quarta riga è [4a riga di A(2) + ] 0 2 −5/3 4 + [(2)× 2a riga di A(2)] 0 −2 2 −2 = [4a riga di A(3)] 0 0 1/3 2 Abbiamo ottenuto la seguente matrice al passo 3: A(3) =     3 3 5 0 0 −1 1 −1 0 0 2 2 0 0 1/3 2     . Calcoliamo l’unico moltiplicatore del terzo passo: l4,3 = − 1 6 . La quarta riga è [4a riga di A(3) + ] 0 0 1/3 2 + [(−1/6)× 3a riga di A(3)] 0 0 −1/3 −1/3 = [4a riga di A(4)] 0 0 0 5/3 CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 51 La matrice triagolarizzata è A(4) =     3 3 5 0 0 −1 1 −1 0 0 2 2 0 0 0 5/3     . Il determinante della matrice è uguale al prodotto degli elementi diagonali della matrice triangolare, ovvero detA = −10. Esempio 3.3.2 Calcolare l’inversa della matrice A =     2 1 0 1 1 1 2 0 −1 0 3 1 1 1 2 2     utlizzando il metodo di eliminazione di Gauss. L’inversa di A è la matrice X tale che AX = I ovvero, detta xi la i−esima colonna di X, questo è soluzione del sistema lineare Axi = ei (3.21) dove ei è l’i−esimo versore della base canonica di Rn. Posto i = 1 risolvendo il sistema Ax1 = e1,     2 1 0 1 1 1 2 0 −1 0 3 1 1 1 2 2         x1 x2 x3 x4     =     1 0 0 0     si ottengono gli elementi della prima colonna di A−1. Posto A(1) = A gli elementi della matrice al passo 2 sono calcolati applicando le formule a (2) ij = a (1) ij − a (1) i1 a (1) 11 a (1) 1j , i, j = 2, 3, 4. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 52 Tralasciando il dettaglio delle operazioni risulta A(2) =     2 1 0 1 0 1/2 2 −1/2 0 1/2 3 3/2 0 1/2 2 3/2     , e (2) 1 =     1 −1/2 1/2 −1/2     Applicando le formula a (3) ij = a (2) ij − a (2) i2 a (2) 22 a (2) 2j , i, j = 3, 4. si ottiene il sistema al terzo passo A(3) =     2 1 0 1 0 1/2 2 −1/2 0 0 1 2 0 0 0 2     , e (3) 1 =     1 −1/2 1 0     . In questo caso non è necessario applicare l’ultimo passo del metodo in quanto la matrice è già triangolare superiore e pertanto si può risolvere il sistema triangolare superiore ottenendo: x4 = 0, x3 = 1, x2 = −5, x1 = 3. Cambiando i termini noti del sistema (3.21), ponendo i = 2, 3, 4 si ottengono le altre tre colonne della matrice inversa. 3.3.1 Costo Computazionale del Metodo di Elimina- zione di Gauss Cerchiamo ora di determinare il costo computazionale (cioè il numero di operazioni aritmetiche) richiesto dal metodo di eliminazione di Gauss per risolvere un sistema lineare di ordine n. Dalle relazioni b (k+1) i = b (k) i − a (k) ik a (k) kk b (k) k , i = k + 1, . . . , n, a (k+1) ij = a (k) ij − a (k) ik a (k) kk a (k) kj , i, j = k + 1, . . . , n CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 55 k i Figura 3.1: Strategia di pivoting parziale. perchè la matrice A (k) 22 ha una colonna nulla. Poichè tutte la matrici A (k) han- no lo stesso determinante di A, dovrebbe essere detA = 0 e questo contrasta con l’ipotesi fatta. Possiamo concludere che se a (k) kk = 0 e detA 6= 0 deve necessariamente esistere un elemento a (k) ik 6= 0, con i ∈ {k + 1, k + 2, . . . , n}. Per evitare che un elemento pivotale possa essere uguale a zero si applica una delle cosiddette strategie di pivoting. La strategia di Pivoting parziale prevede che prima di fare ciò si ricerchi l’elemento di massimo modulo tra gli elementi a (k) kk , a (k) k+1,k, . . . , a (k) nk e si scambi la riga in cui si trova questo elemento con la k-esima qualora esso sia diverso da a (k) kk . In altri termini il pivoting parziale richiede le seguenti operazioni: 1. determinare l’elemento a (k) rk tale che |a(k)rk | = max k≤i≤n |a(k)ik |; 2. effettuare lo scambio tra la r-esima e la k-esima riga e tra la s-esima e la k-esima colonna. in alternativa si può adottare la strategia di pivoting totale che è la seguente: 1. determinare gli indici r, s tali che |a(k)rs | = max k≤i,j≤n |a(k)ij |; CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 56 k i k j Figura 3.2: strategia di pivoting totale. 2. effettuare lo scambio tra la r-esima e la k-esima riga e tra la s-esima e la k-esima colonna. La strategia di pivoting totale è senz’altro migliore perchè garantisce mag- giormente che un elemento pivotale non sia un numero piccolo (in questa eventualità potrebbe accadere che un moltiplicatore sia un numero molto grande) ma richiede che tutti gli eventuali scambi tra le colonne della matri- ce siano memorizzati. Infatti scambiare due colonne significa scambiare due incognite del vettore soluzione pertanto dopo la risoluzione del sistema trian- golare per ottenere il vettore soluzione del sistema di partenza è opportuno permutare le componenti che sono state scambiate. Esempio 3.3.3 Risolvere il sistema lineare Ax = b dove A =     1 2 −1 0 2 −1 −1 1 3 0 −1 1 1 −3 1 1     , b =     2 1 4 2     utlizzando il metodo di eliminazione di Gauss con strategia di pivoting par- ziale. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 57 Posto A(1) = A, osserviamo che l’elemento pivotale della prima colonna si trova sulla terza riga allora scambiamo per equazioni 1 e 3: A(1) =     3 0 −1 1 2 −1 −1 1 1 2 −1 0 1 −3 1 1     , b(1) =     4 1 2 2     calcoliamo i tre moltiplicatori l2,1 = − 2 3 , l3,1 = − 1 3 , l4,1 = − 1 3 . Calcoliamo la seconda riga: [2a riga di A(1) + ] 2 −1 −1 1 1 + [(−2/3)× 1a riga di A(1)] −2 0 2/3 −2/3 −8/3 = [2a riga di A(2)] 0 −1 −1/3 1/3 −5/3 La terza riga è la seguente: [3a riga di A(1) + ] 1 2 −1 0 2 + [(−1/3)× 1a riga di A(1)] −1 0 1/3 −1/3 −4/3 = [3a riga di A(2)] 0 2 −2/3 −1/3 2/3 mentre la quarta riga è [4a riga di A(1) + ] 1 −3 1 1 2 + [(−1/3)× 1a riga di A(1)] −1 0 1/3 −1/3 −4/3 = [4a riga di A(2)] 0 −3 4/3 2/3 2/3 Abbiamo ottenuto la matrice ed il vettore al passo 2: A(2) =     3 0 −1 1 0 −1 −1/3 1/3 0 2 −2/3 −1/3 0 −3 4/3 2/3     , b(2) =     4 −5/3 2/3 2/3     . L’elemento pivotale della seconda colonna si trova sulla quarta riga quindi scambiamo le equazioni 2 e 4: A(2) =     3 0 −1 1 0 −3 4/3 2/3 0 2 −2/3 −1/3 0 −1 −1/3 1/3     , b(2) =     4 2/3 2/3 −5/3     . CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 60 end x(n) = b(n)/A(n,n); for i=n-1:-1:1 x(i) = (b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i); end return function x=Gauss_pp(A,b) % % Metodo di Gauss con pivot parziale % % Parametri di input: % A = Matrice dei coefficienti del sistema % b = Vettore dei termini noti del sistema % % Parametri di input: % x = Vettore soluzione del sistema lineare % n = length(b); x = zeros(n,1); for k=1:n-1 [a,i] = max(abs(A(k:n,k))); i = i+k-1; if i~=k A([i k],:) = A([k i],:); b([i k]) = b([k i]); end for i=k+1:n A(i,k) = A(i,k)/A(k,k); b(i) = b(i)-A(i,k)*b(k); for j=k+1:n A(i,j) = A(i,j)-A(i,k)*A(k,j); end end end x(n) = b(n)/A(n,n); for i=n-1:-1:1 x(i) = (b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i); CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 61 end return function x=Gauss_pt(A,b) % % Metodo di Gauss con pivot totale % % Parametri di input: % A = Matrice dei coefficienti del sistema % b = Vettore dei termini noti del sistema % % Parametri di input: % x = Vettore soluzione del sistema lineare % n = length(b); x = zeros(n,1); x1 = x; indice = [1:n]; for k=1:n-1 [a,riga] = max(abs(A(k:n,k:n))); [mass,col] = max(a); j = col+k-1; i = riga(col)+k-1; if i~=k A([i k],:) = A([k i],:); b([i k]) = b([k i]); end if j~=k A(:,[j k]) = A(:,[k j]); indice([j k]) = indice([k j]); end for i=k+1:n A(i,k) = A(i,k)/A(k,k); b(i) = b(i)-A(i,k)*b(k); for j=k+1:n A(i,j) = A(i,j)-A(i,k)*A(k,j); end end CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 62 end % % Risoluzione del sistema triangolare superiore % x1(n) = b(n)/A(n,n); for i=n-1:-1:1 x1(i) = (b(i)-A(i,i+1:n)*x1(i+1:n))/A(i,i); end % % Ripermutazione del vettore % for i=1:n x(indice(i))=x1(i); end return 3.3.3 La Fattorizzazione LU Supponiamo di dover risolvere un problema che richieda, ad un determinato passo, la risoluzione del sistema lineare Ax = b e di utilizzare il metodo di Gauss. La matrice viene resa triangolare superiore e viene risolto il sistema triangolare A(n)x = b(n). (3.22) Ipotizziamo che, nell’ambito dello stesso problema, dopo un certo tempo sia necessario risolvere il sistema Ax = c i cui la matrice dei coefficienti è la stessa mentre è cambiato il termine noto. Appare chiaro che non è possibile sfruttare i calcoli gia fatti in quanto il calcolo del vettore dei termini noti al passo n dipende dalle matrici ai passi precedenti all’ultimo, quindi la conoscenza della matrice A(n) è del tutto inutile. È necessario pertanto applicare nuovamente il metodo di Gauss e risolvere il sistema triangolare A(n)x = c(n). (3.23) L’algoritmo che sarà descritto in questo paragrafo consentirà di evitare l’e- ventualità di dover rifare tutti i calcoli (o una parte di questi). La Fattoriz- zazione LU di una matrice stabilisce, sotto determinate ipotesi, l’esistenza CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 65 1 2 3 4 5 6 7 1 3 5 7 2 4 6 Tecnica di Crout Tecnica di Doolittle Ogni schema rappresenta in modo schematico una matrice la cui parte trian- golare superiore indica la matrice U mentre quella triangolare inferiore la matrice L mentre i numeri indicano l’ordine con cui gli elementi saranno calcolati. Per esempio applicando la tecnica di Crout si segue il seguente ordine: • 1o Passo: Calcolo della prima riga di U ; • 2o Passo: Calcolo della seconda riga di L; • 3o Passo: Calcolo della seconda riga di U ; • 4o Passo: Calcolo della terza riga di L; • 5o Passo: Calcolo della terza riga di U ; • 6o Passo: Calcolo della quarta riga di L; • 7o Passo: Calcolo della quarta riga di U ; e cos̀ı via procedendo per righe in modo alternato. Nel caso della tecnica di Doolittle si seguono i seguenti passi: • 1o Passo: Calcolo della prima riga di U ; • 2o Passo: Calcolo della prima colonna di L; • 3o Passo: Calcolo della seconda riga di U ; CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 66 • 4o Passo: Calcolo della seconda colonna di L; • 5o Passo: Calcolo della terza riga di U ; • 6o Passo: Calcolo della terza colonna di L; • 7o Passo: Calcolo della quarta riga di U . La fattorizzazione LU è un metodo sostanzialmente equivalente al metodo di Gauss, infatti la matrice U che viene calcolata coincide con la matrice A(n). Lo svantaggio del metodo di fattorizzazione diretto risiede essenzialmente nella maggiore difficoltà, rispetto al metodo di Gauss, di poter programmare una strategia di pivot. Infatti se un elemento diagonale della matrice U è uguale a zero non è possibile applicare l’algoritmo. function [L,U]=crout(A); % % La funzione calcola la fattorizzazione LU della % matrice A applicando la tecnica di Crout % % L = matrice triang. inferiore con elementi diagonali % uguali a 1 % U = matrice triangolare superiore % [m n] = size(A); U = zeros(n); L = eye(n); U(1,:) = A(1,:); for i=2:n for j=1:i-1 L(i,j) = (A(i,j) - L(i,1:j-1)*U(1:j-1,j))/U(j,j); end for j=i:n U(i,j) = A(i,j) - L(i,1:i-1)*U(1:i-1,j); end end return function [L,U]=doolittle(A); CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 67 % % La funzione calcola la fattorizzazione LU della % matrice A applicando la tecnica di Doolittle % % L = matrice triang. inferiore con elementi diagonali % uguali a 1 % U = matrice triangolare superiore % [m n] = size(A); L = eye(n); U = zeros(n); U(1,:) = A(1,:); for i=1:n-1 for riga=i+1:n L(riga,i)=(A(riga,i)-L(riga,1:i-1)*U(1:i-1,i))/U(i,i); end for col=i+1:n U(i+1,col) = A(i+1,col)-L(i+1,1:i)*U(1:i,col); end end return Esempio 3.3.4 Applicare la tecnica di Doolittle per calcolare la fattorizza- zione LU della matrice A =     1 −1 3 −4 2 −3 9 −9 3 1 −1 −10 1 2 −4 −1     . Gli elementi della prima riga di U vanno calcolati utilizzando la formula (3.26) con i = 1: u1j = a1j − 0∑ k=1 l1kukj = a1j, j = 1, 2, 3, 4. Quindi U =     1 −1 3 −4 0 u22 u23 u24 0 0 u33 u34 0 0 0 u44     . CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 70 e U =     1 −1 3 −4 0 −1 3 −1 0 0 2 −2 0 0 0 2     . Esercizio 3.3.1 Risolvere il problema descritto nell’esempio 3.3.2 calcolan- do la fattorizzazione LU della matrice A. 3.4 Condizionamento di sistemi lineari Nel Capitolo 1 è stato introdotto il concetto di rappresentazione in base ed è stata motivata la sostanziale inaffidabilità dei risultati dovuti ad elaborazio- ni numeriche, a causa dell’artimetica finita dell’elaboratore. Appare chiaro come la bassa precisione nel calcolo potrebbe fornire dei risultati numeri- ci molto lontani da quelli reali. In alcuni casi tale proprietà è insita nel problema. Consideriamo il sistema lineare Ax = b (3.29) dove A ∈ Rn×n è la cosiddetta matrice di Hilbert, i cui elementi sono aij = 1 i+ j − 1 , i, j = 1, . . . , n mentre il vettore b è scelto in modo tale che il vettore soluzione abbia tutte componenti uguali a 1, cosicchè si possa conoscere con esattezza l’errore commesso nel suo calcolo. Risolvendo il sistema di ordine 20 con il metodo di Gauss senza pivoting si osserva che la soluzione è, in realtà, molto lontana da quella teorica. Questa situazione peggiora prendendo matrici di dimensioni crescenti ed è legata ad un fenomeno che viene detto malcondizionamento. Bisogna infatti ricordare che, a causa degli errori legati alla rappresentazione dei numeri reali, il sistema che l’elaboratore risolve non coincide con quello teorico, poichè alla matrice A ed al vettore b e necessario aggiungere la matrice δA ed il vettore δb (che contengono le perturbazioni legate a tali errori), e che la soluzione ovviamente non è la stessa, pertanto la indichiamo con x+ δx: (A+ δA)(x+ δx) = b+ δb. (3.30) CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 71 Si può dimostrare che l’ordine di grandezza della perturbazione sulla solu- zione è ‖δx‖ ‖x‖ ≤ ‖A‖‖A −1‖ (‖δA‖ ‖A‖ + ‖δb‖ ‖b‖ ) . Il numero K(A) = ‖A‖‖A−1‖, detto indice di condizionamento del sistema, misura le amplificazioni degli errori sui dati del problema (ovvero la misura di quanto aumentano gli errori sulla soluzione). Il caso della matrice di Hilbert è appunto uno di quelli per cui l’indice di condizionamento assume valori molto grandi (di ordine esponenziale) all’aumentare della dimensione, si parla infatti di matrici malcondizionate. Quando ciò non accade si parla invece di matrici bencondizionate. A volte tale caratteristica può dipendere anche dalla scelta dell’algoritmo di risoluzione, ovvero vi sono algoritmi che forniscono risultati meno influenzati dal condizionamento dei dati (eseguendo il metodo di Gauss con pivoting parziale, per esempio, i risultati sono affetti comunque da errori, ma di meno rispetto al metodo di Gauss senza alcuna strategia di pivoting). 3.5 Metodi iterativi per sistemi lineari 3.5.1 Il Metodo di Jacobi Nella applicazioni numeriche spesso i sistemi lineari che si devono risolvere sono di grandi dimensione e hanno una struttura sparsa, cioè una buona parte degli elementi della matrice dei coefficienti sono nulli. Quando si applicano i metodi diretti a tali sistemi succede che le matrici perdono la struttura di sparsità, cioè molti elementi nulli diventano diversi da zero e inoltre si ha il problema di gestire matrici di grosse dimensioni, il che può causare un notevole degrado delle prestazioni dei metodi usati. Per questi motivi si introduce una nuova classe di metodi, detti metodi iterativi. Supponiamo di dover risolvere il sistema Ax = b, con A matrice non singolare, b 6= 0. Assumiamo inoltre che gli elementi diagonali della matrice aii, i = 1, . . . , n, siano diversi da 0. La i−esima equazione del sistema si scrive ai1x1 + ai2x2 + · · ·+ ainxn = bi e, isolando xi risulta: xi = ( bi − n∑ j=1,j 6=i aijxj ) 1 aii i = 1, . . . , n. CAPITOLO 3. METODI NUMERICI PER SISTEMI LINEARI 72 Queste n equazioni sono del tutto equivalenti al sistema di partenza tuttavia la loro forma suggerisce particolari procedimenti iterativi per cercare la so- luzione. A partire da un’approssimazione iniziale x(0) = (x (0) 1 , x (0) 2 , . . . , x (0) n ) si calcola la successione di vettori {x(k)} ponendo x (k+1) i = ( bi − n∑ j=1,j 6=i aijx (k) j ) 1 aii i = 1, . . . , n; (3.31) con k = 0, 1, 2, . . . . La generica componente i−esima del vettore al passo k + 1 è calcolata per mezzo di tutte le componenti del vettore al passo k eccetto la i−esima. Questo procedimento iterativo prende il nome di metodo di Jacobi. 3.5.2 Il Metodo di Gauss-Seidel Una variante del metodo di Jacobi si ottiene osservando che, quando si calcola x (k+1) i si possono utilizzare le approssimazioni x (k+1) j , con j = 1, . . . , i − 1, ottenendo x (k+1) i = 1 aii ( bi − i−1∑ j=1 aijx (k+1) j − n∑ j=i+1 aijx (k) j ) . (3.32) Si ottiene in questo modo il classico Metodo di Gauss-Seidel, dove la compo- nente i−esima al passo k + 1 è calcolata per mezzo delle componenti dalla prima alla i−1-esima al passo k+1 e dalla i+1-esima alla n−esima al passo k. Si deve osservare che entrambi i metodi appena introdotti non utilizzano la componente i−esima al passo k. Per questo si introduce una nuova va- riante che coinvolge tale valore a partire da un parametro ω 6= 0. Si propone lo schema x (k+1) i = ω aii [ bi − i−1∑ j=1 aijx (k+1) j − n∑ j=i+1 aijx (k) j ] + (1− ω)x(k)i . (3.33) Questa classe di metodi prende il nome di Metodi di Rilassamento. Si osserva facilmente che se si pone ω = 1 il metodo di Rilassamento coincide con il metodo di Gauss-Seidel. Tutti i metodi appena descritti sono molto efficaci quando la matrice dei coefficienti ha pochi elementi diversi da zero (ovvero se la matrice è sparsa) CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 75 quello in cui si cerca la funzione g in (4.1) nella forma g(x; a0, a1, . . . , an) = n∑ i=0 aiΦi(x) dove Φi(x), per i = 0, . . . , n, sono funzioni fissate e i valori di ai, i = 0, . . . , n, sono determinati in base alle condizioni di coincidenza di f con la funzione approssimante nei punti di interpolazione (detti anche nodi), xj, cioè si pone f(xj) = n∑ i=0 aiΦi(xj) j = 0, . . . , n. (4.2) Il processo di determinazione degli ai attraverso la risoluzione del sistema (4.2) si chiama metodo dei coefficienti indeterminati. Il caso più studiato è quello dell’interpolazione polinomiale, in cui si pone: Φi(x) = x i i = 0, . . . , n e perciò la funzione approssimante g assume la forma n∑ i=0 aix i, mentre le condizioni di coincidenza diventano a0 +a1x0 +a2x 2 0 + . . . +an−1x n−1 0 +anx n 0 = f(x0) a0 +a1x1 +a2x 2 1 + . . . +an−1x n−1 1 +anx n 1 = f(x1) ... ... ... ... a0 +a1xn +a2x 2 n + . . . +an−1x n−1 n +anx n n = f(xn) (4.3) Le equazioni (4.3) costituiscono un sistema di n + 1 equazioni nelle n + 1 incognite ai, i = 0, . . . , n : V a = y dove la matrice dei coefficienti è V =            1 x0 x 2 0 . . . x n−1 0 x n 0 1 x1 x 2 1 . . . x n−1 1 x n 1 ... ... ... ... ... 1 xn x 2 n . . . x n−1 n x n n            , CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 76 i vettori dei termini noti e delle incognite sono, rispettivamente, y = [f(x0), f(x1), . . . , f(xn)] T e a = [a0, a1, . . . , an] T . Se i nodi xj sono a due a due distinti allora la matrice dei coefficienti del sistema (4.3), detta matrice di Vandermonde, è non singolare e pertanto il problema dell’interpolazione ammette sempre un’unica soluzione. Il metodo dei coefficienti indeterminati consente di trovare la soluzione del problema so- lo risolvendo un sistema lineare che potrebbe avere grandi dimensioni, essere malcondizionato (soprattutto se due nodi sono molto vicini) e comunque non in grado di fornire un’espressione in forma chiusa del polinomio. Per questi motivi descriviamo un modo alternativo per risolvere il problema di interpolazione in grado di fornire l’espressione esplicita del polinomio cercato. 4.2 Il Polinomio Interpolante di Lagrange Al fine di dare una forma esplicita al polinomio interpolante, scriviamo il candidato polinomio nella seguente forma: Ln(x) = n∑ k=0 lnk(x)f(xk) (4.4) dove gli lnk(x) sono per il momento generici polinomi di grado n. Imponendo le condizioni di interpolazione Ln(xi) = f(xi) i = 0, . . . , n deve essere, per ogni i: Ln(xi) = n∑ k=0 lnk(xi)f(xk) = f(xi) ed è evidente che se lnk(xi) =    0 se k 6= i 1 se k = i (4.5) CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 77 allora esse sono soddisfatte. Infatti calcolando il polinomio (4.4) in un generico nodo xi risulta Ln(xi) = n∑ k=0 lnk(xi)f(xk) = i−1∑ k=0 lnk(xi)f(xk) ︸︷︷︸ =0 + lni(xi) ︸︷︷︸ =1 f(xi) + n∑ k=i+1 lnk(xi)f(xk) ︸︷︷︸ =0 = f(xi). Per determinare l’espressione del generico polinomio lnk(x) osserviamo che la prima condizione di (4.5) indica che esso si annulla negli n nodi x0, x1, . . . , xk−1, xk+1, . . . , xn pertanto deve essere lnk(x) = ck n∏ i=0,i 6=k (x− xi) mentre impondendo la seconda condizione di (4.5) lnk(xk) = ck n∏ i=0,i 6=k (xk − xi) = 1 si trova immediatamente: ck = 1 n∏ i=0,i 6=k (xk − xi) . In definitiva il polinomio interpolante ha la seguente forma: Ln(x) = n∑ k=0 ( n∏ i=0,i 6=k x− xi xk − xi ) f(xk). (4.6) Il polinomio (4.6) prende il nome di Polinomio di Lagrange mentre i polinomi: lnk(x) = n∏ i=0,i 6=k x− xi xk − xi ; k = 0, 1, . . . , n si chiamano Polinomi Fondamentali di Lagrange. CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 80 Esempio 4.2.1 Supponiamo di voler calcolare il polinomio interpolante di Lagrange passante per i punti (−1,−1), (0, 1), (1,−1), (3, 2) e (5, 6). Il grado di tale polinomio è 4, quindi definiamo i nodi x0 = −1, x1 = 0, x2 = 1, x3 = 3, x4 = 5, cui corrispondono le ordinate che indichiamo con yi, i = 0, . . . , 4: y0 = −1, y1 = 1, y2 = −1, y3 = 2, y4 = 6. Scriviamo ora l’espressione del polinomio L4(x): L4(x) = l4,0(x)y0 + l4,1(x)y1 + l4,2(x)y2 + l4,3(x)y3 + l4,4(x)y4 (4.8) e calcoliamo i 5 polinomi fondamentali di Lagrange: l4,0(x) = (x− 0)(x− 1)(x− 3)(x− 5) (−1− 0)(−1− 1)(−1− 3)(−1− 5) = 1 48 x(x− 1)(x− 3)(x− 5) l4,1(x) = (x+ 1)(x− 1)(x− 3)(x− 5) (0 + 1)(0− 1)(0− 3)(0− 5) = − 1 15 (x+ 1)(x− 1)(x− 3)(x− 5) l4,2(x) = (x+ 1)(x− 0)(x− 3)(x− 5) (1 + 1)(1− 0)(1− 3)(1− 5) = 1 16 x(x+ 1)(x− 3)(x− 5) l4,3(x) = (x+ 1)(x− 0)(x− 1)(x− 5) (3 + 1)(3− 0)(3− 1)(3− 5) = − 1 48 x(x+ 1)(x− 1)(x− 5) l4,4(x) = (x+ 1)(x− 0)(x− 1)(x− 3) (5 + 1)(5− 0)(5− 1)(5− 3) = 1 240 x(x+ 1)(x− 1)(x− 3) CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 81 −1 0 1 2 3 4 5 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Figura 4.1: Grafico del polinomio l40(x). Sostituendo in (4.8) il valore della funzione nei nodi si ottiene l’espressione finale del polinomio interpolante: L4(x) = −l4,0(x) + l4,1(x)− l4,2(x) + 2l4,3(x) + 6l4,4(x). Se vogliamo calcolare il valore approssimato della funzione f(x) in un’ascissa diversa dai nodi, per esempio x = 2 allora dobbiamo calcolare il valore del polinomio interpolante L4(2). Nelle figure 4.1-4.5 sono riportati i grafici dei cinque polinomi fondamentali di Lagrange: gli asterischi evidenziano il valore assunto da tali polinomi nei nodi di interpolazione. Nella figura 4.6 è tracciato il grafico del polinomio interpolante di Lagrange, i cerchi evidenziano ancora una volta i punti di interpolazione. 4.2.2 Il fenomeno di Runge Nell’espressione dell’errore è presente, al denominatore, il fattore (n + 1)!, che potrebbe indurre a ritenere che, utilizzando un elevato numero di no- di, l’errore tenda a zero ed il polinomio interpolante converga alla funzione CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 82 −1 0 1 2 3 4 5 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Figura 4.2: Grafico del polinomio l41(x). −1 0 1 2 3 4 5 −1.5 −1 −0.5 0 0.5 1 1.5 Figura 4.3: Grafico del polinomio l42(x). CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 85 −5 −4 −3 −2 −1 0 1 2 3 4 5 −0.5 0 0.5 1 1.5 2 x y Fenomeno di Runge Figura 4.7: Il fenomeno di Runge. Per ovviare al fenomeno di Runge si possono utilizzare insiemi di nodi non equidistanti oppure utilizzare funzioni interpolanti polinomiali a tratti (in- terpolando di fatto su intervalli più piccoli e imponendo le condizioni di continuità fino ad un ordine opportuno). function yy=lagrange(x,y,xx); % % La funzione calcola il polinomio interpolante di Lagrange % in un vettore assegnato di ascisse % % Parametri di input % x = vettore dei nodi % y = vettore delle ordinate nei nodi % xx = vettore delle ascisse in cui calcolare il polinomio % Parametri di output % yy = vettore delle ordinate del polinomio % n = length(x); m = length(xx); CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 86 yy = zeros(size(xx)); for i=1:m yy(i)=0; for k=1:n yy(i)=yy(i)+prod((xx(i)-x([1:k-1,k+1:n]))./... (x(k)-x([1:k-1,k+1:n]))))*y(k); end end return 4.3 Minimizzazione del Resto nel Problema di Interpolazione Supponiamo che la funzione f(x) sia approssimata su [a, b] dal polinomio interpolante Ln(x) e siano x0, x1, . . . , xn i nodi di interpolazione. Come già sappiamo se x ∈ [a, b] risulta e(x) = f(x)− Ln(x) = f (n+1)(ξx) (n+ 1)! ωn+1(x) ξx ∈ [a, b] e dove ωn+1(x) = n∏ i=0 (x− xi). Si noti che variando i nodi xi, i = 0, . . . , n, cambia il polinomio ωn+1(x) e di conseguenze cambia l’errore. Ha senso allora porsi il seguente problema: indicato con Pn+1 l’insieme di tutti i polinomi di grado al più n+1 cerchiamo il polinomio p̃ ∈ Pn+1 tale che: max x∈[a,b] |p̃(x)| = min p∈Pn+1 max x∈[a,b] |p(x)|. (4.9) Per dare una risposta a questo problema è essenziale introdurre i Polinomi di Chebyshev di 1a Specie. 4.3.1 Polinomi di Chebyshev I polinomi di Chebyshev Tn(x), n ≥ 0, sono cos̀ı definiti: Tn(x) = cos(n arccos x) (4.10) CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 87 per x ∈ [−1, 1]. Per esempio: T0(x) = cos(0 arccos x) = cos 0 = 1 T1(x) = cos(1 arccos x) = x e cos̀ı via. È possibile ricavare una relazione di ricorrenza sui polinomi di Chebyshev che permette un più agevole calcolo. Infatti, posto arccos x = θ (ovvero x = cos θ) risulta Tn(x) = cosnθ(x). Considerando le relazioni Tn+1(x) = cos(n+ 1)θ = cosnθ cos θ − sinnθ sin θ Tn−1(x) = cos(n− 1)θ = cosnθ cos θ + sinnθ sin θ e sommandole membro a membro, Tn+1(x) + Tn−1(x) = 2 cos θ cosnθ = 2xTn(x) si ricava la seguente relazione di ricorrenza Tn+1(x) = 2x Tn(x)− Tn−1(x), n ≥ 1 (4.11) che, insieme all’espressione dei primi due polinomi, T0(x) = 1, T1(x) = x. consente di calcolare tutti i polinomi di Chebyshev. L’espressione dei primi polinomi è la seguente T2(x) = 2xT1(x)− T0(x) = 2x2 − 1 T3(x) = 2xT2(x)− T1(x) = 4x3 − 3x T4(x) = 2xT3(x)− T2(x) = 8x4 − 8x2 + 1 T5(x) = 2xT4(x)− T3(x) = 16x5 − 20x3 + 5x Le seguenti proprietà dei polinomi di Chebyshev sono di facile dimostrazione: CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 90 Il teorema di minimax stabilisce che, tra tutti i polinomi di grado n definiti nell’intervallo [−1, 1], il polinomio di Chebyshev monico è quello che ha il massimo più piccolo. Supponendo che l’intervallo di interpolazione della funzione f(x) sia appunto [−1, 1] e scegliendo come nodi gi zeri del polinomio di Chebyshev risulta ωn+1(x) = T̃n+1(x) pertanto e(x) = f (n+1)(ξx) (n+ 1)! T̃n+1(x) e, massimizzando tale errore, risulta max x∈[−1,1] |e(x)| ≤ max x∈[−1,1] ∣ ∣ ∣ ∣ f (n+1)(ξx) (n+ 1)! ∣ ∣ ∣ ∣ max x∈[−1,1] |ωn+1(x)| = 1 2n(n+ 1)! max x∈[−1,1] |f (n+1)(ξx)|. La crescita dell’errore può dipendere solo dalla derivata di ordine n+1 della funzione f(x). Se l’intervallo di interpolazione è [a, b] 6= [−1, 1] allora il discorso può essere ripetuto egualmente effettuando una trasformazione lineare tra i due inter- valli, nel modo riportato in Figura 4.9. Calcolando la retta nel piano (x, t) passante per i punti (−1, a) e (1, b): t = b− a 2 x+ a+ b 2 (4.12) detti xk gli zeri del polinomio di Chebyshev Tn+1(x) allora si possono usare come nodi i valori tk = b− a 2 xk + a+ b 2 , k = 0, 1, . . . , n, ovvero tk = b− a 2 cos (2k + 1)π 2(n+ 1) + a+ b 2 k = 0, 1, . . . , n. (4.13) Il polinomio di Chebyshev, traslato nell’intervallo [a, b], è T [a,b] n+1(x) = Tn+1 ( 2x− (b+ a) b− a ) , CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 91 1−1 t x a b • • Figura 4.9: Trasformazione lineare tra gli intervalli [−1, 1] e [a, b]. il cui coefficiente di grado massimo vale 2n 2n+1 (b− a)n+1 = 22n+1 (b− a)n+1 . Se come nodi di interpolazione scegliamo i punti tk dati da (4.13), cioè gli n+ 1 zeri del polinomio T̃ [a,b] n+1(x), allora abbiamo il polinomio monico è T̃ [a,b] n+1(x) = (b− a)n+1 22n+1 Tn+1 ( 2x− (b+ a) b− a ) , considerato che la trasformazione lineare inversa della (4.12) è t = 2x− (b+ a) b− a , x ∈ [a, b]→ t ∈ [−1, 1] quindi per l’errore dell’interpolazione vale la seguente maggiorazione: max x∈[a,b] |e(x)| ≤ max x∈[a,b] ∣ ∣ ∣ ∣ f (n+1)(ξx) (n+ 1)! ∣ ∣ ∣ ∣ max x∈[a,b] |T̃ [a,b]n+1(x)| = max x∈[a,b] ∣ ∣ ∣ ∣ f (n+1)(ξx) (n+ 1)! ∣ ∣ ∣ ∣ (b− a)n+1 22n+1 . CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 92 −5 −4 −3 −2 −1 0 1 2 3 4 5 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x y Interpolazione su nodi di Cheyshev Figura 4.10: Interpolazione su nodi di Chebyshev. Nella Figura 4.10 sono raffigurati la funzione di Runge ed il polinomio in- terpolante di Lagrange di grado 10 calcolato prendendo come nodi gli zeri del polinomio di Chebyshev di grado 11. Si può osservare la differenza con la Figura 4.7. Di seguito viene riportato il codice per tracciare il grafico del polinomio interpolante la funzione di Runge nei nodi di Chebyshev in un intervallo scelto dall’utente. clear format long e a = input(’Inserire estremo sinistro ’); b = input(’Inserire estremo destro ’); n = input(’Inserire il numero di nodi ’); % % Calcolo del vettore dei nodi di Chebyshev % x = (a+b)/2+(b-a)/2*cos((2*[0:n-1]+1)*pi./(2*n)); xx = linspace(a,b,200); y = 1./(x.^2+1); yy = 1./(xx.^2+1); CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 95 dove yi, i = 0, . . . , n+ 1, sono n+ 2 assegnati valori. Indichiamo con si(x) la restrizione della spline nell’intervallo [xi, xi+1], in cui s′′i (x) è una funzione lineare mentre s (3) i (x) è una costante, quindi s′′i (x) = s ′′ i (xi) + s (3) i (xi)(x− xi) (4.15) ovvero, posto Mi = s ′′ i (xi) ci = s (3) i (xi) abbiamo s′′i (x) = Mi + ci(x− xi). (4.16) Valutando (4.16) in xi+1 si ottiene ci = Mi+1 −Mi hi , hi = xi+1 − xi. (4.17) Scriviamo lo sviluppo in serie di Taylor di si(x) prendendo come punto iniziale xi : si(x) = si(xi) + s ′ i(xi)(x− xi) + s′′i (xi) (x− xi)2 2 + s (3) i (xi) (x− xi)3 6 , (4.18) sostituiamo i valori delle derivate seconda e terza, e calcoliamola in xi+1 si(xi+1) = si(xi) + s ′ i(xi)(xi+1 − xi) +Mi (xi+1 − xi)2 2 + ci (xi+1 − xi)3 6 e, imponendo le condizioni di interpolazione e sostituendo il valore dell’am- piezza dei sottointervalli, si ottiene yi+1 = yi + s ′ i(xi)hi +Mi h2i 2 + ci h3i 6 da cui yi+1 − yi hi = s′i(xi) +Mi hi 2 + ci h2i 6 . (4.19) Scriviamo ora lo sviluppo in serie di Taylor di si−1(x) prendendo come punto iniziale xi : si−1(x) = si−1(xi)+ s ′ i−1(xi)(x−xi)+ s′′i−1(xi) (x− xi)2 2 + s (3) i−1(xi−1) (x− xi)3 6 CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 96 e sostituiamo il valori della derivate seconda (che à uguale a Mi in quanto è continua) e della derivata terza (che invece è uguala a ci−1 in quanto è discontinua), e poniamo x = xi−1 e calcoliamola in xi−1, si−1(xi−1) = si−1(xi)+s ′ i−1(xi)(xi−1−xi)+Mi (xi−1 − xi)2 2 +ci−1 (xi−1 − xi)3 6 . Imponendo le condizioni di interpolazione anche sul nodo xi in modo tale da assicurare la continuità della spline si ottiene yi−1 = yi − s′i−1(xi)hi−1 +Mi h2i−1 2 − ci−1 h3i−1 6 yi−1 − yi hi−1 = −s′i−1(xi) +Mi hi−1 2 − ci−1 h2i−1 6 . (4.20) Osserviamo dalla relazione (4.19) che s′i(xi) può essere calcolata se sono noti i valori Mi. Di conseguenza la spline è completamente determinata se si conoscono i valori M0, M1, . . . , Mn+1 (che sono detti momenti). A questo punto imponendo le condizioni di continuità della derivata prima, ovvero s′i−1(xi) = s ′ i(xi) sommando le equazioni (4.19) e (4.20) le derivate prime si semplificano rica- vando l’equazione Mi hi + hi−1 2 − ci−1 h2i−1 6 + ci h2i 6 = yi+1 − yi hi + yi−1 − yi hi−1 , sostituendo l’espressione delle derivate terze nei due intervalli ci−1 = Mi −Mi−1 hi−1 , ci = Mi+1 −Mi hi 3Mi(hi+hi−1)+(Mi+1−Mi)hi−(Mi−Mi−1)hi−1 = 6 ( yi+1 − yi hi + yi−1 − yi hi−1 ) ottenendo, per i = 1, . . . , n, le equazioni hi−1Mi−1 + 2(hi + hi−1)Mi + hiMi+1 = 6 ( yi+1 − yi hi + yi−1 − yi hi−1 ) (4.21) che rappresentano un sistema lineare di n equazioni nelle n+2 incognite M0, M1, . . . ,Mn+1. Per risolvere il problema delle condizioni mancanti ci sono diverse possibili scelte: CAPITOLO 4. INTERPOLAZIONE DI DATI E FUNZIONI 97 1. Il modo piu semplice e quello di imporre che la spline abbia momenti nulli agli estremi dell’intervallo, ovvero M0 = Mn+1 = 0 ottenendo la cosiddetta spline cubica naturale; 2. In alternativa si possono fissano i valori della derivata prima negli estremi, imponendo che sia s′(x0) = k0 s′(xn+1) = kn+1 con k0 e kn+1 valori assegnati: in tal caso la spline interpolante prende il nome di spline cubica completa; 3. Imponendo s′(x0) = s ′(xn+1), s ′′(x0) = s ′′(xn+1) si ottiene invece la cosiddetta spline cubica periodica. Considerando la spline cubica naturale interpolante il sistema (4.21) si scrive: Am = b dove A =        2(h0 + h1) h1 h1 2(h1 + h2) h2 . . . . . . . . . hn−2 2(hn−2 + hn−1) hn−1 hn−1 2(hn−1 + hn)        (4.22) e m =          M1 M2 M3 ... Mn−1 Mn          b =          b1 b2 b3 ... bn−1 bn          con bi = 6 ( yi+1 − yi hi + yi−1 − yi hi−1 ) i = 1, 2, . . . , n.