INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore...

150

Transcript of INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore...

Page 1: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Indice

��� Premessa � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ���� Introduzione � � � � � � � � � � � � � � � � � � � � � � � � � � � � �

����� Un semplice esempio � � � � � � � � � � � � � � � � � � � �

� Teoria dell� errore �M�Frontini� ���� De�nizioni �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ��� Numeri macchina ��oating point �M�F� � � � � � � � � � � � � ����� Propagazione degli errori �M�F� � � � � � � � � � � � � � � � � � ��

����� Somma algebrica � � � � � � � � � � � � � � � � � � � � � ������� Prodotto � � � � � � � � � � � � � � � � � � � � � � � � � � ������� Divisione � � � � � � � � � � � � � � � � � � � � � � � � � � ������� Radice quadrata � � � � � � � � � � � � � � � � � � � � � �

��� Alcuni esempi �M�F� � � � � � � � � � � � � � � � � � � � � � � � � ����� Cancellazione numerica � � � � � � � � � � � � � � � � � � � ����� Instabilit�a algoritmica � � � � � � � � � � � � � � � � � � ������� Sensitivit�a del problema � � � � � � � � � � � � � � � � � ��

�� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ���� �� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ��

� Algebra lineare numerica �M�Frontini� ����� Sistemi lineari �M�F� � � � � � � � � � � � � � � � � � � � � � � � �

����� Richiami di algebra lineare � � � � � � � � � � � � � � � � ������ Operazioni sulle matrici � � � � � � � � � � � � � � � � � ��

��� Metodi diretti �M�F� � � � � � � � � � � � � � � � � � � � � � � � ������� Il metodo di eliminazione di Gauss � � � � � � � � � � � ������� Decomposizione LU � � � � � � � � � � � � � � � � � � � � ������� Decomposizione di Cholesky � � � � � � � � � � � � � � � �

��� Analisi dell�errore �M�F� � � � � � � � � � � � � � � � � � � � � � ������� Richiami sulle norme di vettore � � � � � � � � � � � � � ��

Page 2: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� INDICE

����� Stima dell�errore � � � � � � � � � � � � � � � � � � � � � ����� Metodi iterativi �M�F� � � � � � � � � � � � � � � � � � � � � � � ��

����� Ra�namento iterativo di Wilkinson � � � � � � � � � � � � ����� Metodo di Jacoby � � � � � � � � � � � � � � � � � � � � � � ����� Metodo di Gauss�Seidel � � � � � � � � � � � � � � � � � ������� Metodi di rilassamento SOR � � � � � � � � � � � � � � � ����� Metodi non stazionari � � � � � � � � � � � � � � � � � � ��

�� Sistemi sovradeterminati �M�F� � � � � � � � � � � � � � � � � � ����� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ��

����� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ����� Autovalori di matrici �M�F� � � � � � � � � � � � � � � � � � � � ��

����� Richiami e de�nizioni � � � � � � � � � � � � � � � � � � � ������� Metodi locali � � � � � � � � � � � � � � � � � � � � � � � ������ Metodi globali � � � � � � � � � � � � � � � � � � � � � � � ������ Analisi dell�errore � � � � � � � � � � � � � � � � � � � � � �

��� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ����� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ��

� Zeri di funzioni �M�Frontini� ���� Metodo di bisezione �M�F� � � � � � � � � � � � � � � � � � � � � ��

����� Falsa posizione � � � � � � � � � � � � � � � � � � � � � � ����� Metodi di punto �sso �M�F� � � � � � � � � � � � � � � � � � � � ����� Metodo di Newton �M�F� � � � � � � � � � � � � � � � � � � � � ��

����� Metodo delle secanti � � � � � � � � � � � � � � � � � � � ������� Metodo di Ste�ensen � � � � � � � � � � � � � � � � � � � ��

��� Sistemi non lineari �M�F� � � � � � � � � � � � � � � � � � � � � ������� Metodi di punto �sso � � � � � � � � � � � � � � � � � � � ������� Metodo di Newton per sistemi � � � � � � � � � � � � � � ��

�� Zeri di polinomi �M�F� � � � � � � � � � � � � � � � � � � � � � � ���� �� Schema di Horner � � � � � � � � � � � � � � � � � � � � � ���� �� Radici multiple � � � � � � � � � � � � � � � � � � � � � � ��

��� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ������� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ��

Teoria dell�approssimazione �M�Frontini� ����� Interpolazione �M�F� � � � � � � � � � � � � � � � � � � � � � � � ��

����� Polinomi di Lagrange � � � � � � � � � � � � � � � � � � � ������� Sistema di Vandermonde � � � � � � � � � � � � � � � � � ��

Page 3: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

INDICE �

����� Stima dell�errore � � � � � � � � � � � � � � � � � � � � � � ����� Di�erenze divise � � � � � � � � � � � � � � � � � � � � � � ������ Interpolazione di Hermite � � � � � � � � � � � � � � � � ������� Spline � � � � � � � � � � � � � � � � � � � � � � � � � � � ������ Derivazione numerica � � � � � � � � � � � � � � � � � � � �

��� Minimi quadrati �M�F� � � � � � � � � � � � � � � � � � � � � � ����� Polinomi trigonometrici � � � � � � � � � � � � � � � � � �

��� Riepilogo � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ������ Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � �

� Formule di Quadratura �M�Frontini� � � �� Formule di Newton�Cotes �M�F� � � � � � � � � � � � � � � � � ���

���� Formule composite � � � � � � � � � � � � � � � � � � � � ��� �� Formule adattive �M�F� � � � � � � � � � � � � � � � � � � � � � �� �� Formule gaussiane �M�F� � � � � � � � � � � � � � � � � � � � � ���

���� Integrali impropri � � � � � � � � � � � � � � � � � � � � � �� �� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � ���

���� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ���

Equazioni di�erenziali ordinarie �M�Frontini� ������ Metodi ad un passo �M�F� � � � � � � � � � � � � � � � � � � � � ���

����� Metodi di Taylor � � � � � � � � � � � � � � � � � � � � � �������� Metodi Runge�Kutta � � � � � � � � � � � � � � � � � � � �������� Sistemi del primo ordine � � � � � � � � � � � � � � � � � �������� Stabilit�a numerica � � � � � � � � � � � � � � � � � � � � ���

��� Metodi a pi�u passi �M�F� � � � � � � � � � � � � � � � � � � � � �������� Metodi di Adams � � � � � � � � � � � � � � � � � � � � � �������� Stabilit�a numerica � � � � � � � � � � � � � � � � � � � � ���

��� Riepilogo �M�F� � � � � � � � � � � � � � � � � � � � � � � � � � �������� Esercizi � � � � � � � � � � � � � � � � � � � � � � � � � � ���

Page 4: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� INDICE

��� Premessa

La continua richiesta da parte degli studenti dei corsi di Calcolo Numericoda me tenuti mi ha spinto� dopo non pochi ripensamenti� alla stesura diqueste brevi note� Sebbene sia sempre convinto che esistono sulla piazzamolti ottimi libri di Calcolo Numerico� mi sono altres�i convinto che per glistudenti pu�o essere di grande aiuto� per la preparazione dell� esame� un testoguida che contenga gli argomenti svolti nel corso dell� anno� Per questaseconda ragione i presenti appunti sono stati stesi� giorno per giorno� durantelo svolgimento del corso di Calcolo Numerico ��� annualita� del secondosemestre ��� Nonostante l�aspetto tipogra�co �grazie al �favoloso� TeXrestano dei semplici appunti per cui si invitano �n d�ora gli studenti adepurare il testo dei vari errori di stampa e dalle eventuali imprecisioni�

Questo sforzo �e stato fatto principalmente per gli studenti che� a causadelle sovapposizioni d�orario� non riescono a seguire le lezioni� Per quelliche seguono un invito a segnalarmi le discrepanze fra appunti e lezioni� Leprime bozze di questi appunti saranno disponibili� per gli studenti del corso�sulla mia pagina Web �http���marfro�mate�polimi�it sotto forma di �le post�script o DVI �sperimentale utilizzo didattico del potente mezzo informatico�Internet�

��� Introduzione

Gli argomenti del corso� divisi in � capitoli� sono i seguenti�

�� Teoria dell� errore

�� Algebra lineare numerica

Sistemi lineari

Autovalori

�� Zeri di equazioni e sistemi non lineari

�� Teoria dell� approssimazione

Interpolazione

Minimi quadrati

Page 5: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� INTRODUZIONE

� Formule di quadratura

�� Equazioni di�erenziali ordinarie

Nel capitolo �� dopo aver introdotto i concetti di errore relativo� di nu�mero macchina e di epsilon macchina� viene studiata la propagazione dell�errore relativo nelle operazioni elementari di macchina� Vengono presentatianche i concetti di stabilit�a algoritmica e di condizionamento di un problemanumerico� Con semplici esempi si introducono anche il concetto di erroredi troncamento� di convergenza ed accuratezza �o stima dell� errore per unprocesso numerico iterativo�

Nel capitolo �� vengono presentati i principali metodi� sia di tipo direttoche di tipo iterativo� per la risoluzione di sistemi lineari� Per i metodi diretti siporr�a l� accento sulla loro complessit�a computazionale� mentre per gli iterativisi analizzeranno in dettaglio le condizioni di convergenza e maggiorazionedell� errore� Il problema del condizionamento di matrici e la sua in�uenzasulla �buona� risolubilit�a di un sistema lineare verr�a trattato in dettaglio�Viene fatto un cenno sulla risoluzione di sistemi lineari nel senso dei minimiquadrati�

Per il calcolo degli autovalori e autovettori di matrici si presenta in det�taglio il metodo delle potenze ed alcune sue varianti� facendo solo un cennoai metodi globali �algoritmo QR�

Nel capitolo �� vengono presentati alcuni metodi iterativi generali peril calcolo degli zeri di funzioni e sistemi non lineari� Dopo aver introdottoil classico metodo di bisezione� viene presentata la famiglia dei metodi dipunto �sso �metodi di ordine k� con le relative condizioni di convergenzaed accuratezza� Il metodo di Newton e le sue varianti vengono studiati indettaglio anche per sistemi� Per le equazioni algebriche vengono presentatimetodi basati sull� uso della matrice di Frobenious o companion matrix� Vie�ne inoltre discusso il problema della sensitivit�a delle radici di un polinomioalle variazioni dei suoi coe�cienti�

Nel capitolo �� si fornisce una breve introduzione alle problematiche dell�approssimazione� Dopo aver richiamato i concetti base di miglior approssi�mazione si d�a spazio all� interpolazione polinomiale con qualche cenno all�interpolazione mediante funzioni spline� Per quanto riguarda l� approssima�zione nel senso dei minimi quadrati verr�a considerato il caso polinomiale e diFourier�

Nel capitolo � vengono presentate le formule di tipo interpolatorio diNewton�Cotes semplici e composte� de�nendone l� accuratezza e provando�

Page 6: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� INDICE

ne la convergenza �composte� Si far�a inoltre un breve cenno alle formuleGaussiane e alla loro propriet�a di ottimalit�a�

Nel capitolo �� si forniscono gli strumenti essenziali per a�rontare in modocritico la risoluzione dei problemi di�erenziali che tanto interessano i moder�ni ingegneri� Vengono presentate in dettaglio le formule ad un passo �qua�li i metodi Runge�Kutta dandone una semplice deduzione e so�ermandosisui concetti di consistenza� stabilit�a e convergenza� Per le formule a pi�upassi verr�a fatta una breve introduzione al �ne di evidenziarne gli aspetticomputazionali�

����� Un semplice esempio

Prima di addentrarci nel vivo dello studio dei metodi numerici mi sembraopportuno dare una pur semplice giusti�cazione dell� importanza di tali me�todi� A tal �ne pu�o essere utile considerare il seguente problema modello�modello di Malthus� �

y� � �yy�� � y�� y� � �

��

la cui soluzione �e y�t � y�e�t� Se il problema �� simula la crescita di una

colonia di batteri� pu�o essere interessante individuare il tasso di crescita ��A tal �ne �e su�ciente conoscere il valore y�t� � y�� per risolvere�

y�e�t� � y�t� ��

mediante la formula�

� ��

t�ln

�y�t�y�

���

��e necessaria una calcolatrice tascabile con la funzione logaritmo�Se si complica un poco il modello includendo il fenomeno di immigrazio�

ne ed emigrazione l� equazione di�erenziale del problema di Cauchy �� sicomplica di poco divenendo

y� � �y � d ��

dove la costante d tiene conto del nuovo fenomeno� La soluzione di �� tenutoconto delle condizioni iniziali date in �� �e quindi�

y�t � e�t�y� �

d

��� e��t

���

Page 7: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� INTRODUZIONE �

dalla � si ha�

e��t�

y�t� � y� �d

��� e��t

��

equazione NON lineare in � per la cui risoluzione sono necessari metodinumerici�

Il semplice esempio riportato suggerisce uno schema generale di lavoroquando si debba risolvere un problema reale�

Schema risolutivo di un problema

Problema realeerrori�

modello

Corsi speci�ci

Modello matematico � Soluzione analitica �errori �troncamento

Modello numerico � � �nito passierrori � arrotondamento

Algoritmo� programma � � �nito bit�

Soluzione numerica�

�l

Analisi risultati

Page 8: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� INDICE

Page 9: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo �

Teoria dell� errore �M�Frontini�

��� De�nizioni �M�F��

Introduciamo alcune de�nizioni che ci serviranno in seguito� Dati due numerireali x e x dove x �e il valore �vero� e x �e il valore approssimato �di x�de�niamo�

De�nizione ����� Errore assoluto la quantit�a�

jx� xj � ex

De�nizione ����� Errore relativo la quantit�a�

jx� xjjxj �

exjxj � �x

De�nizione ����� Numero troncato a t cifre chopping� valore che si ot�tiene ignorando le cifre dopo la t�esima� Es� ��� � �e il troncato alla quartacifra di ��� ��������

De�nizione ���� Numero arrotondato a t cifre rounding� valore che siottiene ignorando le cifre dopo la t�esima se la �t � ��esima �e in base ��� � ed aumentando la t�esima di una unit�a se la �t� ��esima �e in base ��� � Es� ��� � �e l� arrotondato alla quarta cifra di ��� ��������

De�nizione ����� Intervallo d�indeterminazione� intervallo in cui �e com�preso il valore esatto

x� ex � x � x � ex

Page 10: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� TEORIA DELL� ERRORE �M�FRONTINI�

Dalle de�nizioni date per l� errore assoluto e relativo si vede subito che �el� errore relativo quello che pi�u interessa nelle applicazioni �numero di cifreesatte� si pu�o osservare altres�i che se x �e �piccolo� in modulo �prossimo azero l� errore assoluto �e sicuramente pi�u signi�cativo� Per questo motivosi preferisce usare nella pratica la seguente de�nizione di errore �relativo abuon senso��

jx� xj� � jxj �

ex� � jxj � �rx

dove il pedice r sta per �ragionevole�� �E facile osservare che per x �grandi��rx � �x mentre per x �piccoli� si ha �rx � ex� In questo modo un �rx � ����

�e indice di una buona stima di x sia che esso sia �grande� o �piccolo��L� importanza dell� uso dell� errore relativo viene ra�orzata dal fatto che

i numeri che vengono gestiti da un elaboratore elettronico non sono i classicinumeri reali dell� analisi ma sono un sottoinsieme �nito di questi�

��� Numeri macchina ��oating point� �M�F��

�E nota a tutti la rappresentazione posizionale dei numeri per cui�

�� ��� � � ��� � � ��� � ��� � � ���� � � ����

tale rappresentazione non �e utile nel calcolo scienti�co in quanto si devepresupporre a priori quante cifre utilizzare per la parte intera e per la partedecimale� Pi�u conveniente �e la rappresentazione esponenziale per cui�

�� ��� � ���� �� ��� � ����� �� ��� � ��� �� ��� � ���

se si usa la rappresentazione esponenziale non si ha unicit�a a meno di con�siderarla normalizzata ovvero la mantissa deve essere minore di � e la cifrapi�u signi�cativa deve essere diversa da zero �la prima delle � rappresentazioniesponenziali date�

Quanto detto per la base �� si estende in modo naturale alla base ��la base naturale dei calcolatori� Pu�o essere il caso di far osservare chela rappresentazione di un numero razionale decimale �nito pu�o� in base ��generare un numero razionale binario in�nito �periodico� come illustrato dalseguente esempio�

����� ���� ��

���

����� �������

Page 11: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� NUMERI MACCHINA �FLOATING POINT� �M�F�� ��

come �e facile veri�care eseguendo la divisione fra �� e �����Dato un numero reale x indicheremo d�ora in poi con x il suo rappresen�

tante nell� insieme dei numeri macchina� dove

De�nizione ����� Numero macchina �e una stringa di bit della forma�

x � �d��

�d���

�d���

� ��� �dt�t

� �e

dove�� � di � �� i���������t d� �� �� normalizzatoL � e � U� L ed U sono i valori min e max che pu�o assumere e � Z�t �e il numero di cifre utilizzate per la mantissa��

Da quanto sopra si evincono le seguenti propriet�a dell� insieme F deinumeri �oating point�

�� F R

�� F �e limitato tra un minimo ed un massimo �xmin� xmax � F

�� F �e formato da un numero �nito di elementi pari a�

Card�F � �t �U � L � � � �

�� Gli elementi dell� insieme F non sono equidistanziati sulla retta reale�

Queste propriet�a di F generano le seguenti conseguenze

�� Due numeri reali distinti x ed y possono avere lo stesso rappresentantein F �

�� Esistono numeri reali per cui non esiste il rappresentate in F �Over��ow�

�� �x � � �� � � �e lo zero macchina� x � R �e rappresentato in F da x � ��Under�ow�

�� L� errore relativo di rappresentazione di un numero reale mediante ilsuo rappresentante in F �e sempre lo stesso�

Page 12: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� TEORIA DELL� ERRORE �M�FRONTINI�

� Dati due numeri in F la loro somma pu�o non appartenere ad F �

L� aritmetica di un elaboratore digitale pu�o quindi essere caratterizzatadalla terna �t� L� U� Di seguito rappresentiamo l� insieme F che si ottieneconsiderando un elaboratore che ha le seguenti caratteristiche�

�t� L� U � ������ �

il numero degli elementi di F �e quindi ��� le quattro di�erenti rappresenta�zioni della mantissa sono�

f���� � ��� ���� � �� ���� � ��� ���� � ��ged i quattro esponenti� �

��

�� �� �

�per cui i �� numeri positivi rappresentabili sono�

f���� ��� ���� ���� ���� ��� ���� ���� ��� �� ��� ��� ��� �� ��� ��ga cui vanno aggiunti i �� opposti e lo zero�

Dalla de�nizione di insieme F �e facile dedurre l� errore che si commettequando si inserisce un numero reale in un elaboratore� precisamente proviamoche�

Teorema ����� L�errore dovuto alla memorizzazione all�interno di un ele�boratore binario �e�

jx� xjjxj �

����t chopping��t rounding

dove t �e il numero di cifre utilizzato per la mantissa�

Dimostrazione �����

x � �d�d����dt��� �e mentre x � �d�d����dt �e

con d� �� � �oating point normalizzato� per cui

jx� xj ��

��t�e chopping��t�e�� rounding

ed essendojxj � �e��

segue la tesi� �

Page 13: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� PROPAGAZIONE DEGLI ERRORI �M�F�� ��

Essendo l� errore relativo legato al numero di cifre utilizzate per la man�tissa �e chiaro perch�e l� uso della doppia precisione possa risultare utile nellacomputazione�

La precisione di macchina pu�o essere caratterizzata dalla seguente de��nizione

De�nizione ����� Si chiama epsilon macchina o precisione di macchinail pi�u piccolo numero positivo di F che sommato ad � fornisce un risultatomaggiore di ��

eps �� fminx � F j � � x � �g ��

Un semplice programma �MATLAB like che fornisce eps �e il seguente

eps���while �eps�� � ��

eps���eps

end

eps���eps

�E importante non confondere il concetto di epsilon macchina con quellodi zero macchina de�nito da�

De�nizione ����� Lo zero macchina �e il pi�u piccolo numero macchina invalore assoluto

zero macchina �� fmin jxj � x � Fg

e dipende dal minimo valore attribuibile all� esponente��

�� Propagazione degli errori �M�F��

Ci proponiamo ora di analizzare come l� errore relativo �errore di rappresen�tazione si propaga nelle operazioni elementari� Siano dati due numeri realix� y � R �non nulli ed i loro �rappresentanti in F�� x� y � F � si ha�

x � x�� � �x� j�xj � epsy � y�� � �y� j�yj � eps

Page 14: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� TEORIA DELL� ERRORE �M�FRONTINI�

dove eps �e la precisione di macchina� Ricordando lo sviluppo in serie diTaylor di una funzione di due variabili arrestato al primo ordine si ha

f�x� y � f�x�� � �x� y�� � �y � f�x� y � fx�x� yx�x � fy�x� yy�y

per cui

�f �f�x� y� f�x� y

f�x� y� xfx�x � yfy�y

f�x� y����

Mediante la ���� �e facile dedurre le seguenti valutazioni per le operazionielementari

����� Somma algebrica

Essendo f�x� y � x y� supposto x y �� �� si ha�

�x y� �x y

x y� �x�y �

x�x y�yx y

�x

x y�x y

x y�y ����

dalla ���� �e facile osservare che l� errore nella somma algebrica viene am�pli�cato �se x y � � dai fattori x

x�y e yx�y � che possono essere �grandi�

�cancellazione numerica�

����� Prodotto

Essendo f�x� y � x y� si ha�

�x y� �x y

x y � �x�y � xy�x � xy�yx y � �x � �y ����

dalla ���� �e facile osservare che l� errore relativo nel prodotto non vieneampli�cato�

����� Divisione

Essendo f�x� y � xy� supposto y �� �� si ha�

�xy� �xy

xy� �x�y �

�yx�x � x

y�y�y

xy� �x � �y ����

dalla ���� �e facile osservare che l� errore relativo nella divisione non vieneampli�cato�

Page 15: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� ALCUNI ESEMPI �M�F�� �

����� Radice quadrata

Sia f�x �px� supposto x �� �� si ha�

px�pxp

x� �px �

��pxx�xpx

��

��x ���

e la ��� ci mostra come l� errore relativo venga dimezzato� �x �� �� quandosi estrae la radice quadrata�

L� analisi dei risultati precedenti ci fa capire che non tutte le operazionidi macchina si comportano allo stesso modo per quanto riguarda la propaga�zione degli errori� ma esistono operazioni pi�u� �delicate� �somma algebricaed altre pi�u �robuste� �prodotto� estrazione di radice�

�E doveroso osservare che le stime fornite sono �deterministiche�� Nellapratica� quando si eseguono milioni di operazioni macchina� gli arrotonda�menti nelle operazioni migliorano le cose �non abbiamo qui il tempo di con�siderare la propagazione degli errori da un punto di vista statistico per cui�a parte casi patologici� le computazioni forniscono �mediamente� risultatiragionevoli�

�� Alcuni esempi �M�F��

Diamo ora alcuni esempi di come le operazioni di macchina e gli algoritmi dicalcolo possano in�uenzare la bont�a di un risultato�

����� Cancellazione numerica

Prende il nome di cancellazione numerica quel fenomeno che si manifestaquando si sommano due numeri macchina quasi uguali in modulo ma di segnoopposto� In questo caso �cfr� ���� la somma di macchina �e �instabile� enon gode della proprieta associativa� Consideriamo il seguente esempio� sianoa � �������� �e� �� b � ���������e� �� c � �����������e� � tre numeridati� Supponiamo di operare con una calcolatrice che lavori in base �� con �cifre signi�cative �cifre della mantissa� Veri�chiamo che

a � �b � c �� �a � b � c

infatti si ha�

a � �b � c � �������� �e� � � ����������e� � � ����������e� ��a � b � c � �������� �e� � � ������������e � � � ����������e� �

Page 16: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� TEORIA DELL� ERRORE �M�FRONTINI�

����� Instabilit�a algoritmica

�E quel fenomeno per cui un algoritmo di calcolo propaga �ampli�ca gli errorisui dati� Supponiamo di dover fornire una tabella contenente i primi �� valoridei seguenti integrali

En �

Z �

xnex��dx� n � �� �� ��� ����

integrando per parti si ottiene la relazione ricorrente

En � �� nEn��� n � �� �� ��� ����

sapendo che E� � �e� �e facile ricavare i valori cercati dalla ����� Sfortuna�

tamente� pur dovendo essere En � �� �n� la ���� fornisce� utilizzando uncomune personal computer� valori negativi gi�a per n � � �

Se si inverte la ricorrenza ���� nella forma

En�� ��� En

n� n � ���� �� � ����

ponendo E�� � �� si possono ottenere i primi �� valori di En alla precisione dimacchina� Una pi�u attenta analisi delle ���� e ���� mostra come nel primocaso l� errore su E� venga ampli�cato� ad ogni passo di un fattore che� dopon passi� �e pari ad n�� mentre nel secondo l� errore su E�� � �� viene ridottoad ogni passo di un fattore che� dopo n passi� �e pari ad �

n�

����� Sensitivit�a del problema

Esistono problemi che sono� per loro natura� sensibili alle variazioni sui dati�Un semplice esempio �e fornito dalla ricerca delle radici del seguente polinomiodi secondo grado

�x� �� � ���� ���

che puo�� ovviamente� essere scritto nella forma

x� � �x � ��� ���� � � �����

le radici di ����� sono�

x � � ����

Page 17: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� RIEPILOGO �M�F�� ��

se nella ��� sostituiamo ���� con ����� introduciamo un errore di �����

sul termine noto della ����� le cui radici divengono

x � � � ����

per cui si induce un errore dell� ordine di ���� sulle radici� Tutto questo ope�rando �con carta e matita� �non ci sono errori dovuti alla rappresentazionedei numeri macchina�

��� Riepilogo �M�F��

Dalle precedenti note emerge che�

�� gli errori sono sempre presenti nella computazione ��nitezza della pa�rola di un elaboratore�

�� gli elaboratori che arrotondano sono meglio di quelli che troncano�

�� bisogna utilizzare algoritmi stabili per non propagare gli errori�

�� esistono problemi sensibili alle variazione sui dati �mal condizionatiche vanno trattati con molta attenzione�

� l� uso della doppia �multipla precisione �e utile come veri�ca dei risul�tati pi�u che come strumento naturale di lavoro �aumento del tempo dicalcolo� ine�cacia sui problemi mal condizionati�

����� Esercizi

Esercizio ����� Qual �e il risultato della seguente sequenza di istruzioni MA�TLAB �

x��

while x � ��

x�x�����

x sqrt�x�

end

Esercizio ����� Qual �e il risultato della seguente sequenza di istruzioni MA�TLAB �

Page 18: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� TEORIA DELL� ERRORE �M�FRONTINI�

x��

while x � ��

x�x����

x sqrt�x�

end

Esercizio ����� Dire come si opererebbe per calcolare la seguente espressio�ne

f�x � x�px � ��px

per valori di x � ��n� n � � � �

Esercizio ���� Dire come si opererebbe per calcolare la seguente espressio�ne

f�x ��� cos�x

x�

per valori di x � ���n� n � � � �

Esercizio ����� Commentare il gra�co che si ottiene eseguendo il plot inMATLAB della funzione

f�x � x� � �x � � x� � ��x� � � x� � �x � �

per ��� � x � ������

Esercizio ���� Dovendo calcolare z �px� � y� si suggeriscono le seguenti

formulejxjp� � �yx�� � � jyj � jxjjyj

p� � �xy�� � � jxj � jyj

giusti�care la bont�a delle formule proposte�

Page 19: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo �

Algebra lineare numerica�M�Frontini�

��� Sistemi lineari �M�F��

����� Richiami di algebra lineare

Per sistema lineare si intende un insieme di m equazioni in n incognite dellaforma �����

a��x�� a��x�� ���� a�nxn � b�a��x�� a��x�� ���� a�nxn � b���� ��� ��� ��� ���am�x�� am�x�� ���� amnxn � bm

����

dove aij � R �eventualmente a C� bi � R �eventualmente a C� In formamatriciale ���� si pu�o scrivere

Ax � b

dove A �e una matrice in Rm�n di elementi aij� b �e un vettore in Rm ed x �e inRn�

Richiamiamo le seguenti de�nizioni� perch�e utili in seguito�

De�nizione ����� Vettore �e una n�pla ordinata di numeri in colonna�

De�nizione ����� Matrice �e una tabella rettangolare o quadrata di numeriinsieme di vettori colonna�

Page 20: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

De�nizione ����� Versore i�esimo ei �e un vettore colonna formato da tutti� tranne l�elemento i�esimo che vale ��

ei �

����������

�����������

������������ i

De�nizione ���� Matrice identit�a I �e la matrice formata dagli n versoriei� i���������n�

De�nizione ����� Matrice diagonale �e una matrice con gli elementi aii �� ��i���������n e gli altri nulli�

De�nizione ���� Matrice tridiagonale �e una matrice per cui solo aii �� ��ai��i �� �� aii�� �� ��

T �

�������

a�� a�� � � � � �

a�� a�� a��� � �

���

�� � � � � � � � � �

���� � � an��n�� an��n�� an��n

� � � � � ann�� ann

��������De�nizione ����� Matrice bidiagonale inferiore superiore �e una matriceper cui solo aii �� �� ai��i �� �� aii �� �� aii�� �� ���������

a�� � � � � � � � �

a�� a��� � � � � �

���

�� � � � � � � � �

������

� � � an��n�� an��n�� �� � � � � ann�� ann

�������� �

�������

a�� a�� � � � � �

� a�� a��� � �

������

� � � � � � � � � ����

� � � � � � an��n�� an��n� � � � � � � � ann

��������

Page 21: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� SISTEMI LINEARI �M�F�� ��

De�nizione ����� Matrice triangolare inferiore superiore �e una matriceper cui aij � �� j � i aij � �� j � i��������

a�� � � � � � � � �

a�� a��� � �

� � ����

���� � � � � � � � �

���

an���� � � an��n�� an��n�� �

an� an� � � � ann�� ann

�������� �

�������

a�� a�� � � � a�n�� a�n

� a�� a��� � � a�n

���� � � � � � � � �

������

� � � � � � an��n�� an��n� � � � � � � � ann

��������De�nizione ����� Matrice in forma di Hessemberg inferiore superiore �euna matrice per cui aij � �� j � i � � aij � �� j � i� ���������

a�� a�� � � � � �

a�� a��� � � � � �

������

� � � � � � � � � �

an���� � � an��n�� an��n�� an��n

an� an� � � � ann�� ann

�������� �

�������

a�� a�� � � � a�n�� a�n

a�� a�� a��� � � a�n

�� � � � � � � � �

������

� � � � � � an��n�� an��n� � � � � ann�� ann

��������De�nizione ����� Matrice trasposta �e una matrice che si ottiene scam�biando le righe con le colonne si indica con AT �

De�nizione ������ Matrice simmetrica �e una matrice ad elementi reali� percui A � AT �

De�nizione ������ Matrice de�nita positiva �e una matrice ad elementi rea�li� per cui xTAx � �� �x ����

De�nizione ������ Matrice inversa di una matrice quadrata A �e una ma�trice per cui A��A � AA�� � I�

De�nizione ����� Matrice ortogonale �e una matrice ad elementi reali� percui ATA � AAT � I AT � A���

De�nizione ������ Matrice a diagonale dominante �e una matrice per cui

jaiij �nX

j���j ��ijaijj � �i�

Page 22: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

����� Operazioni sulle matrici

Richiamiamo brevemente la de�nizione delle operazioni de�nite fra vettori ematrici�

De�nizione ����� Somma� di�erenza di due vettori x� y � Rn e� un vettorez� Rn i cui elementi sono dati da zi � xi yi

De�nizione ������ Somma� di�erenza di due matrici A�B � Rm�n e� unamatrice C � Rm�n i cui elementi sono dati da cij � aij bij�

De�nizione ������ Prodotto scalare fra due vettori x� y � Rn �e lo scalare

s � yTx �nPi��

xiyi� Il costo computazionale per ottenere s �e di n �ops per

�ops si intende un prodotto � una somma�

De�nizione ������ Prodotto fra una matrice A � Rm�n ed un vettore x� Rn

�e un vettore y � Rm i cui elementi sono dati da yi �nP

k��

aik xk� Il costo

computazionale per ottenere y �e di mn �ops�

De�nizione ����� Prodotto fra due matrici A � Rm�n e B � Rn�p �e una

matrice C � Rm�p i cui elementi sono dati da cij �nP

k��

aik bkj� Il costo

computazionale per ottenere C �e di mnp �ops�Il prodotto fra matrici quadrate non �e� in generale� commutativo AB ��

BA�

��� Metodi diretti �M�F��

Data una matrice quadrata A� dovendo risolvere il sistema lineare

Ax � b ����

�e noto� dai corsi di Analisi e Geometria� che la condizione det�A �� � �eequivalente all�esistenza della matrice inversa A�� di A per cui

x � A��b� ����

La determinazioine di x� soluzione di ����� nota la matrice inversa ha uncosto computazionale di n� �ops�

Page 23: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI DIRETTI �M�F�� ��

Prima di addentrarci nel calcolo della soluzione di ����� quando l�inversanon �e nota� �e istruttivo dare una chiara rappresentazione geometrica di cosavuol dire risolvere un sistema lineare� E� chiaro che una matrice A applicataad un vettore x lo trasforma in un vettore y� �e naturale chiedersi� tutti ivettori y sono immagine di qualche vettore x mediante la A� La risposta �e

s�i� se e solo se det�A �� � �o equivalentemente se �A���La ricerca della soluzione di un sistema lineare pu�o quindi essere vista

come la ricerca di una n�pla di numeri �le componenti del vettore x cheentrano nella combinazine lineare delle colonne della matrice A �Ai perrappresentare il termine noto �vettore b � ovvero

x�A� � x�A� � x�A� � ��� � xnAn � b�

L�esistenza di una soluzione di un sistema lineare �e quindi legata al fattoche il vettore b appartenga allo spazio generato dalle colonne della matriceA�

De�nizione ����� Spazio delle colonne range di A

R�A ���y � Rnjy � Ax� �x � Rn

���

Per l�unicit�a della soluzione �e importante la de�nizione di spazio nullo�nucleo di una matrice�

De�nizione ����� Spazio nullo null di A

N�A �� fx � RnjAx � �g ��E� facile far vedere che se ex �e soluzione di un sistema lineare e bx �� �

appartiene allo spazio nullo allora anche il vettore � ex� bx �e soluzione dellostesso sistema lineare �A ex � b� Abx � ���� A�ex� bx � b � � � b�

Esempio ����� Come esempio si considerino la matrice A ed il vettore b

A �

�� �� �

�� b �

���

�la matrice A �e singolare ed il suo spazio delle colonne coincide con i puntidella retta y� � �y�� per cui il sistema Ax � b non ammette soluzione essendob � R�A� Lo spazio nullo di A �e formato dai punti della retta x� � �x��infatti A

�x��x�

�� �� ovvio��

Page 24: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Se la matrice A �e diagonale con elementi diagonali tutti diversi da zero �eimmediato risolvere il corrispondente sistema lineare �disaccoppiato con ndivisioni�

Se la matrice A �e triangonale superiore �inferiore con elementi diagonalitutti diversi da zero �e possibile risolvere il sistema lineare associato con lasostituzione all�indietro �in avanti �back�forward�substitution con un costodell�ordine di n�

��ops�

Osservazione ����� Per �ops si intende una moltiplicazione divisione pi�uuna somma sottrazione in virgola mobile �oating�point� In algebra lineare�e naturale questa unit�a di misura in quanto il prodotto scalare� che �e alla basedel calcolo matriciale� costa proprio n prodotti ed n� � addizioni n �ops�

Si fa osservare che in entrambi i casi il sistema �e non singolare essendo ildeterminante diverso da zero �prodotto degli elementi diagonali�

Il problema �e quindi passare da un sistema con matrice piena ad unsistema equivalente con matrice diagonale o triangolare� Gauss forn�i perprimo un tale algoritmo combinando in modo opportuno le righe della matricee il termine noto�

����� Il metodo di eliminazione di Gauss

L�idea di Gauss si basa sulla semplice osservazione che il vettore soluzio�ne veri�ca tutte le equazioni del sistema lineare e quindi veri�ca anche un�equazione ottenuta combinando linearmente fra loro due di tali equazioni� E�quindi possibile passare da un sistema �pieno� ad uno equivalente triangolarein n�� passi eliminando al j�esimo passo l�incognita xj da tutte le equazionidalla �j���esima �no alla n�esima� Schematicamente scriviamo la ���� nellaforma� �����

a��x�� a��x�� ���� a�nxn � b�a��x�� a��x�� ���� a�nxn � b���� ��� ��� ��� ���

an�x�� an�x�� ���� annxn � bn

�� passso Si moltiplica la prima riga per��a��

a��

�e si sostituisce la seconda riga

con la somma delle prime due� si moltiplica la prima riga per��a��

a��

�e si sostituisce la terza riga con la somma fra la prima e la terza� ���� si

Page 25: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI DIRETTI �M�F�� �

moltiplica la prima riga per��an�

a��

�e si sostituisce la n�esima riga con

la somma fra la prima e la n�esima� giungendo a�����a��x�� a��x�� ���� a�nxn � b�

a���x�� ���� a��nxn � b����� ��� ��� ���

a�n�x�� ���� a�nnxn � b�n

�� passo Si opera analogamente sul sottosistema�a���x�� ���� a��nxn � b����� ��� ��� ���

a�n�x�� ���� a�nnxn � b�n

usando come moltiplicatori��a�

��

a���

����

��a�n�

a���

�si otterr�a

�����a���x�� a���x�� ���� a��nxn � b��

a���x�� ���� a��nxn � b����� ��� ��� ���

a�n�x�� ���� a�nnxn � b�n

j� passo Si opera analogamente sul sottosistema�aj��jj xj� ���� aj��jn xn � bj��j

��� ��� ��� ���

aj��nj xj� ���� aj��nn xn � bj��n

usando come moltiplicatori

��aj��j��j

aj��jj

����

��aj��nj

aj��jj

��

In n�� passi si arriva �se tutti gli aj��jj �� � ad un sistema equivalente

triangolare superiore� Il costo computazionale dell�algoritmo di Gauss �e n�

�ops� infatti sono necessarie �n���n�� �ops al primo passo� �n���n��al secondo� �n��n al terzo� �no a �� al passo �n���esimo� per cui in totale

flops �n��Xi��

i�i � � �n��Xi��

i� � �n��Xi��

i �

�n�n� ���n� �

�� �

n�n� �

��

n�

�� O�n��

Page 26: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Vale la pena di ricordare che la risoluzione di un sistema lineare mediantela regola di Cramer �calcolo di n� � determinanti utilizzando la de�nizionenel calcolo dei determinanti �det�A �

P���sa�s�a�s� ���ansn� �s�� s�� ���� sn

permutazione degli indici ��� �� ���� n � s numero di scambi fra �s�� s�� ���� sn e��� �� ���� n ha una complessit�a dell�ordine di �n � ��n� �n� �ops�

Su un elaboratore che esegue ���Mega�ops al secondo la risoluzione diun sistema lineare di �� equazioni in �� incognite richiede circa

��� ��� sec con l�algoritmo di Gauss�� ���� sec con la regola di Cramer

sfortunatamente ���� sec sono circa ���� anni ��Nell�algoritmo di Gauss abbiamo supposto che ad ogni passo fosse ben

de�nito il fattore moltiplicativo �aj��rj

aj��jj

� �r � j � �� ���� n il che equivale a

supporre aj��jj �� �� �j� Se aj��jj � �� per qualche j� l�algoritmo di Gaussnaturale si blocca� E� possibile modi�care l�algoritmo eseguendo la ricercadell�elemento di massimo modulo ��� � fra i restanti elementi della colonna�pivoting parziale e� se esiste� scambiare fra loro le righe j ed s �se s �e lariga contenente il massimo� Questa variante non richiede ulteriori �ops percui non aumenta la complessit�a computazionale� Va osservato che se tuttigli elementi della j�esima colonna sono nulli �partendo dalla j�esima riga allora il sistema non �e determinato �il determinante di A �e uguale a zero�La tecnica del pivoting parziale �e stabile computazionalmente per cui vienesempre e�ettuata la ricerca del massimo anche se aj��jj �� ��

����� Decomposizione LU

Il metodo di Gauss permette� se tutti i minori principali di testa della ma�trice A sono diversi da zero �il che �e equivalente a dire aj��jj �� �� �j� pro�varlo per esercizio� di decomporre la matrice A nel prodotto di due matrici�L triangolare inferiore �con � sulla diagonale principale ed U triangolaresuperiore

A � LU� ����

La matrice U coincide con la vecchia matrice triangolare superiore delsistema equivalente� mentre

L�� � Ln��Ln��Ln�����L�L�

Page 27: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI DIRETTI �M�F�� ��

�e ottenuta moltiplicando le matrici Lj del j�esimo passo del processo di Gauss

Lj �

������������

� � � � � � � � � � � �

�� � � �

���

�� � � �

� � ����

� �aj��j��j

aj��jj

� ����

�� � � � � � �

� � � � �aj��nj

aj��jj

� � � � �

��������������

Con il metodo di Gauss si ha quindi

Ln��Ln��Ln�����L�L�A � L��A � U ���

dalla ��� segue immediatamente la ����� Come si calcola la matrice L�Essendo� veri�carlo per esercizio� L��j � �I � Lj� si ha che la matrice L �

L��� L��� L��� ���L��n�� �e� senza eseguire alcuna operazione�

L �

������������

� � � � � � � � � � � �

�a��a��

� � � ����

�a��a��

� � � �� � �

������ �

aj��jj��

aj��j��j��

�aj��j��j

aj��jj

� ����

�an���a��

���� � � � � � �

�an�a��

� � � �aj��nj

aj��jj

� � � �an��nn��

an��nn�

�����������������

�la ���� si ottiene semplicemente tenendo conto della particolare strutturadelle matrici L��j �

La seguente matrice

A �

�� �� �

�non ammette decomposizione LU �a�� � �� primo minore principale di testanullo� �e facile osservare che� se si scambiano la prima e la seconda riga� la

matrice A diviene la matrice U �

�� �� �

�� per cui la decomposizione �e

implicita �L � I� Possiamo quindi formulare il seguente

Page 28: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Teorema ����� Data una matrice A� se det�A �� �� esiste una matrice Pdi permutazione per cui

PA � LU�

Dimostrazione ����� La dimostrazione del ����� �e costruttiva coincidecon il metodo di Gauss con pivoting parziale��

L�importanza della decomposizione PA � LU nasce dal fatto che non�e necessario conoscere il vettore dei termini noti quando si esegue la fasedi decomposizione� La risoluzione di un sistema lineare pu�o quindi essereschematizzata nel modo seguente

�� No pivoting

Ax � bA � LU

��� LUx � b ��

�Ly � bUx � y

�� pivoting

Ax � bPA � LU

��� LUx � Pb ��

�Ly � PbUx � y

Si fa osservare che la parte computazionalmente onerosa �e la decompo�sizione LU �n� �ops� mentre la risoluzione dei � sistemi triangolari� partedestra dello schema� costa solo n� �ops� La fase di determinazione del vettorePb� nel secondo schema� costa zero �ops �sono solo scambi di righe��

Anche la costruzione della matrice P �e immediata e pu�o essere fattautilizzando solo un vettore di n interi� Precisamente� de�nito

sp � ��� �� �� ���� n� �� n ����

il vettore dei primi n interi �sp�k � k� per costruire la matrice P bastaspostare� al j�esimo passo dell�eliminazione di Gauss� il contenuto della cella jcon quello della cella r del vettore sp �dove r �e l�indice di riga del pivot al passoj� Il vettore de�nito in ���� rappresenta� in modo compatto� la matriceidentit�a �� in posizione �k� sp�k� k � �� �� ���n� � altrove� analogamentela matrice P avr�a � ovunque e � solo nelle posizioni �k� sp�k� k � �� �� ���n�

Page 29: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI DIRETTI �M�F�� �

����� Decomposizione di Cholesky

Se la matrice A �e simmetrica e de�nita positiva �e possibile dimostrare l�esi�stenza della seguente decomposizione �decomposizione di Cholesky

A � V TV ����

dove V �e una matrice triangolare superiore con elementi diagonali positivi�La dimostrazione dell�esistenza di tale decomposizione �e costruttiva� e

viene qui riportata perch�e evidenzia come sia proprio la decomposizione diCholesky la via per determinare se una matrice data A �e de�nita positiva�Essendo

V T �

������ v��v�� v�����

���� � �

vn��� vn��� � � � vn��n��vn� vn� � � � vnn�� vnn

������� �

V �

������ v�� v�� � � � vn��� vn�

v�� � � � vn��� vn�� � �

������

vn��n�� vnn��vnn

������� �dalla ���� si ha� eseguendo il prodotto righe per colonne� partendo dallaprima riga

v��v�� � v��� � a�� � v�� �pa��

v��vj� � a�j � aj� � vj� �aj�v��

� j � �� �� ���� n

per la seconda riga� vj� sono ora noti �j � �� si ha

v��� � v��� � a�� � v�� �qa�� � v���

v��vj� � v��vj� � a�j � aj� � vj� �aj� � v��vj�

v��� j � �� �� ���� n

in generale al passo k�esimo si ottiene

vkk �

vuutakk �k��Xm��

v�mk� k � �� ���� n ���

Page 30: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

vki �

aki �k��Pm��

vkmvmi

vkk� k � �� ���� n� �� i � k � �� ���� n �����

nelle �������� la sommatoria va ignorata se il valore d�arrivo �e minore delvalore di partenza�

Il costo computazionale della decomposizioni di Cholesky �e n�

��ops �oltre

a n estrazioni di radice quadrata� Proprio la necessit�a di dover estrarre laradice quadrata permette di veri�care se la matrice di partenza era de�nitapositiva� infatti se un radicando risulta negativo l�algoritmo si interrompe�non esiste la decomposizione per cui la matrice di partenza non era de�nitapositiva� E� facile dimostrare che� ad ogni passo k� il radicando �e dato dalquoziente fra il minore principale di testa� della matrice di partenza� di ordinek e quello di ordine k � ��

Si pu�o osservare che nelle �������� i prodotti scalari indicati possonoessere eseguiti in doppia precisione� pur lasciando le variabili in sempli�ce� aumentando cos�i l�accuratezza del risultato �algoritmo stabile oltre chee�ciente�

E� evidente che una volta nota la decomposizione di Cholesky il sistemalineare

Ax � b

pu�o essere risolto mediante il seguente schema

Ax � bA � V TV

��� V TV x � b ��

�V Ty � bV x � y

�� Analisi dell�errore �M�F��

La soluzione del sistema lineare Ax � b calcolata con uno dei metodi proposti�chiamiamola xc di�erir�a� per l�inevitabile presenza di errori �cfr� capitolo�� dalla soluzione vera �chiamiamola xv� quanto grande sar�a questa di�eren�za� Dipender�a solo dal metodo usato� Dipender�a dalla matrice del sistemalineare�

Prima di rispondere a queste domande �e necessario de�nire una misuradi distanza fra vettori e matrici�

Page 31: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� ANALISI DELL�ERRORE �M�F�� ��

����� Richiami sulle norme di vettore

De�nizione ����� Dato un vettore x si de�nisce norma del vettore x kxkuna funzione a valori non negativi che gode delle seguenti propriet�a

�� kxk � � �x � Rn� kxk � � � x � �

�� kxk � jj kxk �x � Rn e � � R

�� kx � yk � kxk� kyk �x� y � Rn

Tre classici esempi di norma di vettore sono

�� kxk� � jx�j� jx�j� ��� � jxnj�� kxk� �

px�� � x�� � ��� � x�n

�� kxk � maxifjxijg

Analogamente per le matrici

De�nizione ����� Data una matrice A si de�nisce norma della matrice AkAk una funzione a valori non negativi che gode delle seguenti propriet�a

�� kAk � � �A � Rn�n� kAk � � � A � �

�� kAk � jj kAk �A � Rn�n e � � R

�� kA � Bk � kAk� kBk �A�B � Rn�n

�� kABk � kAkkBk �A�B � Rn�n

Norme di matrice possono essere ricavate dalle norme di vettore �normaindotta o naturale ponendo

kAk � maxx

kAxkkxk �A � Rn�n� x � Rn�

norme di vettore e norme di matrice si dicono compatibili se

kAxk � kAkkxk�Esempi di norme di matrice sono

Page 32: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

�� kAk� � maxj

�nPi��

jaijj�

�� kAk� �p��ATA �dove ��A �e il raggio spettrale di A�

�� kAk � maxi

�nP

j��

jaijj�

�� kAkE �

snP

j��

nPi��

a�ij �questa norma �e compatibile con la norma � di

vettore ma non �e indotta�

����� Stima dellerrore

Dato il sistema lineare Ax � b� detta xc la soluzione calcolata si pu�o deter�minare il residuo

r � Axc � b�

essendo xv � A��b� posto e � xc � xv� si ha

r � Axc � b � Axc � Axv � A�xc � xv � Ae

per cuie � A��r

da cui seguekek � kxc � xvk � kA��rk � kA��kkrk �����

la ����� non �e molto signi�cativa in quanto anche se krk �e piccola l�errorepu�o essere grande se kA��k �e grande� Quello che si pu�o dire �e che se krk �egrande allora la soluzione calcolata non �e �buona��

Possiamo dimostrare il seguente

Teorema ����� Dato il sistema Ax � b� con A matrice non singolare� ed ilsistema A�x � �x � b � �b� si ha

K�A

k�bkkbk � k�xk

kxk � K�Ak�bkkbk �����

dove K�A � kAkkA��k �e il numero di condizionamento della matrice A�

Page 33: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� ANALISI DELL�ERRORE �M�F�� ��

Dimostrazione ����� Valgono le seguenti � uguaglianze

�� Ax � b �� �� x � A��b�� A�x � �b �� �� �x � A���b

per cui� usando la � e la �� si ha

kxk � kA��bk � kA��kkbkk�bk � kA�xk � kAkk�xk

��� kxkk�bk � K�Akbkk�xk

da cui�

K�A

k�bkkbk � k�xk

kxk �

Analogamente� usando la � e la �� si ha

kbk � kAxk � kAkkxkk�xk � kA���bk � kA��kk�bk

��� k�xkkbk � K�Ak�bkkxk

da cuik�xkkxk � K�A

k�bkkbk

che prova la ����� �

Si pu�o dimostrare un risultato pi�u generale se si considera una variazioneanche sugli elementi della matrice�

Osservazione ����� Valgono le seguenti relazioni

�� K�A � K�A���

�� K�A dipende dalla norma scelta�

�� K�A � �� infatti � � kIk � kAA��k � kAkkA��k � K�A��

�� Le matrici ortogonali �QTQ � QQT � I hanno K��A � �� Infattiessendo kAk� �

p��ATA si ha

kQk� �p��QTQ �

p��I � ���

Page 34: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

� Se A �e simmetrica �A � AT allora K��A �����max�A �min�A

��� � Infatti si pu�o

facilmente dimostrare che kAk� � j�max�Aj

kAk� �p��ATA �

p��A� �

q��max � j�max�Aj

per cui� essendo ��A�� � ���A

kA��k� ��

j�min�Aj ��

�� Un classico esempio di matrice mal condizionata �e la matrice di Hilbertdi ordine n

Hn �

������ � �

�� � � �

n���n

��

��

� � � �n

�n��

������ � � � ���

����

n���n

� � � ��n��

��n��

�n

�n��

� � � ��n��

��n��

������� �

�� Metodi iterativi �M�F��

Quando la matrice A del sistema lineare ���� risulta sparsa e non strutturata�quindi con un numero di elementi diversi da zero dell�ordine di n e nondisposti a banda o con altra particolare struttura pu�o essere convenientel�uso di metodi iterativi� Il vantaggio nasce dal fatto che i metodi iterativi�sfruttando la sparsit�a della matrice� eseguono il prodotto matrice vettorein meno di n� �ops� per cui possono essere assai e�cienti� Per contro imetodi diretti� applicati a matrici sparse non strutturate� generano il �ll�in�riempimento della matrice durante la fase di eliminazione� per cui hannouna complessit�a computazionale che �e sempre di n� �ops� Poich�e la soluzionedel sistema di partenza �e ottenuta come limite di una successione

xk�� � Bxk � g �����

dove la matrice d�iterazione B �e scelta in modo che il sistema

x � Bx � g �����

sia equivalente al sistema ����� bisogner�a garantire non solo la convergenzadella ����� alla soluzione x� ma anche essere in grado di stimare l�errore chesi commette arrestandosi ad una data iterazione�

Page 35: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI ITERATIVI �M�F�� �

����� Ranamento iterativo di Wilkinson

Un primo metodo� proposto da Wilkinson� viene utilizzato anche quando lamatrice A non �e sparsa per ottenere una miglior stima della soluzione delsistema ���� calcolata con un metodo diretto�

Sia xc la soluzione del sistema ���� calcolata� per esempio� mediantel�algoritmo di Gauss� calcoliamo il residuo

r � Axc � b ����

in doppia precisione �l�uso di una precisione superiore nel calcolo del residuo�e fondamentale per non rischiare di ottenere un residuo identicamente nullo�sia ec la soluzione del sistema

Ae � r

calcolata utilizzando la stessa decomposizione LU della matrice del sistema�Si pu�o osservare che� se i calcoli sono eseguiti �senza errori� allora

xv � xc � ec � xc � A��r � xc � A���Axc � b � A��b

di fatto� a causa degli inevitabili errori dovuti alla �nitezza della macchina�si ha

kxvera � xck � kxvera � xvked� in generale� sono su�cienti � o � iterazioni per raggiungere la soluzionenella precisione di macchina�

Se keck non decresce vuol dire che la matrice del sistema �e mal condizio�nata� Se �e nota la decomposizione LU di A il metodo non �e molto costoso�bastano solo n� �ops per iterazione�

����� Metodo di Jacoby

Si trasforma il sistema ���� nella forma equivalente ����� utilizzando ilseguente splitting �partizionatura della matrice A � L � D � U � dove L �etriangolare inferiore con zeri sulla diagonale principale� D �e diagonale e U �etriangolare superiore con zeri sulla diagonale principale� La matrice B ed ilvettore g sono quindi dati da

B � �D���L � U� g � D��b�

Page 36: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Il metodo diventa quindi

xk�� � �D���L � Uxk � D��b� �����

Per l�equivalenza con il sistema ����� detta xvera la soluzione del sistema�si avr�a� dalla �����

xvera � �D���L � Uxvera � D��b

o� equivalentemente�xvera � �I �B��g�

La convergenza del metodo di Jacoby� ed in generale di tutti i metodidella forma ������ dipender�a quindi dalle propriet�a di B� pi�u precisamentevale il seguente

Teorema ���� Condizione necessaria e su�ciente per la convergenza delmetodo ���� �e che ��B � ��

Dimostrazione ����� Preso un vettore iniziale x� arbitrario la ���� cifornisce

x� � Bx� � g �����

x� � Bx� � g � B�x� � �B � Ig

x� � Bx� � g � B�x� � �B� � B � Ig

� � � � � �xk � Bxk�� � g � Bkx� � �Bk�� � � � �� B � Ig

vale il seguente lemma la cui dimostrazione �e lasciate come esercizio�

Lemma ����� Se e solo se ��B � � allora

limk

Bk � �

limk

�kXi��

Bi

�� I � B � � � � � �I � B����

per cui� dalle ��� � utilizzando il lemma ����� si ha

k ��� xvera � �I � B��g��

Page 37: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI ITERATIVI �M�F�� ��

Osservazione ���� Ricordando che� data una matrice A� per ogni de�ni�zione di norma si ha

j��Aj � kAkvale il seguente

Teorema ���� Condizione su�ciente per la convergenza della sucessione���� �e che esista una norma di B per cui kBk � ���

Oltre a garantire la convergenza del metodo �e importante poter stimarel�errore che si commette fermandosi alla k�esima delle iterazioni ����� �erroredi troncamento� Vale il seguente risultato

Teorema ��� Se kBk � � allora per l�errore alla k�esima iterazione si ha

kxvera � xkk � kBkk kxvera � x�k

kxvera � xkk � kBkk�� kBk kx� � x�k

Dimostrazione ����� Per la prima disuguaglianza si ha

xk � Bxk�� � g

xvera � Bxvera � g

sottraendoxvera � xk � B�xvera � xk��

e quindixvera � xk � Bk�xvera � x�

passando alle norme

kxvera � xkk ���Bk�xvera � x�

�� � ��Bk�� kxvera � x�k � �����

Per la seconda disuguaglianza� osserviamo che

kxvera � x�k � kxvera � x� � x� � x�k � kxvera � x�k� kx� � x�k� kBk kxvera � x�k� kx� � x�k

essendo kBk � �� si ha

��� kBk kxvera � x�k � kx� � x�k ����

kxvera � x�k � �

�� kBk kx� � x�k

Page 38: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

dalle ��������� si ha

kxvera � xkk �kBkk

�� kBk kx� � x�k �� �����

Osservazione ���� Nella ���� la quantit�a kx� � x�k �e nota dopo il primopasso per cui si pu�o calcolare quanti passi saranno necessari per ottenere unapre�ssata accuratezza essendo kBk � ��

����� Metodo di Gauss�Seidel

Sia dato il sistemaAx � b

e la matrice A sia partizionata come

A � L � D � U

la j�esima componente del vettore xk �xjk fornita dal metodo di Jacoby �edata da

Djjxjk � �Ljxk�� � Ujxk�� � bj� k � �� �� � � � �����

dove Lj e Uj sono le righe j�esime delle matrici L e U � Essendo le componentidel vettore Lj

Lji � �� i � j

nella ����� al generico passo k�esimo� il primo vettore xk�� pu�o essere so�stituito dal vettore xk in quanto sono gi�a state calcolate le prime �j�� com�ponenti del vettore xk stesso� In questo modo si tiene subito conto dellevariazioni apportate nell�iterazione� Si osservi che le restanti n� j � � com�ponenti� non ancora note� vengono moltiplicate per zero� Si ottiene quindi ilmetodo �di Gauss�Seidel

Djjxjk � �Ljxk � Ujxk�� � bj� k � �� �� � � � � n

che� con ovvio signi�cato dei simboli� pu�o essere scritto in forma matriciale

Dxk � �Lxk � Uxk�� � b�

La matrice d�iterazione del metodo ed il vettore g sono quindi

B � ��D � L��U � g � �D � L��b

Page 39: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI ITERATIVI �M�F�� �

per cui sono applicabili� per la convergenza del metodo� i teoremi �������������Si pu�o osservare� analizzando la struttura dei due metodi� che se entrambii metodi convergono� in generale� il metodo di Gauss�Seidel converger�a pi�urapidamente del metodo di Jacoby�

Diamo ora due teoremi di convergenza

Teorema ���� Se la matrice A del sistema lineare ��� �e fortemente diago�nalizzata allora i metodi di Jacoby e Gauss�Seidel convergono per ogni sceltadel vettore iniziale�

Dimostrazione ����� La dimostrazione� per il metodo di Jacoby� �e imme�diata ricordando che� essendo A fortemente diagonalizzata� allora

jaiij �nX

j��j ��i

jaijj � �i � �� �� � � � � n

per cui kBk � ���

Teorema ��� Se la matrice A del sistema lineare ��� �e simmetrica ede�nita positiva allora il metodo di Gauss�Seidel converge per ogni scelta delvettore iniziale��

����� Metodi di rilassamento SOR

Si basano sulle seguenti due osservazioni�

�� la velocit�a di convergenza di un metodo iterativo della forma �����dipende dal raggio spettrale della matrice d�iterazione B� pi�u preci�samente pi�u ��B �e prossimo a zero tanto pi�u velocemente convergeil metodo ����� �se ��B � � il metodo converge al pi�u in n passi�diventa un metodo diretto�

�� se la matrice B dipende da un parametro �B � B� si pu�o pensaredi trovare un valore ottimale di � opt per cui il raggio spettrale risultiminimo ���B� opt � ��B� � ��

Partendo dal metodo di Gauss�Seidel� scritto nella forma

xk�� � D�� ��Lxk�� � Uxk � b�

Page 40: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

sottraendo ad entrambi i membri xk si ha

xk�� � xk � D�� ��Lxk�� � Uxk � b�� xk �����

de�nito il vettore vk

vk � D�� ��Lxk�� � Uxk � b�� xk

la ����� si pu�o scrivere� introducendo il parametro � �� nella forma

xk�� � xk � vk� �����

Nella ����� il vettore variato �e ottenuto dal vettore precedente muoven�dosi nella direzione di vk con un passo pari a kvkk� l�importante sar�a trovareil passo in modo che xk�� sia pi�u vicino a xvera di quanto non lo sia xk� Informa matriciale la ����� �e equivalente a

xk�� � xk � D�� ��Lxk�� � �D � Uxk � b�

ovvero�D � Lxk�� � ���� D � U xk � b

e in de�nitiva

xk�� � �D � L�� ���� D � U xk � b� �����

Nella ����� la matrice d�iterazione �e quindi

B� � �D � L�� ���� D � U

che si riduce� per � �� alla matrice del metodo di Gauss�Seidel� mentreper la convergenza del metodo �e necessario �non su�ciente che � � � ��Se �e costante ad ogni passo si ha un metodo stazionario� se viene variatoad ogni iterazione si ha un metodo non stazionario �si cambia il passo� Laricerca del valore ottimale di dipende da problema a problema e non pu�oessere qui trattata�

����� Metodi non stazionari

L�idea alla base dei metodi non stazionari �e quella di cercare ad ogni iterazionedi muoversi lungo un�opportuna direzione con un opportuno passo in modo

Page 41: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI ITERATIVI �M�F�� ��

che l�iterata �k���esima sia �pi�u vicina� alla soluzione dell�iterata k�esima�Dedurremo i metodi nel caso di matrici simmetriche e de�nite positive�

De�niamo l�errore alla generica iterazione k

ek � xvera � xk

ed associamogli la seguente forma quadratica de�nita positiva

��ek � �Aek� ek � eTkAek � �xvera � xkTA�xvera � xk ����

� xTveraAxvera � �xTveraAxk � xTkAxk

in modo analogo si de�nisce ��ek��� Cerchiamo� analogamente a quantofatto in ������

xk�� � xk � vk� �����

Utilizzando le ���� ����� si pu�o cercare e vk in modo che

��ek�� � ��ek�

e la succesione ����� converga� per k ���� a xvera�Per la determinazione di vale il seguente

Teorema ���� Se A �e simmetrica e de�nita positiva

��ek�� � ��ek�

per ogni dato da

� �rTk vkvTkAvk

� � � � � �� �����

dove rTk �e il residuo �����

Dimostrazione ����� Dalle ��������� si ha

��ek�� � �xvera � xk��TA�xvera � xk��

� xTveraAxvera � �xTveraAxk�� � xTk��Axk��

� xTveraAxvera � �xTveraA�xk � vk � �xk � vkTA�xk � vk

� ��ek� �xTveraAvk � �xTkAvk � �vTkAvk� ��ek� �eTkAvk � �vTkAvk�

Page 42: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Perch�e sia��ek�� � ��ek�

deve essere�vTkAvk � �eTkAvk � � �����

ovvero� essendo Aek � rk�

� �eTkAvkvTkAvk

� �rTk vkvTkAvk

��

Osservazione ���� La condizione su � �e la stessa data su nei metodiSOR� E� necessaria per la convergenza� ma non su�ciente bisogna ancorascegliere la direzione lungo cui muoversi� Il valore di � � � �e quello che o�re�in generale� il miglior guadagno essendo corrispondente al valore minimodella �����

Metodo delle coordinate

Il vettore vk che compare nella ����� viene scelto �metodo delle coordinateciclico ciclicamente uguale al versore ej

vk � ej� k � j modn

mentre � dato dalla ������ risulta uguale a

� �rjkajj

����

dove rjk �e la componente j�esima di rk� Con queste scelte �e garantita laconvergenza del metodo�

Se il vettore vk che compare nella ����� viene scelto uguale al versoreej� dove j �e l�indice della massima componente in modulo di rk� ed �edato dalla ����� questo �e equivalente ad annullare� al passo k�esimo� lamassima delle componenti del residuo� si ha il metodo delle coordinate NONciclico� La convergenza NON �e pi�u garantita �potrebbero esserci componentidi xk che non vengono mai modi�cate per ripristinare la convergenza bisognagarantire che dopo m passi �ovviamente m � n tutte le componenti di xksiano state modi�cate� Se ci�o non avviene si devono modi�care ciclicamentele componenti non variate�

Page 43: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� SISTEMI SOVRADETERMINATI �M�F�� ��

Metodo del gradiente

Si sceglie come direzione vk che compare nella ����� il vettore residuo rk cherisulta essere il gradiente della forma quadratica ���� �da qui il nome delmetodo� Lo scalare � dato dalla ������ risulter�a quindi uguale a

� �rTk rkrTkArk

Vale la pena di osservare che la velocit�a di convergenza del metodo delgradiente dipende dal valore

K�A� �

K�A � �

per cui la convergenza pu�o essere molto lenta se A �e mal condizionata�

��� Sistemi sovradeterminati �M�F��

In svariati problemi si ha a che fare con modelli lineari in cui �e importanteindividuare dei parametri �vettore delle incognite x� potendo disporre diun numero di osservazioni �vettore dei termini noti b superiore al numerodei parametri� per cui la matrice del sistema lineare �A risulta rettangolare�sistemi sovradeterminati� Le osservazioni sono sovente a�ette da errori�errori di misura per cui non esiste una soluzione� in senso classico� delproblema �il vettore termini noti non appartiene allo spazio delle colonne diA�

Dato il sistema lineareAx � b

con A � Rmn� x � Rn� e b � Rm� � m � n� si cerca il vettore x che veri�ca

minxkAx� bk�� �����

� x �e soluzione nel senso dei minimi quadrati�Vale il seguente

Teorema ����� Se �y�� con y� � A x�� tale per cui

�b� y�Ty � �� �y � R�A�

il vettore �b�y� risulta ortogonale a tutti i vettori di R�A allora x� �e unasoluzione di �����

Page 44: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Dimostrazione ����� Utilizzeremo la tecnica del completamento del qua�drato

� � ��b� y����

� �b� yT �b� y � bT b� �bTy � yTy

aggiungendo ��b� y�Ty � �� per ipotesi� si ha

� � bT b� �bTy � yTy � ��b� y�Ty

� bT b � �bTy � �y�T y � �bTy � yTy

� bT b� �y�Ty � yTy

aggiungendo e togliendo y�Ty�

� � bT b� y�T y� � �y�Ty� � �y�T y � yTy

� kbk�� ���y����

�� �y� � yT �y� � y

in conclusione

� � ��b� y����

� kbk�� ���y����

�� �y� � yT �y� � y

che risulta minima quando y � y���

Osservazione ����� Risolvere un sistema lineare nel senso dei minimi qua�drati �e quindi equivalente a trovare un vettore x la cui immagine� mediante A��e il vettore y proiezione ortogonale del vettore termine noto b sullo spaziodelle colonne di A� Se A ha rango massimo tale soluzione sar�a unica�

Per determinare il vettore x� soluzione� osserviamo che� in virt�u delteorema ��� �� e della propriet�a commutativa del prodotto scalare� si ha

�b� y�Ty � yT �b� y� � �� �y � R�A�

ed� essendoy� � Ax�� y � Ax

si ha

yT �b� y� � xTAT �b� Ax� � �� �x � Rn

� � xT �AT b� ATAx�

da cuiAT b� ATAx� � �

Page 45: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� SISTEMI SOVRADETERMINATI �M�F�� �

e quindiATAx� � AT b �����

che prende il nome di equazione normale� Se A �e a rango massimo alloraATA �e invertibile� per cui

A� � �ATA��AT

ex� � A�b

�e l�unica soluzione �la matrice A� prende il nome di pseudo inversa di Moore�Penrose�

Se A �e mal condizionata peggio sar�a il condizionamento di ATA� percui la risoluzione di ����� pu�o risultare poco accurata� Facciamo un brevecenno ad una tecnica risolutiva� utile quando A �e a rango massimo ma malcondizionata� che si basa sul seguente teorema

Teorema ����� Data una matrice A � Rmn� a rango massimo� esistonodue matrici� Q � Rmm ortogonale e R � Rmn pseudo�triangolare superioreil pre�sso pseudo per ribadire che �e rettangolare� con elementi rii �� �� taliper cui

A � QR�� �����

�Forniremo un algoritmo per la determinazione della decomposizione QRnel paragrafo ������ relativo al calcolo degli autovalori�

Utilizzando il teorema ��� �� si ha

minxkAx� bk�� � min

xkQRx� bk�� � min

x

��Q�Rx�QT b����

� minx

��Rx�QT b����

� minx

���Rx�bb������

Scriviamo la matrice R ed il vettore bb � QT b� nella forma

R �

������ r�� � � � r�n

�� � �

������

� � � rnn� � � � �� � � � �

������� �

�R��

�� bb �

�������

bb����bbn���bbm

�������� �

� bb�bb�

Page 46: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

dove R�� � Rnn� bb� � Rn e bb� � Rm�n� �e su�ciente quindi risolvere il sistema�quadrato non singolare

R��x � bb�mantenendo il condizionamento in norma � uguale �K��A � K��R �K��R���

Se A non �e a rango massimo allora ATA �e singolare� si pu�o ancoraparlare di soluzione nel senso dei minimi quadrati e di pseudo�inversa� mal�argomento va oltre i con�ni di questo corso�

�� Riepilogo �M�F��

Sintetizziamo quanto esposto relativamente alla risoluzione dei sistemi lineariabbiamo

�� Matrici piene �elementi diversi da zero dell�ordine di n�

�a NON simmetriche� metodo di Gauss e sue varianti �A�LU� PA�LUcon costo di n�

��ops�

�b matrici sparse e strutturate� Gauss� A�LU� PA�LU con costodipendente dalla ampiezza della banda�

�c simmetriche de�nite positive� decomposizione di Choleshy concosto di n�

��ops�

�� Matrisci sparse �numero di elementi diversi da zero �� n�

�a NON de�nite positive� Jacoby� Gauss�Seidel e SOR� Convergenzae velocit�a di convergenza dipendono da ��B�

�b simmetriche e de�nite positive� Coordinate e Gradiente� Velocit�adi convergenza dipende da K�A�

�� Accuratezza della soluzione dipende dal numero di condizionamento�K�A � kAk kA��k �

�� Ra�namento iterativo di Wilkinson per migliorare la soluzione� quandoK�A non �e troppo grande�

� Soluzione nel senso dei minimi quadrati per sistemi sovradeterminati�Decomposizione QR�

Page 47: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� RIEPILOGO �M�F�� ��

����� Esercizi

Esercizio ���� Calcolare� se esiste� la decomposizione LU della seguentematrice

A �

��� � � � �� � � �� � � �� � � �

���� �Giusti�care la risposta�

Esercizio ���� Risolvere� mediante il metodo di Jacoby� il seguente sistemalineare Ax � b� dove

A �

����� � � � �

� � � � �

������ � b �

����� ����

������usando diversi valori per la tolleranza richiesta nel test di arresto e partendosempre da x� � b� Giusti�care i risultati ottenuti�

Esercizio ���� Se una matrice �e triangolare allora �e ben condizionata �Giusti�care la risposta proponendo anche un esempio�

Esercizio ��� Sia A una matrice simmetrica e fortemente diagonalizzatadimostrare che pu�o essere decomposta nel prodotto

A � LDLT

dove L �e triangolare inferiore con � sulla diagonale principale e D �e diago�nale�

Esercizio ���� Dati i punti del piano xy� le cui coordinate sono date nellaseguente tabella

xi ��� �� ��� ��� �yi ���� ���� � ���� ��

calcolare la retta di regressione�

Page 48: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

��� Autovalori di matrici �M�F��

�� �� Richiami e de�nizioni

Data una matrice quadrata A si de�niscono autovalori e autovettori rispet�tivamente quei valori reali o complessi e quei vettori non nulli che veri�canol�equazione

Au � �u� �����

Il sistema ����� �e equivalente al sistema omogeneo

�A� �Iu � �

che ammette soluzioni non banali �u �� � se e solo se � �e tale per cui

det�A� �I � pn�� � �

dove pn�� �e un polinomio di grado n�

Teorema ����� Se la matrice A �e diagonalizzabile allora ammette n auto�vettori linearmente indipendenti�

Dimostrazione ����� Se A �e diagonalizzabile allora

Au � C��DCu � �u

DCu � �Cu

posto

Cu � w

segue

Dw � �w

ed essendo D diagonale �i � Dii sono autovalori e w � ei sono i corri�spondenti autovettori� per cui� essendo u � C��w anche gli u risulterannolinearmente indipendenti��

Teorema ����� Se la matrice B �e trasformata per similitudine della matriceA B � C��AC� allora ��A � ��B e u�A � Cu�B�

Page 49: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� � AUTOVALORI DI MATRICI �M�F�� �

Dimostrazione ����� Essendo B � C��AC� da Bu � �u segue

C��ACu � �u

ACu � �Cu

Aw � �w��

Teorema ����� Se la matrice A ha autovalori ��A �� �� allora la matri�ce A�� ha autovalori ��A�� � �

��A e gli autovettori sono uguali u�A �

u�A����

Teorema ���� Se la matrice A ha autovalori ��A� allora la matrice �A�I ha autovalori ��A�I � ��A� e gli autovettori sono uguali u�A �u�A� I��

Si lascia al lettore� come esercizio� la dimostrazione di questi due ultimiteoremi�

Per il calcolo degli autovalori ed autovettori considereremo due famigliedi metodi� i metodi locali� che permettono il calcolo di un particolare auto�valore e del corrispondente autovettore� ed i metodi globali che� basandosisu trasformazioni per similitudine� permettono di stimare simultaneamentetutti gli autovalori ed autovettori�

Vale la pena di notare che� se si conosce un autovettore della matrice A��e immediato il calcolo del corrispondente autovalore� vale infatti �quozientedi Rayleigh

Ax � �x

xTAx � �xTx

� �xTAx

xTx��

al contrario� noto � il calcolo di x mediante la de�nizione ����� risulta dif��coltoso� se non impossibile� in quanto anche un piccolo errore su � rende ilsistema ����� non pi�u omogeneo� per cui solo x � � veri�ca la ������

�� �� Metodi locali

Fra i metodi locali presenteremo il metodo delle potenze� che permette ilcalcolo dell�autovalore di modulo massimo� ed il metodo delle potenze inverse�

Page 50: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Metodo delle potenze

Per sempli�care la dimostrazione della convergenza del metodo delle potenzefaremo l�ipotesi che

�� la matrice A abbia n autovettori linearmente indipendenti�

�� che esista un solo autovalore di modulo massimo ���

�� il vettore iniziale x�� della successione xk�� � Axk� dipenda da u�autovettore associato a ���

Vale il seguente

Teorema ����� Data una matrice A con un solo autovalore di modulo mas�simo

j��j � j��j � j��j � � � � � j�nje n autovettori ui linearmente indipendenti� posto

x� � �u� � �u� � � � �� nun �nXi��

iui

con � �� �� la successionexk�� � Axk �����

converge ad u� mentre

�jk �xjk��

xjk� xjk �� �

converge a ���

Dimostrazione ����� Dalla ����� per la de�nizione di autovalore e au�tovettore� si ha

x� � Ax� �nXi��

iAui �nXi��

i�iui

x� �nXi��

i�iAui �nXi��

i��iui

� � �xk�� �

nXi��

i�kiAui �

nXi��

i�k��i ui

Page 51: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� � AUTOVALORI DI MATRICI �M�F�� �

ovvero� raccogliendo �k���

xk�� � �k���

��u� �

�����

�k��

�u� � � � ����n��

�k��

nun

�����

dalla ���� sia halimk

xk�� � u�

si ricordi che gli autovettori sono de�niti a meno di una costante�La ����� scritta per k� diviene

xk � �k�

��u� �

�����

�k

�u� � � � ����n��

�k

nun

facendo il quoziente fra le j�esime componenti non nulle�

�jk �xjk��

xjk� ��

�uj� �

�����

�k��

�uj� � � � �

�uj� �

�����

�k�u

j� � � � �

per cuilimk

�jk � ����

Dalla ���� si vede che la velocit�a di convergenza a u� e �� dipende dallavelocit�a con cui il quoziente ��

��va a zero �convergenza lineare�

Il fattore �k pu�o essere calcolato in generale nella forma

�k �vTk xk��

vTk xk

dove vk �e un generico vettore �la scelta precedente equivale a vk � ej� Se sipone vk � xk� si ottiene il quoziente di Rayleigh

�k �xTkAxkxTk xk

�xTk xk��

xTk xk

che garantisce una convergenza quadratica� dipendente da�����

��

� se la ma�

trice A �e simmetrica�

Page 52: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Vediamo ora se si possono rilassare� o modi�care� le ipotesi fatte nelprecedente teorema� La lineare indipendenza degli autovettori pu�o esseresostituita con la diagonalizzabilit�a di A o dall�ipotesi che gli autovalori sianodistinti �perch�e �� Il vincolo � �� �� se rimosso fa si che inizialmente siabbia convergenza su �� �ovviamente se � �� �� ma con il crescere delleiterazioni� a causa degli inevitabili errori di arrotondamento� si converge an�cora a ��� Se esistono due autovalori di modulo massimo reali il metodo dellepotenze funziona ancora pur di applicarlo alla matrice A� �basta ricordareche ��A� � ��A��

Operativamente� nella pratica� converr�a normalizzare i vettori xk al �nedi evitare l�over�ow�

Un semplice programma in MATLAB per il metodo delle potenze �e ilseguente

x�rand�max�size�A�� ���

x�x�norm�x ���

betaold���

beta���

for i�������

if abs��beta�betaold��beta� � eps�e��

break

end

betaold�beta�

y�Ax�

beta�y�x�

x�y�norm�y ���

end

Metodo delle potenze inverse

E� una variante del metodo delle potenze che si basa sull�osservazione che� seuna matrice �e invertibile� allora

��A�� ��

��A

ed i corrispondenti autovettori coincidono�Si pu�o quindi calcolare l�autovalore di modulo minimo di A� ed il corri�

spondente autovettore� calcolando l�autovalore di modulo massimo di A���

Page 53: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� � AUTOVALORI DI MATRICI �M�F�� �

Per evitare il calcolo della matrice inversa presente in

xk�� � A��xk

basta risolvere �per esempio mediante la decomposizione LU il sistema equi�valente

Axk�� � xk�

Il valore

�k �xTk xk��

xTk xkconverger�a al valore

� � limk

�k � �max�A�� �

�min�A

mentre xk�� converger�a al corrispondente autovettore�Ricordando inoltre che

��A� I � ��A�

mentre i corrispondenti autovettori coincidono� il metodo delle potenze in�verse pu�o essere utilizzato per stimare l�autovalore pi�u prossimo ad un valorepre�ssato ed il suo corrispondente autovettore� Posto

B � �A� I��

il metodo delle potenze applicato alla matrice B �e equivalente a

�A� Ixk�� � xk

per cui il valore

�k �xTk xk��

xTk xkadesso converger�a al valore

� � limk

�k � �max�A� I�� ��

���A�

dove con ���A si �e indicato l�autovalore di A pi�u prossimo ad � ed allora

���A ��

�� �

Per poter applicare il metodo delle potenze inverse con shift �e necessarioavere una idea di dove sono dislocati� nel piano complesso� gli autovalori diA� A tale �ne richiamiamo� senza dimostrazione� alcuni classici risultati�

Page 54: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Teorema ���� Se una matrice A �e simmetrica allora i suoi autovalori sonoreali�

Teorema ����� Se una matrice A �e simmetrica e de�nita positiva allora isuoi autovalori sono reali e positivi�

Teorema ����� Per ogni autovalore di A vale la limitazione

j��Aj � kAk

per ogni norma di A�

Teorema ����� di Gerschgorin Data la matrice A di elementi aij� de��niamo

rk �nX

j��j ��k

jakjj � ck �nX

j��j ��k

jajkj

Rk � fzj jz � akkj � rkgCk � fzj jz � akkj � ckg

gli autovalori di A appartengono ai cerchi Rk e Ck� Pi�u precisamente appar�tengono alla regione del piano complesso che �e data dall�intersezione delleunioni dei due insiemi di cerchi��

Sebbene grossolane� le limitazioni fornite dai precedenti teoremi� in alcunicasi risultano su�cienti per stimare alcuni autovalori ed autovettori dellamatrice A�

�� �� Metodi globali

Sono quei metodi che permettono di calcolare tutti gli autovalori della matriceA trasformando� mediante trasformazioni per similitudine� la matrice in unamatrice in forma diagonale o triangolare �forme per le quali �e immediatoil calcolo degli autovalori� E� importante osservare che tali trasformazioniporteranno� in generale per n � � ad una matrice triangolare o diagonalesolo dopo un �numero in�nito� di passi� �E� una naturale conseguenza delteorema di Abel sul calcolo delle radici di equazioni algebriche�

Page 55: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� � AUTOVALORI DI MATRICI �M�F��

Le matrici di trasformazione� che stanno alla base di questi metodi� so�no le stesse che permettono la decomposizione QR di una matrice o la suatrasformazione in forma di Hessemberg�

Fra le diverse tecniche di trasformazione presenteremo in dettaglio quel�la basata sulle matrici di ri�essione di Householder che� a di�erenza dellematrici di rotazione� ri�ettono i vettori a cui sono applicate rispetto degliiper�piani�

Matrici di Householder

De�nizione ����� Dato un vettore u de�niamo matrice di Householder lamatrice

Q � I � �uuT

kuk���

Teorema ����� La matrice Q �e simmetrica� ortogonale ed e�ettua la ri��essione di un vettore rispetto all�iper�piano ortogonale al vettore u�

Dimostrazione ����� La simmetria segue dalla de�nizione di Q� Perl�ortogonalit�a

QTQ �

�I � �uuT

kuk��

�T �I � �uuT

kuk��

��

�I � �uuT

kuk��

��

� I � �uuT

kuk���

�u�uTu

�uT

kuk��� I�

essendo�uTu

�� kuk�� �

Sia z un vettore dell�iper�piano ortogonale ad u� ovvero zTu � �� allora

Qz �

�I � �uuT

kuk��

�z � z

mentre

Qu �

�I � �uuT

kuk��

�u � �u

essendo x � u � �z� �x � Rn� il vettore x �e decomposto nelle sue duecomponenti rispetto a u e z� allora

Qx � Q�u � �z � �u � �z

risulta il ri�esso di x rispetto all�iper�piano ortogonale ad u��

Page 56: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Nasce spontanea la domanda� �e possibile costruire delle matrici Q cheapplicate a sinistra ad una matrice A �eventualmente anche rettangolaregenerino una matrice triangolare superiore � Questa operazione �e simile alladecomposizione LU ma adesso la matrice sinistra �e ortogonale� La risposta�e a�ermativa basta osservare che �e possibile �mediante una opportuna Qrendere nulle tutte le componenti di un dato vettore da un certo indice inpoi� mantenendone inalterata la norma �� In de�nitiva vale il seguente

Teorema ������ Dato un vettore v arbitrario� per ottenere il vettore w �Qv della forma

w �

�wj � vj� j � �� � � � � k � �

wkj kwk � kvk �wj � �� j � k � �� � � � � n

�����

basta prendere u della forma

u �

�������uj � � j � �� � � � � k � �

uk � vk s

nPj�k

v�j �

uj � vj� j � k � �� � � � � n� �����

Dimostrazione ����� Si ha

kuk�� � uTu � v�k �nX

j�k

v�j �vk

vuut nXj�k

v�j �nX

j�k��

v�j �

� �nX

j�k

v�j �vk

vuut nXj�k

v�j �

e

uTv � v�k vk

vuut nXj�k

v�j �nX

j�k��

v�j �

�nX

j�k

v�j vk

vuut nXj�k

v�j �kuk��

��

Page 57: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� � AUTOVALORI DI MATRICI �M�F�� �

essendo w � Qv� si ha

w � Qv �

�I � �uuT

kuk��

�v � v � �u

�uTv

�kuk��

� v � u�

che �e della forma ������

La fattorizzare QR di una matrice A si pu�o quindi ottenere in n�� passicostruendo le opportune matrci Qj di Householder che� al j�esimo passo� an�nullano gli elementi della j�esima colonna della matrice Aj sotto la diagonaleprincipale� dove con Aj � Qj��Qj�� � � �Q�A� abbiamo indicato la matrice alj�esimo passo e A� � A�

Qn��Qn�� � � �Q�Q�A � R�

Osservazione ����� Il prodotto Qy pu�o essere eseguito con n �ops al postodi n� osservando che

Qy �

�I � �uuT

kuk��

�y � y � �uTy

kuk��u � y � �u

essendo

� ��uTy

kuk���

Osservazione ����� Nella ��� il segno viene scelto in modo da ridurrela cancellazione numerica � se vk � �� � se vk � ��

Algoritmo QR

Data la fattorizzazioneA � A� � Q�R� �����

si consideri la matriceA� � R�Q� ����

�e immediato dimostrare che la matrice A� �e simile alla A� infatti

A� � R�Q� ��QT

�Q�

�R�Q� � QT

� �Q�R�Q� � QT�A�Q���

Iterando la decomposizione ����� ed il prodotto ���� si genera unasuccesione di matrici Ak per la quali vale il seguente

Page 58: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

Teorema ������ Se la matrice A ha autovalori reali distinti� j��j � j��j �� � � � j�nj� allora

limk

Ak � R

dove R �e triangolare superiore� se A �e simmetrica R � D diagonale��

Poich�e l�algoritmo QR salva la forma di Hessemberg di una matrice� �e con�veniente trasformare preventivamente la matrice A in forma di Hessemberg�tridiagonale se A �e simmetrica� Tale trasformazione pu�o essere e�ettuata�mantenendo la similitudine� ancora con matrici di Householder in n�� passi�Questa volta si dovr�a moltiplicare� sia a sinistra che a destra� per le matriciQj costruite in modo di azzerare gli elementi della colonna j�esima della ma�trice Aj sotto la sotto�diagonale principale� Senza entrare nei dettagli �chevengono lasciati al lettore come esercizio di veri�ca

H � Qn��Qn�� � � �Q�Q�AQ�Q� � � �Qn��Qn��� �����

Osservazione ����� La ���� �e ovviamente una trasformazione per simili�tudine essendo le Qj ortogonali e simmetriche�

La velocit�a di convergenza dell�algoritmo QR dipende dallo spettro del�la matrice A per cui� nella pratica� si preferisce utilizzare al posto delle��������� le

A� I � A� � I � Q�R�

eA� � R�Q� � I

dove �e un opportuno valore di shift� scelto per accelerare la velocit�a diconvergenza del metodo�

�� �� Analisi dellerrore

Analogamente a quanto fatto per la risoluzione di un sistema lineare �e natu�rale chiedersi di quanto pu�o variare la stima di un autovalore a fronte di unavariazione degli elementi della matrice� Vale il seguente

Teorema ������ di Bauer�Fike Se � �e un autovalore della matrice �A � Ee

X��AX � D � diag���� ��� � � � � �n

Page 59: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� RIEPILOGO �M�F��

alloramin j�� �j � Kp�X kEkp

dove k�kp denota una norma p di matrice��

E� interessante osservare che �e il condizionamento della matrice X chediagonalizza la A che in�uenza il calcolo degli autovalori e non il condiziona�mento di A�

Se A �e simmetrica� essendo diagonalizzabile mediante una matrice orto�gonale� non crea problemi per il calcolo degli autovalori �si pensi alla matricedi Hilbert� cos�i per�da per la risoluzione dei sistemi lineari� ma di cui siriescono a calcolare bene gli autovalori�

Essendo le colonne della matrice X gli autovettori della matrice A� ilteorema di Bauer�Fike evidenzia come siano gli autovettori ad in�uire sulcalcolo degli autovalori �pi�u gli angoli formati da questi sono prossimi a �

�e

pi�u le cose vanno bene�

��� Riepilogo �M�F��

Sintetizzando quanto esposto relativamente al calcolo di autovalori ed auto�vettori abbiamo

�� Metodi locali� permettono il calcolo di un autovalore� Metodo dellepotenze e sue varianti

�a autovalori distinti non creano problemi �autovettori linearmenteindipendenti�

�b velocit�a di convergenza dipende da ����� matrici simmetriche

�����

��

�c shift per accellerare convergenza e calcolare autovettore associatoad un dato autovalore�

�� Metodi globali� tutti gli autovalori� Trasformazioni per similitudine�Algoritmo QR�

�a velocit�a di convergenza legato allo spettro della matrice�

�b shift per accellerare convergenza�

�� Condizionamento� dipende dagli autovettori� non ci sono problemi perle matrici simmetriche essendo gli autovettori ortogonali�

Page 60: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO �� ALGEBRA LINEARE NUMERICA �M�FRONTINI�

����� Esercizi

Esercizio ����� Dato il vettore x � ��� �� �� �� � � costruire la matrice Q diHouseholder in modo che y � Qx � ��� �� y�� �� �� � �

Esercizio ����� Utilizzare il metodo delle potenze per calcolare gli autovaloridi modulo massimo e minimo della matrice di Hilbert di ordine �� e quindi ilsuo condizionamento in norma ��

Esercizio ����� Sapendo che e� � ��� �� �� �� � T �e autovettore della matricedell�esercizio ����� calcolare il corrispondente autovalore utilizzando� nel quo�ziente di Rayleigh� e� ed e� ��e�� dove �e� �e una perturbazione di e�� Provareper varie perturbazioni �e� e commentare i risultati ottenuti�

Page 61: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo �

Zeri di funzioni�M�Frontini�

La risoluzione di molti problemi conduce al calcolo degli zeri di una funzione�si veda per esempio il problema della determinazione del tasso di crescita nelmodello di Maltius presentato nell�introduzione� oppure il problema dellaricerca del minimo �massimo di una funzione e tutti quei problemi legatialla massimizzazione del pro�tto o minimizzazione della spesa�

Considereremo in dettaglio il caso monodimensionale estendendo� ove pos�sibile� alcuni risultati al caso di sistemi non lineari� Partiremo quindi dalproblema� data f�x � C��a� b con f�af�b � �� trovare � � �a� b tale chef�� � ��

Un primo algoritmo per la determinazione di � ci viene fornito dal teoremadegli zeri visto nel corso di Analisi�

Teorema �� �� degli zeri Sia f�x � C��a� b con f�af�b � �� alloraesiste almeno uno � � �a� b tale che f�� � ���

Se alle ipotesi del teorema degli zeri si aggiunge la monotonia della f�xsi ha anche l�unicit�a del punto ��

�� Metodo di bisezione �M�F��

Sotto le ipotesi precedenti de�niamo la seguente sequenza di operazioni�

x� � a� x� � b� i � �

��

Page 62: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

xi�� �xi � xi��

�� valuta f�xi��

f�xi�� � � � � � xi�� FINE�

altrimenti f�xif�xi�� � � � xi�� � xi��� i � i � �

altrimenti f�xi��f�xi�� � � � xi � xi��� i � i � �

continua�

Nella pratica non si arriver�a mai ad avere f�xi�� � �� per cui nellasequenza precedente bisogna introdurre un opportuno test di arresto in modod�ottenere un algoritmo utilizzabile� Un buon test �e il seguente

jxi � xi��j � toll ����

dove toll �e la precisione richiesta sul risultato� Essendo � � �xi� xi�� la ����fornisce l�intervallo di indeterminazione di ��

Osservazione ����� Nel calcolo del punto medio �e pi�u stabile l�uso dellaseguente formula

xi�� � xi �xi�� � xi

��

che permette di evitare l�over�ow ed assicura� come deve essere� che xi�� �xi�� � xi xi � xi�� � xi���

E� possibile sapere a priori il numero di iterazioni necessarie per raggiun�gere la tolleranza richiesta essendo

xi � xi�� �b� a

�i��

per cui� dalla �����

i � � �log� b�a

toll

log ��

Il costo computazionale del metodo �e dato dal calcolo� ad ogni passo� dellaf�xi��� La convergenza �e lenta� si dimezza l�intervallo ad ogni passo� ma siha il vantaggio d�avere sempre una stima per difetto e una per eccesso dellaradice� Il metodo non pu�o essere usato per problemi in dimensione maggioredi ��

Page 63: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODO DI BISEZIONE �M�F�� ��

����� Falsa posizione

E� una variante del metodo di bisezione� l�unica di�erenza �e che il punto xi��

viene calcolato come intersezione della retta congiungente i punti �xi� f�xie �xi��� f�xi�� ovvero

xi�� �xif�xi��� xi��f�xi

f�xi��� f�xi

e� per il resto� si procede analogamente� Un inconveniente di questa variante�e che� in generale� non fornisce pi�u l�intervallo di indeterminazione della � e�se la f�x �e monotona� uno dei due estremi �a o b non viene mai modi�cato�confronta �gura ����

−1 −0.5 0 0.5 1 1.5 2 2.5 3−4

−3

−2

−1

0

1

2

3

4

5

6

f(x)

^

>

Corde

x

f(x)

o

Falsa Posizione

Figura ����

Anche il test di convergenza deve essere modi�cato per contemplare sia ilcaso in cui la derivata prima di f�x sia �piccola� �funzioni piatte vicono a �che quando �e �grande� �funzioni ripide vicono a �� Un buon test di arresto�e dato quindi dalla contemporanea veri�ca delle seguenti disuguaglianze

jf�xi��j � zero

jxi � xi��j � toll

dove zero e toll sono in generale diversi �confronta �gure ��� e ����

Page 64: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

0 ξa

b

Figura ����

�� Metodi di punto �sso �M�F��

Sono i metodi che si ottengono trasformando il problema del calcolo dellozero di f�x nel problema equivalente della costruzione di una successione�convergente

xk�� � g�xk ����

dove g�x �e scelta in modo che

� � g�� ����

sef�� � ��

Osservazione ����� Si invita lo studente a so�ermarsi sulle analogie con imetodi iterativi stazionari� presentati per la risoluzione dei sistemi lineari�

Vediamo un semplice esempio

Esempio �����

f�x � x� � �x � �

g��x �x� � �

g��x �p

�x� �� x � �

�g��x � x� � x � �

Page 65: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI DI PUNTO FISSO �M�F�� �

ξ0

a

b

Figura ����

Le tre gi�x ammettono tutte � � � come �punto unito�� ma non tuttee tre vanno bene per generare una successione convergente partendo� peresempio� da x� � �

��

La convergenza della ���� a � dipender�a quindi dalla g�x�Valgono iseguenti teoremi�

Teorema ����� Se la funzione g�x � C��a� b e se g�a� b � �a� b � alloraesiste almeno un punto unito � � �a� b per la trasformazione g�x�

Dimostrazione ����� Essendo

g�a� b � �a� b

si ha

a � g�a � b

a � g�b � b�

Se g�a � a oppure g�b � b si ha la tesi con � � a oppure � � b�altrimenti si ha

g�a� a � �

g�b� b � ��

Page 66: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

PostoF �x �� g�x� x

allora F �x � C��a� b ed inoltre F �aF �b � �� per cui� applicando il teoremadegli zeri F �� � � � g��� � � � � g����

Teorema ����� Se la funzione g�x � C��a� b e se g�a� b � �a� b � ed inol�tre jg��xj � L � � allora esiste un solo punto unito � � �a� b per latrasformazione g�x�

Dimostrazione ����� per assurdo Supponiamo che esistano due distintipunti uniti �� e ��

�� � g���� �� � g���� �� �� ��

per il teorema della media �� � ���� �� per cui

j�� � ��j � jg���� g���j � jg����� � ��j� L j�� � ��j � j�� � ��j

che �e assurdo��

Teorema ����� Se la funzione g�x � C��a� b e se g�a� b � �a� b � conjg��xj � L � � allora la sucessione xk�� � g�xk converge a � ed inoltre�de�nito en � xn � �� si ha

jenj � Ln

�� Ljx� � x�j ��

Osservazione ����� Sebbene le condizioni indicate nel teorema ����� pos�sano essere di�cili da veri�care a priori� il teorema fornisce una chiara ri�sposta riguardo alle condizioni di convergenza e� se �e nota L� il numero diiterazioni necessarie per ottenere � con una pre�ssata tolleranza�

Osservazione ����� Vale la pena di osservare che l�equazione ��� �e l�e�quazione risolvente il sistema �

y � xy � g�x

per cui il metodo ��� pu�o essere interpretato geometricamente come indicatoin �gura ����

Page 67: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� METODI DI PUNTO FISSO �M�F�� ��

y=x

0x

y

g(x)

Figura ����

Per quanto riguarda la velocit�a di convergenza e l�ordine di un metododel tipo ���� diamo la seguente de�nizione

De�nizione ����� Un metodo del tipo ��� si dice di ordine k se

limn

en��

ekn� c costante�

dove en � xn � �� g�i �� � � per i � k e g�k �� �� ��

I metodi per cui g��� �� �� vengono detti del primo ordine �a convergenzalineare� infatti

en�� � xn�� � � � g�xn� g�� � � � �Taylor � � � ����

� g���en �g�� ��

��e�n � � � �� g�k �n

k�ekn

dove n � �xn� � � Per k � � si ha

en�� � g��nen

e per il teorema ������ essendo

limn

xn � �

sar�a anchelimn

n � �

Page 68: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

per cui

limn

en��

en� g��� � L � ���

In generale se g��� � g�� �� � � � � � g�k�� �� � � dalla ���� segue

limn

en��

ekn�

g�k ��

k�� c�

Osservazione ���� Mentre per i metodi del primo ordine il guadagno adogni passo �e costante tanto pi�u si guadagna quanto pi�u L �e minore di ��nei metodi di ordine superiore pi�u si �e vicini alla soluzione maggiore �e ilguadagno al passo sucessivo poco importa il valore di c� ci�o che conta �e chek � ��

� Metodo di Newton �M�F��

Si pu�o ottenere un metodo del secondo ordine dalla seguente osservazione�sia

g�x � x � h�xf�x

dove la nuova funzione h�x verr�a scelta in modo che g��� � �� E� immediatoosservare che� se f�� � � allora � � g��� costruiamoci la h�x� derivando siha

g��x � � � h��xf�x � h�xf ��x

per x � �� essendo f�� � �� si ha

g��� � � � h��� f����

�h��f ��� � � � h��f ���

da cui

h�� � � �

f ����

Otteniamo quindi il metodo� almeno del secondo ordine �metodo di New�ton� nella forma

xn�� � xn � f�xn

f ��xn�

Vale il seguente

Page 69: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� METODO DI NEWTON �M�F�� �

Teorema ����� Se la f �� �x �e continua e la f ��x �� � in un aperto con�tenente �� allora esiste un � � � per cui � jx� � �j � � il metodo di Newtonconverge quadraticamente a ���

Osservazione ����� Il precedente teorema a�erma che� pur di partire �ab�bastanza vicino� alla soluzione� la convergenza �e quadratica� Nulla dice se siparte �lontano��

Per garantire la convergenza �comunque si parta� bisogna richiedere chela f�x non cambi concavit�a in �a� b � ovvero

Teorema ����� Se f�x � C��a� b � f�af�b � �� f �� �x �e di segno costantein �a� b e f ��x �� �� allora il metodo di Newton converge quadraticamente a� se si parte da un x� � �a� b tale per cui f�x�f

�� �x� � ���

Una chiara interpretazione del precedente teorema �e data nella �gura �� �

0 x

y

ξ

f(x)x

x1

x2

Figura �� �

Il metodo di Newton �e� in generale� un metodo del secondo ordine ma adogni iterazione richiede la valutazione di due funzioni la f�x e la f ��x� Sela f ��x non �e nota� o se �e molto oneroso il suo computo� il metodo non pu�oessere usato o pu�o risultare lento� Vedremo pi�u avanti come si pu�o ovviare aquesti inconvenienti�

Se la radice � �e multipla �per esempio doppia il metodo non �e pi�u delsecondo ordine ma decade al primo �con costante L � �

�� infatti

g��x � �� f ��x� � f�xf �� �x

f ��x��

f�xf �� �x

f ��x�

Page 70: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

per cui� utilizzando il teorema di de L�Hopital�

limx�

g��x � limx�

f ��xf �� �x � f�xf �� �x

�f ��xf �� �x�

� limx�

f �� �x� � �f ��xf �� �x � f�xf �� �x

��f �� �x� � f ��xf �� �x�

���

E� facile dimostrare che se p �e la molteplicita� della radice �� allora

xn�� � xn � pf�xn

f ��xn

�e ancora del secondo ordine� mentre il metodo classico di Newton sarebbe delprimo con costante L � �� �

p�

Presentiamo due varianti al metodo di Newton che permettono d�evitareil calcolo della f ��x�

����� Metodo delle secanti

Si approssima la derivata prima con il rapporto incrementale per cui

f ��xn � f�xn� f�xn��xn � xn��

ed il metodo diviene

xn�� � xn � f�xn�xn � xn��f�xn� f�xn��

Il metodo delle secanti richiede ad ogni passo il calcolo della funzione inun solo punto� La convergenza �e iperlineare� ovvero

limn

en��

epn� c� p �

p � �

�� �����

La successione generata dal metodo delle secanti NON coincide con quellagenerata dal metodo della falsa posizione �cfr� �gura ����

Page 71: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� SISTEMI NON LINEARI �M�F�� ��

0 x

0

y f(x)x

ξ1

23

4c

4f

4c:=corde

4f:=falsa posizione

Figura ����

����� Metodo di Ste�ensen

Si approssima la f ��x con

f ��xn � f�xn � f�xn� f�xn

f�xn

per cui il metodo diviene

xn�� � xn � f�xn�

f�xn � f�xn� f�xn�

Il metodo di Ste�ensen richiede il calcolo della funzione in due punti adogni passo �come Newton ed �e un metodo del secondo ordine�

� Sistemi non lineari �M�F��

����� Metodi di punto �sso

Dato il sistema non linearef�x � � ���

lo si trasforma nel problema equivalente

xk�� � ��xk ����

Page 72: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

per cui se� �jf�� � � allora � � ����Indicata con

J���xk �

�������

����x�

����x�

� � � ����xn��

����xn

����x�

����x�

� � � ����xn��

����xn

������

� � ����

�����n���x�

��n���x�

� � � ��n���xn��

��n���xn

��n�x�

��n�x�

� � � ��n�xn��

��n�xn

��������x�xk

����

la matrice jacobiana di ��x valutata nel punto xk� vale il seguente

Teorema ���� Condizione su�ciente per la convergenza del metodo ����e che ��J���xk

�� � L � �� �k��

Osservazione ���� Si confronti il teorema ����� con il teorema ������

����� Metodo di Newton per sistemi

E� la generalizzazione ai sistemi non lineari del metodo di Newton visto alparagrafo ����� Per calcolare la soluzione del sistema non lineare ��� sicostruisce la successione

xk�� � xk � hk

dove il vettore hk �e ottenuto come soluzione del sistema lineare

J�f�xkhk � �f �xk ����

dove con J�f�xk si �e indicata �cfr� ���� la matrice jacobiana di f�xvalutata nel punto xk� Il metodo pu�o essere dedotto facilmente sviluppandola f�xk in un intorno di xk� Dopo aver posto � � xk � hk� si ha

f�� � � � f�xk � hk � ���sviluppando������������f��xk � hk�

��f��x�

�x�xk

� � � �� hkn

��f��xn

�x�xk

� � � � � �

� � � � � � �fn�xk � hk�

��fn�x�

�x�xk

� � � �� hkn��fn�xn

�x�xk

� � � � � �

���

Page 73: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� SISTEMI NON LINEARI �M�F�� ��

dove con hki si �e indicata la i�esima componente di hk� La ���� trascurando itermini di ordine superiore al primo� pu�o essere scritta nella forma matriciale�����

Vale la pena d�osservare che nella determinazione del vettore xk�� la partepi�u onerosa computazionalmente �e la risoluzione del sistema lineare ����� percui pu�o essere utile tenere �ssa la matrice J�f�xk per un numero �ssato�diciamo r di iterazioni� in modo da poter utilizzare la decomposizione LUgi�a e�ettuata al passo k� e riaggiornare la matrice solo al passo �k � r�

Osservazione ���� Anche nel caso monodimensionale esiste una varianteanaloga che consiste nel non calcolare ad ogni passo la f ��xk mandando� adogni passo k� la retta parallela alla precedente non pi�u la tangente passanteper il punto �xk� f�xk cfr� �gura �� �

0 x

y f(x)x

NON tangenti

x1

xξ x23

Rette parallele

Figura ����

Molto costoso �e anche il preventivo calcolo simbolico della matrice jaco�biana �n� derivate parziali da calcolare � per cui� quando queste derivatesono molto complicate o non sono computabili �non si hanno le f�x in for�ma analitica� si pu�o operare con una variante� simile al metodo delle corde�Precisamente si approssimano gli elementi della matrice jacobiana con deirapporti incrementali nella forma��

�fj�xi

�x�xk

� Bji �fj�xk � ekhki� fj�xk

hki� ekj � �kj

dove� per esempio� hki � xki � xk��i�

Page 74: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

�� Zeri di polinomi �M�F��

Una trattazione a parte pu�o essere fatta per le equazioni algebriche in quan�to� oltre ai metodi visti in precedenza� si possono considerare altri metodinumerici sfruttando il legame che esiste fra un polinomio e la sua matrice diFrobenious �o Companion matrix � vale infatti il seguente

Teorema ����� Dato il polinomio monico

pn�x � xn � a�xn�� � � � �� an��x � an �����

questi risulta essere il polinomio caratteristico della matrice

Cn �

�������

� � � � � � �an�

� � �� � �

��� �an��� �

� � � ����

���� � �

� � � � �a�� � � � � � �a�

�������� �� �����

In virt�u del risultato precedente �e quindi possibile calcolare gli zeri di����� mediante il calcolo degli autovalori di ����� �cfr� capitolo sugli auto�valori�

����� Schema di Horner

Il polinomio ����� pu�o essere scritto nella forma �dovuta ad Horner

pn�x � an � x�an�� � x�an�� � x�� � � � x�a� � x � � � �����

�si lascia al lettore la facile veri�ca�La forma ����� permette di valutare il polinomio pn�x nel punto

mediante il seguente schema��� � ��k � �k�� � ak

� k � �� �� � � � � n �����

per cui

pn� � �n�

Page 75: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� ZERI DI POLINOMI �M�F�� �

E� facile veri�care che i valori �k �k � �� � � � � n� � sono i coe�cienti delpolinomio quoziente pn���x della divisione per �x� di pn�x� mentre �n�e il resto r�� Lo schema ����� rappresenta quindi un modo compatto pere�ettuare la divisione sintetica

pn�x � �x� pn���x � r�� �����

Dimostrazione ����� Scriviamo il polinomio

pn���x � xn�� � ��xn�� � � � �� �n��x � �n��

dalla ����� eguagliando i coe�cienti dei termini di ugual grado� tenendoconto della ���� si ha

��� ����� � � a� �� � � a��� � �� � a� �� � �� � a�� � � � � ��n�� � �n�� � an�� �n�� � �n�� � an��r� � �n�� � an r� � �n�� � an � �n

��

Analogamente si pu�o calcolare la derivata prima di pn�x in x � operando la stessa divisione con pn���x� infatti� posto

pn���x � �x� pn���x � r�

si hapn�x � �x� �pn���x � �x� r� � r�

e derivando

p�n�x � ��x� pn���x � �x� �p�n���x � r�

da cui�p�n� � r���

Si lascia al lettore� per esercizio� la veri�ca che

p�m n � � m�rm� m � �� �� � � � � n�

Lo schema di Horner� o divisione sintetica� fornisce un algoritmo e�cienteper valutare il polinomio e la sua derivata prima in un punto permettendo�con un costo di �n �ops� di eseguire un passo del metodo di Newton

xn�� � xn � pn�xn

p�n�xn�

Page 76: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

����� Radici multiple

Abbiamo gi�a visto che� nel caso di radici multiple� il metodo di Newton nonconverge pi�u quadraticamente� Per quanto riguarda i polinomi �e semprepossibile operare con polinomi che hanno radici semplici� infatti il polinomio

pn�r�x �pn�x

qr�x

ha tutte e sole le radici di pn�x ma semplici� essendo

qr�x � MCD fpn�x� p�n�xg

dove con MCDfg si �e indicato il massimo comun divisore fra i polinomi inparentesi fg�

Ovviamente se le radici di pn�x sono semplici allora qr�x � q��x � k�costante� Ricordiamo che il MCD pu�o essere calcolato agevolmente con ilnoto algoritmo euclideo�

� Riepilogo �M�F��

Possiamo sintetizzare i risultati presentati nel modo seguente

�� f�x � C��a� b �funzione poco regolare

�a metodo di bisezione�

�b metodo della falsa posizione�

�� f�x � Ck�a� b �k � �� funzione regolare

�a metodi di ordine k �punto �sso�

�� Complessit�a computazionale

�a � valutazione per passo� bisezione� falsa posizione e corde�

�b � valutazioni per passo� Newton e Ste�ensen�

�� Velocit�a di convergenza

Page 77: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� RIEPILOGO �M�F�� ��

�a metodi del primo ordine �lineare� dipende da L �jg��xj � L � ��

�b ���� corde�

�c � Newton e Ste�ensen�

� Equazioni algebriche�

�a tutti i metodi precedenti�

�b metodi dell�algebra lineare �calcolo degli autovalori della compa�nion matrix�

I metodi di Newton� corde ed i metodi di punto �sso possono essere estesial caso di sistemi non lineari�

����� Esercizi

Esercizio ���� Quanti passi sono necessari per stimare� con errore minoredi ����� la radice di

ex � � � �

utilizzando il metodo di bisezione partendo dall�intervallo ��� � �

Esercizio ���� Utilizzare� per calcolare la radice dell�equazione

�e�x � ��ex � � � �

�� il metodo di Newton

�� il metodo

xn�� �xne

xn � exn � �

exn�

partendo sempre da x� � �� Quale metodo converge prima� Perch�e�

Esercizio ���� Utilizzare� se possibile� i seguenti metodi

�� xn�� � � lnx�

�� xn�� � e��x

�� xn�� � xn�e��xn

per calcolare la radice di�x � lnx � ��

Quale metodo converge pi�u velocemente� Perch�e� Si fornisca un metodomigliore di quelli proposti�

Page 78: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � ZERI DI FUNZIONI �M�FRONTINI�

Esercizio ��� Dato il sistema non lineare�x � y�

��

y � x�

��

studiare la convergenza del seguente metodo iterativo�xn�� � �y�n

��

yn�� � �x�n�

� �

Esercizio ���� Calcolare le radici del seguente polinomio

x� � �x� � x� �� � �

�� utilizzando il metodo di Newton�

�� il calcolo degli autovalori della companion matrix associata�

Page 79: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo

Teoria dell�approssimazione�M�Frontini�

Il problema che vogliamo a�rontare ora �e quello della ricostruzione di unafunzione� o semplicemente la determinazione del suo valore in un punto� notii valori che la funzione stessa assume in un numero �nito di punti assegnati�E� evidente che cos�i formulato il problema non �e ben posto� nel senso che pu�oavere un numero in�nito di soluzioni o nessuna soluzione� Si pensi al caso incui la funzione deve essere continua e si hanno due punti con uguale ascissa ediversa ordinata� Per ripristinare l�esistenza ed unicit�a bisogner�a aggiungeredelle condizioni

�� sulla regolarit�a della funzione

�� sul tipo di approssimazione da usare

�� sulle funzioni approssimanti�

In generale� data una funzione f�x nota in n punti yi � f�xi �i ��� ����� n� si cerca un�approssimante della forma

g�x �nXi��

i�i�x

dove i pesi i sono scelti secondo un dato criterio di merito� mentre le�i�x sono funzioni opportune che veri�cano date condizioni di regolarit�ae sono facilmente manipolabili �ovvero facilmente computabili� derivabili edintegrabili�

Page 80: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

Come esempi di funzioni �i�x si possono considerare

�� i polinomi di grado n���

�� i polinomi a tratti di grado m�

�� le funzioni spline di grado m�

�� i polinomi trigonometrici di Fourier�

� le funzioni razionali�

Come criteri di merito per determinare gli i si possono considerare

�� l�interpolazione �g�xi � f�xi � yi� i � �� ����� n�

�� i minimi quadrati �norma L� �minai

�Pni�� �g�xi� yi

���

�� minimo massimo errore �minimax� norma L �minai

nmax

ijg�xi� yij

o�

�� minima somma dei moduli �norma L� �minaifPn

i�� jg�xi� yijg�

che danno una misura della distanza fra la funzione incognita f�x el�approssimante g�x� A seconda dei campi d�applicazione �e meglio utilizzareun criterio piuttosto che un altro

L�interpolazione polinomiale �e utile quando si hanno pochi dati �n piccoloper cui si ottiene un polinomio di grado basso� La regolarit�a del polinomiointerpolante �e utile poi nella deduzione di formule approssimate per il calcolointegrale e per l�intergrazione di equazione di�erenziali �cfr� capitoli e ��

L�interpolazione mediante spline �e utile nelle applicazioni di gra�ca �CAD�non a caso il termine �spline� indica in inglese il �curvilineo mobile� deldisegno tecnico�

I minimi quadrati sono utilizzati nei problemi di identi�cazione di para�metri e modelli� quando si hanno a disposizione pi�u dati ma a�etti da errori�di misura casuali ecc��

Il criterio del minimo massimo errore �e essenziale quando si deve garantireuna certa accuratezza �routine di calcolo delle funzioni elementari�

Contrariamente a quanto avviene in Rn� nel caso dell�approssimazionedi funzioni le norme non sono pi�u equivalenti come mostrato dal seguenteesempio

Page 81: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��

Esempio � �� Data la funzione f�x identicamente nulla in ��� � conside�riamo la successione di funzioni gk�x cfr� �gura ��� cos�i de�nite

gk�x �

�k�k�x� � �

k�� x � �

k�

�k�k�x� � �k�� x � �

k�

� altrove�

0 x

g(x)k

1/k

k

2 2/k 2 3/k 2 3

k

g(x)

Figura ����

Considerando le tre classiche de�nizioni di norma di funzioni� si ha�

kf � gkk� � �k

kf � gkk� �q

��

kf � gkk � k

per cui la succesione gk�x converge a f�x in norma L� ma non converge innorma L� e L��

Valgono i seguenti teoremi

Teorema � �� La convergenza in norma L garantisce la convergenza del�le altre due norme�

Dimostrazione ���� Dalla de�nizione di norma e dal teorema della me�dia si ha

kf � gk� �

Z b

a

jf�x� g�xj dx � ���teorema media���

Page 82: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

� jf��� g��jZ b

a

dx �� max

xjf�x� g�xj � �b� a �

� C kf�x� g�xk �

essendo kf�x� g�xk � maxxjf�x� g�xj e C � b� a� Analogamente per

la norma L���

De�nito con !N lo spazio dei polinomi di grado N valgono i seguentiteoremi

Teorema � �� di Weierstrass Data una funzione f�x � C��a� b allora��� � � esiste un N tale per cui �e possibile trovare un polinomio p�x � !N�

per cui

kf � pk � ���

Teorema � �� di miglior approssimazione Data una funzione f�x �C��a� b e un numero positivo N � esiste un unico polinomio p��x � !N percui kf � p�k � kf � pk � �p � !N ��

Osservazione � �� Questo teorema vale per ogni norma� anche se p��xvaria al variare della norma�

Osservazione � �� Il teorema di Weierstrass a�erma l�esistenza del poli�nomio p�x di grado N ma non dice come il grado N dipenda da ��

�� Interpolazione �M�F��

Il problema che vogliamo a�rontare ora �e la costruzione di una funzione g�xinterpolante la f�x in n � � nodi xi� La g�x sar�a quindi caratterizzata dalfatto che

g�xi � f�xi � yi� i � �� �� �� ���� n�

Fra tutte le possibili g�x interpolanti considereremo solo il caso deipolinomi� delle funzioni spline e dei polinomi trigonometrici�

Page 83: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� ��

����� Polinomi di Lagrange

Per costruire g�x � !n de�niamo i seguenti n � � polinomi di Lagrange �digrado n lj�x de�niti sui nodi d�interpolazione e caratterizzati dal nodo xj

lj�x ��x� x��x� x� � � � �x� xj���x� xj�� � � � �x� xn���x� xn

�xj � x��xj � x� � � � �xj � xj���xj � xj�� � � � �xj � xn���xj � xn�

De�nito il polinomio monico �con coe�ciente della potenza di gradomassimo uguale ad � di grado n � �

�n���x �nYi��

�x� xi

si ha

lj�x ��n���x

�x� xj��n���xj�

Osservazione ���� Ogni polinomio lj�x gode delle seguenti propriet�a

�� �e un polinimio di grado n�

�� lj�xi � �ij simbolo di Kroneker vale � nel nodo x � xj e vale � neinodi x �� xj�

Si ottiene quindi �esistenza il polinomio interpolante

g�x � pn�x �nXi��

yili�x

per il quale si veri�ca banalmente che

yi � pn�xi�

Se i nodi xi sono distinti vale il seguente

Teorema ���� unicit�a Se i nodi xi� i � �� �� ���� n� sono distinti il polino�mio pn�x � !n interpolante �e unico�

Page 84: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

Dimostrazione ���� per assurdo supponiamo che esista qn�x � !n taleche

qn�x �� pn�x� qn�xi � pn�xi � yi� i � �� �� ���� n

de�niamorn�x � pn�x� qn�x � !n

si avr�arn�xi � pn�xi� qn�xi � yi � yi � �

per cui il polinomio rn�x di grado n avendo n � � zeri �e identicamentenullo� da cui

qn�x � pn�x��

����� Sistema di Vandermonde

Il polinomio interpolante pu�o essere scritto nella forma

pn�x � a� � a�x � a�x� � � � �� an��xn�� � anx

n

gli n� � parametri ai possono essere determinati mediante la risoluzione delsistema lineare ����� che si ottiene imponendo il passaggio del polinomio pergli n � � punti d�interpolazione �xi� yi�����

a��a�x� � a�x��� � � � �anx

n� � y�

a��a�x� � a�x��� � � � �anx

n� � y�

� � � � � � � � � � � � �a��a�xn � a�x

�n� � � � �anx

nn � yn

����

che in foma matriciale diviene

V a � y

dove

V �

���� � x� � � � xn�� x� � � � xn����

������

���� xn � � � xnn

������e la matrice di Vandermonde �che �e noto essere non singolare quando xi �� xjper i �� j�

Page 85: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� �

Esempio ���� Costruiamo i polinomi di Lagrange ed il polinomio interpo�lante i seguenti punti

xi � � �yi �� � �

si ha

l��x ��x� ��x� �

��� ���� �� l��x �

�x� ��x� �

��� ���� �� l��x �

�x� ��x� �

��� ���� ��

e quindi

p��x � ��x� � �x � �

�� �

x� � �x

��� �

x� � x

�� x� � �x� ��

Si lascia al lettore la costruzione dello stesso polinomio interpolante con latecnica di Vandermonde��

����� Stima dellerrore

E� possibile valutare� sotto opportune ipotesi di regolarit�a della funzione f�x�l�errore che si commette sostituendo per ogni x il polinomio interpolante allafunzione f�x stessa� vale infatti il seguente

Teorema ���� Se la funzione f�x � Cn���a� b � de�nito �n���x �nQi��

�x�xi� con xi � �a� b � allora �x � �a� b

en�x � f�x� pn�x �f �n�� ��

�n � ���n���x� ����

Dimostrazione ����� De�niamo la funzione u�t � Cn���a� b come

u�t � f�t� pn�t� f�x� pn�x

�n���x�n���t

la u�t ha n � � zeri t � x� t � xi i � �� �� ���� n� applicando il teorema diRolle generalizzato alla derivata n���esima di u�t si ha

u�n�� �� � f �n�� ��� f�x� pn�x

�n���x�n � �� � � ����

pn�t �e un polinomio di grado n per cui p�n�� n �t � �� �n���t �e un polinomio

di grado n � � monico per cui ��n�� n�� �t � �n � �� mentre � � ��x �e un

punto dell�intervallo determinato da x ed i nodi xi� dalla ��� segue la tesi��

Page 86: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

Osservazione ���� La ��� �e una valutazione �locale� dell�errore dipendedal punto x� Se per x � xi si ottiene come ovviamente deve essere errorenullo� possono esistere valori di x per cui la ��� pu�o fornire valori troppograndi si pensi ad x lontani dall�intervallo de�nito dai nodi xi� Questo nondeve stupire in quanto �e logico che lontano dai punti dati dove si hanno leinformazioni su f sia pi�u di�cile ricostruire la funzione�

Se��f �n�� ��

�� � Kn�� costante� �� � �a� b allora dalla ���� si ottiene

maxxjen�xj � max

xjf�x� pn�xj � Kn��

�n � ��maxxj�n���xj �

e� se i nodi sono equidistanti in �a� b � maxx j�n���xj �e ottenuto per x prossi�mo al bordo di �a� b �in generale l�approssimazione �e pi�u accurata al centrodell�intervallo che ai bordi� perch�e��

Osservazione ���� L�uso dei polinomi di Lagrange ha l�inconveniente chel�aggiunta di un nodo d�interpolazione comporta il ricalcolo di tutti i polinomili�x� Inoltre anche il calcolo della derivata del polinomio interpolante� scrittonella forma

pn�x �nXi��

yili�x

risulta poco agevole�

����� Di�erenze divise

Un generico polinomio di grado n� che assume pre�ssati valori yi in n � �nodi xi �i � �� �� ���� n� pu�o essere scritto nella forma

pn�x � a� � a��x� x� � a��x� x��x� x� � � � � ����

� � �� an�x� x��x� x� � � � �x� xn��

dove i coe�cienti ai vengono determinati imponendo i pre�ssati valori yi�E� immediato osservare che

pn�x� � a�� a� � y��

ed inoltre

pn�x� � a� � a��x� � x�� a� �y� � y�x� � x�

Page 87: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� ��

mentre� dopo alcuni calcoli�

pn�x� � a� � a��x� � x� � a��x� � x��x� � x�� a� �

y��y�x��x� �

y��y�x��x�

x� � x��

Senza entrare nei dettagli si pu�o dimostrare che

ak � f �x�� x�� ���� xk

dove con f �x�� x�� ���� xk si �e indicata la di�erenza divisa k�esima sui nodix�� x�� ���� xk de�nita dalle seguenti relazioni

f �xi � � f�xi

f �xi� xi�� � �f �xi�� � f �xi

xi�� � xi���

����

f �xi� xi��� ���� xi�k � �f �xi��� ���� xi�k � f �xi� ���� xi�k��

xi�k � xi�

Le ��� permettono di costruire� ricorsivamente� tutti i coe�cienti ak� Ilvantaggio di scrivere un polinomio nella forma delle ���� risiede nel fattoche l�aggiunta di un nuovo nodo xn�� non richiede il ricalcolo dei primi ncoe�cienti� Si ha infatti

pn���x � pn�x � an���x� x��x� x� � � � �x� xn���x� xn

e baster�a calcolare an�� mediante le ��� �

����� Interpolazione di Hermite

Se oltre ai valori della funzione in dati punti si conoscono anche i valoridella�e derivate nei punti si pu�o cercare �se esiste il polinomio che interpolasia la funzione che le sue derivate� Un esempio in questo senso �e dato dalpolinomio di Taylor che� come noto� ha gli stessi valori della funzione e dellesue prime n derivate in un dato punto x�

pn�x �nX

j��

f �j �x�

j��x� x�

j

che esiste ed �e unico essendo ricavabile dalle condizioni

p�j n �x � f �j �x�� j � �� �� ���� n

Page 88: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

che portano al sistema lineare triangolare superiore� non singolare���������a� �a�x� �a�x

�� � � � �� � � � �anx

n� � f�x�

�a� ��a�x� � � � �� � � � �nanxn��� � f ��x�

� � � � � � �� � � � � � � � � � ��n� ��an�� �n�anx� � f �n�� �x�

�n�an � f �n �x�

Osservazione ��� Val la pena di osservare che nodi� valori della funzio�ne e derivate non possono essere assegnati arbitrariamente� pena la nonesistenza o la non unicit�a del polinomio interpolante� come dimostra ilseguente

Esempio ����� Costruire se esiste�� la parabola passante per i punti��� e ��� con derivata �� nel punto di ascissa �� Dalle condizioni date siha

p��x � a� � a�x � a�x�

e quindi il sistema lineare �a� � ��a� ��a� � �a� ��a� � ��

che �e impossibile�� Cosa succede se la derivata in � vale ���

Se conosciamo la funzione �yi e la derivata prima �y�i negli n � � nodixi �e possibile costruire il polinomio interpolante di Hermite di grado �n � �dato da

p�n���x �nX

j��

�yjAj�x � y�jBj�x

�dove Aj�x e Bj�x sono polinomi di grado �n � � de�niti da

Aj�x ���� ��x� xjl

j�xj�l�j �x

Bj�x � �x� xjl�j �x

e lj�x sono i polinomi di Lagrange sui nodi xi�Si pu�o dimostrare �si lasciano i conti al lettore che

Aj�xi � �ji� Bj�xi � �A

j�xi � �� B�

j�xi � �ji�i� j

Page 89: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� �

per cuip�n���xi � yi� p��n���xi � y�i� i � �� �� ���� n�

Con tecnica analogo a quella usata per l�interpolazione di Lagrange si pu�odimostrare che per l�errore si ha

e�n���x � f�x� p�n���x �f ��n�� ��

��n � ����n���x�

����� Spline

Quando il numero n di nodi in cui �e nota la funzione f�x �e elevato �problemidi gra�ca ecc�� non �e conveniente utilizzare il polinomio interpolante inquanto

�� �e di grado troppo elevato�

�� pu�o presentare troppe oscillazioni �gra�co non �liscio�� non smooth�

�� la costruzione diviene instabile �matrice di Vandermonde mal condizio�nata�

Per ovviare a questi inconvenienti �e preferibile l�interpolazione polino�miale a tratti o l�uso delle spline� L�idea dell�interpolazione a tratti nascedall�esigenza di mantenere basso il grado del polinomio� per cui gli n nodidi interpolazione vengono divisi in k gruppi di m �n � k � m in modo dacostruire k polinomi di grado m� � che si raccordano� a due a due� nel nodocomune �la funzione risultante �e quindi a priori solo continua� Per le splinesi aggiunge la richiesta di regolarit�a della funzione �in generale C�m�� cos�ida ottenere una curva approssimante �pi�u liscia��

Per la costruzione dell�approssimante data dall�interpolazione polinomialea tratti basta applicare k volte quanto visto nei paragra� precedenti� Perquanto riguarda le spline so�ermeremo la nostra attenzione su quelle linearie le cubiche�

De�nizione ���� Dati i nodi x��x�� ���� xn � �a� b si de�nisce spline inter�polante di grado k una funzione Sk�x che gode delle seguenti propriet�a

�� Sk�x � Ck���a� b � regolarit�a

�� Sk�xj � f�xj� j � �� �� ����� n� interpolazione

Page 90: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

�� Skj�x � !k �x � �xj� xj�� � j � �� �� ����� n� �� forma�

E� immediato osservare che esistono k�� gradi di libert�a per la costruzionedelle Sk�x� infatti si devono determinare n�k � � incognite �i coe�cientidegli n polinomi di grado k de�niti sugli n intervalli dei nodi� mentre sipossono scrivere solo k�n � � � �n � � equazioni imponendo la regolarit�adella funzione e l�interpolazione�

Spline lineari

E� l�equazione della �spezzata� �funzione solo continua passante per gli n��punti dati� E� immediata la dimostrazione dell�esistenza ed unicit�a� la co�struzione pu�o essere e�ettuata banalmente utilizzando la tecnica di Lagrange�retta per due punti��

Consideriamo

S�j�x �

�������� x � xj��

x�xj��xj�xj�� � xj�� � x � xjxj���xxj���xj � xj � x � xj��

� x � xj��

il cui gra�co �e dato in �gura ���

0x x x x x x x n

1

S(x) S(x) S(x) S(x)

x

S(x)

0 j j+1 n

1 j-1 j j+1 j+2 n-1

Figura ����

si pu�o allora rappresentare S��x nella forma

S��x �nX

j��

f�xjS�j�x�

Page 91: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� �

Le S�j�x godono delle seguenti propriet�a

�� sono linearmente indipendenti� �la funzione identicamente nulla in �a� b si ottiene solo combinando con pesi tutti nulli le S�j�

��Pn

j�� S�j�x � �� �sono una partizione dell�unit�a�

�� in ogni intervallo �xj� xj�� ci sono solo due S�j�x diverse da zero�

�� la generica S�j�x �e diversa da zero solo in �xj��� xj�� �bi�spline linea�ri�

Per quanto riguarda l�errore d�interpolazione che si commette sostituendoalla f�x la spline lineare �e possibile dimostrare il seguente risultato

Teorema ���� Se f�x � C��a� b allora

jf�x� S��xj � h�

��f �� �x�� � �x � �a� b � ����

Dimostrazione ����� Dalla de�nizione di S�j�x� �x � �xj� xj�� si ha�ricordando l�errore nell�interpolazione in Lagrange�

f�x� S��x � �x� xj�x� xj��f �� ��

��

� � ��x � �xj� xj�� � per cui

jf�x� S��xj � h�j�

maxx

��f �� �x��

��

e� posto h � maxj

hj� allora

jf�x� S��xj � h�

��f �� �x�� ��

Dalla ����� al crescere del numero di nodi �quindi per h� �� si ottienela convergenza uniforme delle spline lineari�

Page 92: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

Spline cubiche

Per k � � si hanno le spline cubiche che possono essere costruite parten�do dall�interpolazione d�Hermite nel modo seguente� Dati gli n � � nodix�� x�� ���� xn �presi equidistanti� per sempli�care la notazione� si ha

�� S�i�x � !� �forma

�� S�i�xi�� � yi��� i � �� �� ���� n

�� S�i�xi � yi� �interpolazione

�� S ��i�xi � S ��i���xi� i � �� �� ���� n� �

� S��

�i�xi � S��

�i���xi� �regolarit�a

Postof ��xi � mi� h � xi � xi��

dove i valori mi NON sono noti� ma verranno determinati per costruire laspline� si ha� per l�interpolazione d�Hermite�

S�i�x �

�� �

h�x� xi��

��x� xi�h

��

yi�� � ����

��� �

h�x� xi

��x� xi��

h

��

yi �

��x� xi���x� xi�h

��

mi�� � �x� xi

�x� xi��

h

��

mi

dalle ���� segue che le S�i�x veri�cano le prime � condizioni date� per laquinta condizione� essendo

S��

�i�xi�� � �h�

���yi � yi��� h�mi � �mi��

S��

�i�xi � �h�

���yi�� � yi � h��mi � mi��

si ottiene il sistema lineare� di ordine n � � ����������������

�m� �m� � �h�y� � y�

m� ��m� �m� � �h�y� � y�

�m� ��m� �m� � �h�y� � y�

� � � � � � � � ����

�mn�� ��mn�� �mn � �h�yn � yn��

�mn�� ��mn � �h�yn � yn��

����

Page 93: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� INTERPOLAZIONE �M�F�� �

dove la prima e l�ultima equazione sono ottenute �spline naturali imponendo

S��

���x� � S��

�n���xn � ��

Il sistema ���� �e tridiagonale e fortemente diagonalizzato per cui esisteun�unica soluzione calcolabile con un numero di �ops dell�ordine di n� Per lespline con data tangente �m� e mn date� spline vincolate il sistema ���� siriduce a n� � equazione in n� � incognite �scaricando� m� e mn nel vettoretermini noti�

Una giusti�cazione della liscezza �smoothness delle spline cubiche natu�rali �e data dal seguente

Teorema ��� Fra tutte le funzioni f�x � C��a� b con f�xi � yi e xi ��a� b � la spline cubica naturale S��x interpolante �e tale per cuiZ b

a

hS

��

� �xi�dx �

Z b

a

hf��

�xi�dx��

Il teorema ����� pu�o essere parafrasato dicendo che la spline cubica natu�rale passa per tutti i punti �ssati �oscillando� il meno possibile�

Per quanto riguarda l�errore �e possibile dimostrare i seguenti risultati

Teorema ���� Se f�x � C��a� b la spline cubica interpolante S��x �e taleper cui

kf�x� S��xk � �

�h���f �� �x

�� ��

Teorema ��� Se f�x � C��a� b la spline cubica interpolante S��x �e taleper cui� per r � �� �� �� ����f �r �x� S

�r � �x

���� crh

��r ��f �� �x�� ��

Gli ultimi due teoremi non solo ci garantiscono la convergenza uniformedella spline S

�r � �x alla f �r �x� ma ci evidenziano come questa sia pi�u rapida

al crescere della regolarit�a di f�x �pi�u lenta al crescere dell�ordine delladerivata��

Page 94: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

���� Derivazione numerica

Con il termine di derivazione numerica si intende il calcolo della derivataprima �o superiore di un funzione f�x nota solo in un numero �nito di puntixi� Avendo calcolato il polinomio o la spline interpolanti �e naturale pensaredi utilizzare le loro derivate come stime delle derivate della funzione f�x�Mentre per le funzioni spline il teorema ����� ci garantisce la convergenzauniforme �sulle prime � derivate altrettanto non si pu�o dire per le formulederivanti dall�interpolazione polinomiale� In generale �e anzi vero il contrario�al crescere del numero dei nodi �tendere a zero della distranza fra due nodiconsecutivi� a parit�a d�errore sui dati� l�errore totale cresce�

In de�nitiva derivare numericamente con polinomi di grado elevato pu�oessere assai pericoloso� Per meglio comprendere la precedente a�ermazioneconsideriamo la seguente formula di derivazione �di�erenza centrale

f ��xi �f�xi��� f�xi��

�h� h�

�f �� ��

dove h � xi���xi e � � �xi��� xi�� � L�errore di troncamento che tale formulacomporta va a zero come h�� si pu�o per�o osservare che se i dati f�xi sononoti con un errore pre�ssato � nel calcolo di

f�xi��� f�xi���h

si commente un errore di arrotondamento dell�ordine di h

che �e sempremaggiore al diminuire di h� Pi�u precisamente posto M � max

��f �� ���� si

ha

E�h �

����f ��xi� f�xi��� f�xi���h

���� � �

h� M

h�

����

che �e minimo� come �e facile veri�care derivando rispetto ad h� quando

h � h� ��

r��

M� �����

Il valore h� prende il nome di h�ottimo in quanto minimizza l�errore totaledato dalla ���� Dalla ����� si vede bene che se � �e grande �errori grandi suidati anche h deve essere preso grande �punti �distanti� per ridurre l�e�ettodell�errore d�arrotondamento �anche se l�errore di troncamento crescer�a��

Page 95: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� MINIMI QUADRATI �M�F��

Lasciamo al lettore la veri�ca dell�errore di troncamento introdotto dalleseguenti formule di derivazione

f ��xi �f�xi��� f�xi

h� h

�f �� ��

f ��xi �f�xi� f�xi��

h�h

�f �� ��

f��

�xi �f�xi��� �f�xi � f�xi��

h�� h�

��f �� ��� �����

Altre formule di derivazione� anche per derivate di ordine pi�u elevato�possono essere reperite in letteratura�

�� Minimi quadrati �M�F��

Analogamente a quanto fatto nel paragrafo �� � relativamente alla risoluzionedi sistemi lineari nel senso dei minimi quadrati� ci proponiamo ora di costruireil polinomio di grado m approssimante� nel senso dei minimi quadrati� unadata funzione� Questo polinomio �e pi�u generale del polinomio interpolante �idue coincidono quando m coincide con il numero n dei punti d�interpolazioneed esiste anche se m � n�

Siano dati gli n punti �xi� yi i � �� �� ���� n� dove yi � f�xi �e il valoreassunto dalla funzione f in xi� consideriamo il polinomio di grado m

pm�x �mXj��

jxj� m � n

e cerchiamo gli m � � valori j che minimizzano

J� �nXi��

�f�xi� pm�xi � �

nXi��

�f�xi�

mXj��

jxji

� �����

Derivando la ����� rispetto ad k �k � �� �� ���� m ed uguagliando a zerole m � � equazioni che si ottengono si ha

�J�

�k� ��

nXi��

��f�xi�

mXj��

jxji

xki

�� �� k � �� �� ���� m�

Page 96: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

che� scritto in forma matriciale diviene�������������������������

n� ��

Pi

xi � � � �m

Pi

xmi �Pi

yi

Pi

xi ��

Pi

x�i � � � �m

Pi

xm��i �

Pi

yix�i

Pi

x�i ��

Pi

x�i� � � �m

Pi

xm��i �

Pi

yix�i

������

� � ����

����

Pi

xm��i ��

Pi

xmi � � � �m

Pi

x�m��i �Pi

yixm��i

Pi

xmi ��

Pi

xm��i � � � �m

Pi

x�mi �Pi

yixmi

� �����

De�niamo la matrice �di Vandermonde rettangolare

V �

������ � x� x�� � � � xm�� x� x�� � � � xm����

������

������

� xn�� x�n�� � � � xmn��� xn x�n � � � xmn

������� �����

ed i vettori

������ �

����

m��m

������� � y �

������ y�y����

yn��yn

�������il sistema ����� si pu�o scrivere nella forma

V TV � V Ty�

�si confronti il paragrafo �� sui sistemi sovradeterminati�

Osservazione ���� La matrice ���� �e stata ottenuta considerando la basedi polinomi xk k � �� �� ���� m� che compaiono nella rappresentazione dipm�x �

Pmj�� jx

j� valutati nei nodi xi i � �� �� ���� n� Se si cambia la baseo se si considerano di�erenti funzioni approssimanti� per esempio splinebasta cambiare la matrice ����� restando tutto il resto inalterato�

Page 97: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� MINIMI QUADRATI �M�F�� �

����� Polinomi trigonometrici

Consideriamo brevemente il caso dei polinomi trigonometrici di Fourier ilcui utilizzo� ed importanza nelle applicazioni� non �e certo il caso di ricordarequi�

Data una funzione f�x� periodica in ��� �� � in n nodi equidistanti xi�i � �� �� ���� n si vuole determinare il polinomio trigonometrico

Fm�x �a��

�mXj��

aj cos�jx �mXj��

bj sin�jx

che minimizza

J�aj� bj �nXi��

�f�xi� Fm�xi � �

�nXi��

�f�xi�

�a��

�mXj��

aj cos�jx �mXj��

bj sin�jx

� �

Per quanto precisato nell�osservazione ����� il calcolo dei coe�cienti aj ebj richiede la risoluzione� nel senso dei minimi quadrati� del sistema lineare

F TFc � F Ty ����

dove F � c ed y sono dati da

F �

������ � cos�x� � � � cos�mx� sin�x� � � � sin�mx�� cos�x� � � � cos�mx� sin�x� � � � sin�mx����

������

������

� cos�xn�� � � � cos�mxn�� sin�xn�� � � � sin�mxn��� cos�xn � � � cos�mxn sin�xn � � � sin�mxn

�������

c �

����������

a��

a����amb����bm

������������ y �

��������

y�y�y����

yn��yn

����������

Page 98: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

La risoluzione del sistema ���� � se i nodi xi sono equidistanti� non ri�chiede� come ci si potrebbe aspettare� un numero di operazioni dell�ordinen� �mn �ops� ma solo mn� in quanto la matrice F risulta �quasi�ortogonale�per cui basta calcolare il vettore F Ty�

Osservazione ���� Si veri�chi� utilizzando il programma MATLAB� che

F TF � diag�m�m

�� ����

m

��

Osservazione ���� Se si cerca il polinomio di Fourier interpolante m �n la matrice F risulta quadrata per cui basta risolvere il sistema

Fc � y�

che� sfruttando le propriet�a di F pu�o essere risolto in soli n log� n �opstrasformata veloce di Fourier FFT�

� Riepilogo

Sintetizzando quanto visto nel presente capitolo possiamo dire che

�� L�interpolazione polinomiale di Lagrange �e utile quando si hanno pochipunti d�interpolazione�

�� Le di�erenze divise sono pi�u e�cienti computazionalmente quando sidevono aggiungere dei nodi �si salvano i calcoli gi�a fatti�

�� L�nterpolazione d�Hermite �e utile quando si hanno informazioni anchesulle derivate della funzione�

�� Le spline sono preferibili quando si devono ricostruire linee �lisce��computer graphics� CAD ecc��� sono molto accurate se la funzioned�approssimare �e regolare�

� La derivazione numerica �e pericolosa in quanto gli errori d�arrotonda�mento crescono al tendere a zero della distanza fra i nodi�

�� I minimi quadrati sono utili per individuare �modelli� quando i datisono a�etti da errore� La scelta delle funzioni di base dipendono dalmodello� mentre i coe�cienti lo caratterizzano�

Page 99: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

�� RIEPILOGO

�� L�utilizzo dell�approssimazione mediante polinomi di Fourier� indipen�dentemente dalle innumerevoli applicazioni pratiche� �e anche computa�zionalmete poco oneroso grazie alle peculiarit�a delle funzioni di base eall�algoritmo FFT�

����� Esercizi

Esercizio ���� Costruire il polinomio interpolante i dati

xi � � � �yi ���� ��� ���� ���� �

Esercizio ���� Sapendo che un certo fenomeno si evolve nel tempo con laseguente legge

f�t � at � b log t

determinare� nel senso dei minimi quadrati� i parametri a e b sapendo che

t � �� �� �f�t ��� � ����� ����� � �� �

Esercizio ���� Dati tre punti nel piano� non allineati� costruire la cubicapassante per il primo ed il terzo ed ivi tangente alle � rette uscenti dal secondopunto� E� unica � Perch�e �

Esercizio ��� Scrivere un programma in MATLAB che risolva l�esercizioprecedente e disegni i � punti� le � rette e la cubica� Meglio se il programmaprevede l�input gra�co dei � punti�

Esercizio ���� Ricavare la formula di derivazione �����

Page 100: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO � TEORIA DELL�APPROSSIMAZIONE �M�FRONTINI�

Page 101: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo

Formule di Quadratura�M�Frontini�

Data una funzione f�x continua nell�intervallo �a� b dovendo calcolare

I�f �

Z b

a

f�xdx � ��

�e ben noto �teorema fondamentale del calcolo integrale che esiste F �x �C��a� b funzione primitiva di f�x �F ��x � f�x� ed inoltre

I�f � F �b� F �a� � ��

E� interessante osservare che mentre la � �� rappresenta un numero reale�la � �� a�erma che tale numero pu�o essere determinato se si conosce unadelle in�nite funzioni primitive di f�x� E� chiara la di�erenza di �quantit�adi informazioni� in gioco � �La conoscenza di un numero contro la conoscenzadi una� anzi in�nite� funzioni���

Non sempre la funzione primitiva� che pur esiste se f�x �e continua� �enota o �e facilmente computabile� per cui ci si propone il compito di calcolareI�f utilizzando solo le informazioni date in � �� �intervallo d�integrazione�a� b e funzione integranda f�x�

Le formule di quadratura che considereremo� ricavabili dall�interpolazionedi Lagrange� sono della forma

I�f � Qn�f � R�Qn� f �nX

j��

A�n j f�xj � R�Qn� f � ��

���

Page 102: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

dove gli n pesi A�n j e gli n nodi xj � �a� b � della formula di quadratura Qn�f�

sono scelti in modo che il resto R�Qn� f sia �piccolo� �in dipendenza dellaregolarit�a di f�x e dal numero di nodi n�

��� Formule di Newton�Cotes �M�F��

Dall�interpolazione di Lagrange abbiamo che� �ssati i nodi xi � �a� b � �epossibile costruire il polinomio interpolante pn�x tale per cui

f�x � pn�x � en�x � ��

dove en�x �e l�errore d�interpolazione�Integrando la � �� si ottieneZ b

a

f�xdx �

Z b

a

pn�xdx �

Z b

a

en�xdx

essendo

pn�x �nX

j��

f�xjlj�x

dove lj�x sono i polinomi di Lagrange sui nodi xi� si ha

Qn�f �nX

j��

f�xj

Z b

a

lj�xdx �nX

j��

A�n j f�xj

dove

A�n j �

Z b

a

lj�xdx � �

sono i pesi della formula di quadratura e

R�Qn� f �

Z b

a

en�xdx

�e l�errore che si commette sostituendo ad I�f il valore Qn�f�

Osservazione ����� Dalla ��� �e immediato osservare che i pesi dipendo�no dai nodi d�interpolazione zeri della formula di quadratura ed inoltre laformula sar�a esatta se la f�x �e un polinomio di grado al pi�u n�

Page 103: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� FORMULE DI NEWTON�COTES �M�F�� ���

De�nizione ����� Una formula del tipo ��� si dice di ordine n se �e esattaper polinomi di grado al pi�u n�

Dall�osservazione ���� e dalla precedente de�nizione segue che le formuledi quadratura interpolatorie si possono costruire utilizzando il concetto diordine di precisione determinando i pesi in modo che l�ordine risulti massimo�

Posto

Qn�f �nX

j��

A�n j f�xj

consideriamo i monomi xk �k � �� �� ���� n ed imponiamo che

I�xk �

Z b

a

xkdx ��

k � ��bk�� � ak�� �

nXj��

A�n j xkj � Qn�xk

otteniamo il sistema lineare� nelle incognite A�n j �

nXj��

A�n j xkj �

k � ��bk�� � ak��� � ��

Esempio ����� Dato l�intervallo ��� h � posto x� � � e x� � h� costruiamola formula del primo ordine�

Avremo

Q��f � A�f�x� � A�f�x�

da cui � R h

�dx � h � A� � A�R h

�xdx � h�

�� A�� � A�h

e quindi �A� � A� � h

A�h � h�

�� A� � A� �h

che �e la formula dei Trapezi

Q��f �h

��f�x� � f�x� �� � ��

Page 104: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

Analogamente �e possibile costruire la formula di Simpson su tre nodiequidistanti nell�intervallo ��h� h �x� � �h� x� � �� x� � h

Q��f �h

��f�x� � �f�x� � f�x� �� � ��

La formula dei trapezi poteva essere ottenuta� mediante l�interpolazionedi Lagrange� costruendo la retta passante per i punti ��� f�� e �h� f�h�

Ricordando che� in questo caso� l�errore d�interpolazione �e

e��x � f�x� p��x �x�x� h

�f��

��� � � ��x � ��� h

abbiamo

R�Q�� f �

Z h

e�xdx �

Z h

x�x� h

�f��

��dx�

Applicando il teorema della media �essendo x�x�h di segno costante in��� h si ha

R�Q�� f �

Z h

x�x� h

�f��

��dx � f��

��

Z h

x�x� h

�dx �

�x�

�� hx�

�h�

f��

�� � �h�

��f��

���

che ci fornisce l�errore per la formula dei trapezi�

Osservazione ����� La presenza� nella formula dell�errore� della derivataseconda della funzione integranda riconferma che la formula �e esatta perpolinomi di primo grado�

Per ottenere l�errore nella formula di Simpson osserviamo cheZ h

�hx�dx � � �

h

��h� � � � � � h��

per cui la formula risulta esatta �sorpresa� anche per polinomi di terzogrado�

Se applichiamo la formula di Simpson al monomio x�� abbiamoZ h

�hx�dx �

h �� h

���h� � � � � � h�

��

�h

Page 105: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� FORMULE DI NEWTON�COTES �M�F�� ��

per cui� in questo caso� si commette un errore pari a � ��h�

Se la funzione f�x �e su�cientemente regolare possiamo considerarne losviluppo di Mac�Laurin arrestato al quarto ordine� per cui

f�x � f�� � xf ��� �x�

�f ���� �

x�

��f���

�� �x�

��f iv��

ed integrando termine a termine si commetter�a errore solo sull�ultimo adden�do �i primi � termini sono polinomi di grado minore od uguale al terzo� percui

R�Q�� f � � �

� h

��f iv�� � � �

�hf iv���

Osservazione ����� Il guadagno di un ordine di precisione con � nodi siha ordine � e non � non �e peculiare del metodo di Simpson ma �e proprio ditutte le formule di Newton�Cotes su un numero dispari di nodi equidistanti�

A sostegno della precedente osservazione ricordiamo che per la formuladei rettangoli �si ricordano gli �scaloidi � nella de�nizione di integrabilit�asecondo Reimann� sull�intervallo ��h� h

Q��f � �hf��

comporta un errore

R�Q�� f � �h�

�f��

�� � �

per cui risulta esatta non solo per costanti� ma anche per rette� Una chiarainterpretazione geometrica di questo risultato �e dato dalla seguente �gura ���La scelta� del tutto arbitraria� dell�intervallo ��� h piuttosto che ��h� h � nelladeterminazione delle formule precedenti� �e giusti�cata dal fatto che �e semprepossibile passare� con un opportuno cambio di variabile� da un intervallod�integrazione ad un altro� Analogamente �e possibile� nota una formula diquadratura �nodi e pesi su un certo intervallo� ottenere l�equivalente formulasu un altro intervallo�

Supponiamo di dover calcolare

I�f �

Z b

a

f�xdx

avendo a disposizione la formula di quadratura

Qn�f �nX

j��

Ajf�tj �Z ��

��f�tdt� � ���

Page 106: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

0 0.2 0.4 0.6 0.8 1 1.20

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Punto medio

Trapezi

Uguale Area

Uguale Area

Figura ���

Come �e possibile ottenere la formula

Q�n�f �

nXj��

A�jf�xj �Z b

a

f�xdx � ���

dalla � ��� � Basta considerare il cambiamento di variabile

x �b� a

�t �

b � a

�� ���

per cui si ha Z b

a

f�xdx �

Z ��

��f�b� a

�t �

b � a

�b� a

�dt

� b� a

nXj��

Ajf�tj �nX

j��

A�jf�xj

e quindi i pesi della Q�n�f sono dati da

A�j �b� a

�Aj

mentre i nodi xj� ottenuti dalla � ��� sono

xj �b� a

�tj �

b � a

���

Page 107: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� FORMULE DI NEWTON�COTES �M�F�� ���

Esistono due tipi di formule di Newton�Cotes

�� formule di tipo chiuso� comprendono fra i nodi anche gli estremi del�l�intervallo d�integrazione �vedi Trapezi e Simpson�

�� formule di tipo aperto� utilizzano solo nodi interni all�intervallo d�inte�grazione �vedi punto medio�

Le formule di tipo aperto sono utili� per esempio� quando si hanno singo�larit�a agli estremi�

Per formule di Newton�Cotes di ordine pi�u elevato si rimanda alla biblio�gra�a�

����� Formule composite

Le stime dell�errore fornite sono state ottenute considerando� �ssato l�inter�vallo d�integrazione ed il numero di nodi� il polinomio interpolante globale�Se �e vero che� al crescere del numero dei nodi� cresce l�ordine della derivatadella f�x che compare nella formula dell�errore �maggiore ordine d�accura�tezza� �e altres�i vero che si ha anche una potenza crescente dell�intervallod�integrazione �che �e per�o �ssato� quindi costante� Come nell�interpolazionesi �e passati� per motivi di convergenza� dal polinomio interpolante globale aipolinomi a tratti �spline cos�i ora converr�a suddividere l�intervallo d�integra�zione �a� b �additivit�a dell�integrale in tanti sotto intervalli di ampiezza hed applicare� su ogni sotto intervallo� la formula prescelta �convergenza perh� �� E� l�applicazione del vecchio principio di Cesare Augusto �divide etimpera��

Si ottengono cos�i le formule composite �o generalizzate� di cui proveremola convergenza nel caso dei trapezi� Consideriamo

I�f �

Z b

a

f�xdxadditivit�a

�m��Xk��

Z xk��

xk

f�xdx �m��Xk��

Ik�f

dove

x� � a� xm � b� h �b� a

m

Page 108: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

e

Ik�f �

Z xk��

xk

f�xdx �nX

j��

Ajf�xj � Qn�f � Ik�

con Qn�f � Ik si �e indicata la formula Qn�f� di ordine n� applicata sull�in�tervallo Ik�

Per la formula dei trapezi� applicata a Ik� abbiamoZ xk��

xk

f�xdx �h

��f�xk � f�xk�� � h�

��f ����k

da cui la formula composta

Tm�f �h

m��Xk��

�f�xk � f�xk�� �

Per l�errore abbiamo quindi

ET � I�f� Tm�f � �m��Xk��

h�

��f ����k� � ���

che pu�o essere sempli�cata nella forma

ET � �f ��

��m��Xk��

h�

��� �mh�

��f��

�� � ��b� ah�

��f��

���� � ���

Per ottenere dalla � ��� la � ���� abbiamo fatto uso del seguente

Lemma ����� Data una funzione g�x � C��a� b e dei coe�cienti di segnocostante ak� allora esiste un � � �a� b tale per cui� �xk � �a� b � si ha

m��Xk��

akg�xk � g��m��Xk��

ak��

Osservazione ���� Si osservi che l�errore nella formula dei trapezi com�posita �e quindi della forma ch� dove

c � ��b� a

��f��

��

e va a zero� quando h� �� come h� convergenza quadratica�

Page 109: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� FORMULE ADATTIVE �M�F�� ��

Analogamente per la formula di Simpson composita si pu�o dimostrare chel�errore �e della forma Ch�� infatti� postoZ xk��

xk

f�xdx �h

��f�xk � �f�xk�� � f�xk�� � h

�f iv��k

si ha

Sm�f �m��Xk��

h

��f�x�k � �f�x�k�� � f�x�k��

dove

x� � a� xm�� � b� h �b� a

�m

per cui� sempre in virt�u del Lemma �����

ES � I�f� Sm�f � �m��Xk��

h

�f iv��k � ��b� a

h�

���fiv

���� � ��

Osservazione ����� Dalle ��������� si sarebbe tentati di prendere n �in��nitamente grande� per avere h �in�nitamente piccolo�� ed ottenere quindiun errore �trascurabile�� Nella pratica NON ha senso prendere n �troppogrande� in quanto� oltre ad aumentare il tempo di calcolo troppe valutazionidella funzione integranda� si rischia che� da un certo punto in poi� gli erroridi calcolo epsilon macchina non facciano guadagnare in precisione si pu�oanzi peggiorare il risultato��

��� Formule adattive �M�F��

In tutte le stime dell�errore fornite �e implicitamente sottointesa una certaregolarit�a della funzione integranda f�x �devono almeno esistere le derivatein gioco nelle varie formule ma non si dice quanto valgono le costanti c e Cche dipendono dal punto �� per cui alla domanda� �quanti punti si devonoprendere per avere un errore pre�ssato� � Non possiamo dare una rispostaa priori� senza conoscere ��

Sovente non �e possibile determinare preventivamente � ��e di�cile� se nonimpossibile� maggiorare la derivata di f�x che compare nella formula dell�er�rore� in questi casi sono utili le formule adattive �note anche come integratoriautomatici o programmi con controllo automatico dell�errore�

Page 110: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

Chiariamo con un semplice esempio di cosa si tratta�L�idea di base �e quella di fornire� se possibile� noti la funzione integranda

f�x e l�intervallo d�intregrazione �a� b � il valore dell�integrale I�f con unapre�ssata precisione �toll� Cercheremo di fare questo seguendo le seguentidue linee

�� usare l�errore locale per controllare l�errore globale�

�� ridurre al minimo le valutazioni della funzione integranda�

Esempli�chiamo tale processo mediante la ben nota formula di Simpson�Utilizzando il metodo di Simpson sull�intero intervallo �a� b � si haZ b

a

f�xdx� S��f � � �

�b� a

fiv

�� � ���

mentre� dividendo in due parti l�intervallo� si haZ b

a

f�xdx� S��f � ���

�b� a

fiv

�� � ���

dove � �e� in generale� diverso da �� Supponendo� e questa �e una condizioneforte se �a� b �e grande� mentre �e �ragionevole�� se f�x �e abbastanza regolaree �a� b �e �piccolo�� che

fiv

�x � costante in �a� b �

posto

I�f �

Z b

a

f�xdx

si ottiene� dalle � ���� ���� che

I�f� S��f � �

���I�f� S��f

per cui�� �I�f� S��f � I�f� S��f�

Sottraendo ad entrambi i menbri la quantit�a I�f� S��f� si ha

� �I�f� S��f � S��f� S��f

Page 111: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� FORMULE ADATTIVE �M�F�� ���

e quindi

I�f� S��f � S��f� S��f

� � � ���

Nella � ��� l�errore commesso intregrando con la formula di Simpson su� sotto intervalli �e controllato dalla quantit�a a secondo membro che� essendolegata alla di�erenza fra due integrazioni numeriche� �e computabile�

Osservazione ����� Volendo essere pessimisti� e quindi salvaguardarsi nel�la stima dell�errore� la costante �

�che compare nella ���� pu�o essere sosti�

tuita da ���

Un algoritmo che sintetizza quanto detto pu�o essere cos�i formulato�

�� si parte dall�intervallo �a� b e con la tolleranza toll�

�� si calcolano S��f e S��f�

�a se���S��f �S��f �

��� � tolleranza� allora S��f �e uguale all�integrale

nella precisione richiesta�

�� altrimenti si dimezzano l�intervallo e la tolleranza�

�a se NON si sono fatte �troppe� sottodivisioni� si opera analoga�mente sul nuovo intervallo�

�b altrimenti l�integrale non pu�o essere calcolato con la tolleranza tollrichiesta�

�� calcolato� nella nuova precisione� l�integrale sul sotto intervallo� si tornaal punto � con l�intervallo restante e ridotta tolleranza��

Quanto visto per la formula di Simpson pu�o essere esteso� con le dovutevarianti� ad altre formule� Non potendoci dilungare oltre si rimanda il lettoreinteressato alla bibliogra�a speci�ca�

Quando la funzione f�x ha derivata seconda di segno costante in �a� b �supponiamola per esempio negativa si pu�o fornire una stima dell�errore�che si commette integrando numericamente la f�x in �a� b � utilizzando leformule del punto medio e dei trapezi� E� facile osservare �si confronti la�gura �� che

Page 112: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

0 0.2 0.4 0.6 0.8 1 1.20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

f(x)

Trapezi

Punto medio

Figura ���

IT �f � I�f � IMP �f�

dove con IT �f si �e indicato il valore ottenuto con la formula dei trapezi econ IMP �f quello ottenuto con la formula del punto medio�

Utilizzando le formule composite� suddividendo l�intervallo �a� b in sottointervalli di ampiezza h� si ha� con ovvio signi�cato dei simboli�

IhT �f � I�f � IhMP �f

da cui

� � IhT �f� IhT �f � I�f� IhT �f � IhMP �f� IhT �f

ed inoltre

� �IhMP � IhT �f

� � I�f� IhMP �f � IhMP �f� IhMP �f � �

sommando� termine a termine� la prima relazione e la seconda� si ha

� �IhMP � IhT �f

� � �I�f� �IhT �f � IhMP �f

� � IhMP �f� IhT �f

per cui ����I�f� IhT �f � IhMP �f

���� � IhMP �f� IhT �f

��

dove la quantit�a a secondo membro �e computabile�

Page 113: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� FORMULE GAUSSIANE �M�F�� ���

�� Formule gaussiane �M�F��

Se nella formula � �� i nodi xj non sono �ssati a priori� �equidistanti nelleformule di Newton�Cotes imponendo che la formula abbia il massimo ordinedi accuratezza si ottiene un sistema non lineare nelle �n � � incognite xj eAj� Gauss fu il primo a dimostrare che� sotto opportune ipotesi� tale sistemaammette una ed una sola soluzione� dando origine ad uno dei pi�u proli�cisettori della matematica�

Alla base delle formule di quadratura gaussiane sta la teoria dei polino�mi ortogonali� Non potendo qui dedicarle il dovuto spazio ci limiteremo ariassumere solo alcuni degli aspetti pi�u signi�cativi relativi alle formule diquadratura�

De�nita con

Q�n�f �nX

j��

A�n j f�xj � ��

la formula di quadratura con n nodi liberi� vale il seguente

Teorema ����� Il sistema non lineare ��� dove i nodi xj sono incogniti�con �a� b � ���� � � ammette una ed una sola soluzione ed inoltre nella Q�n�f

�� i pesi A�n j sono positivi�

�� i nodi xj sono gli zeri del polinomio ortogonale di Legendre di grado nLn�x�

�� tale formula di quadratura gaussiana �e esatta per ogni polinomio digrado minore o uguale a �n � ���

Che la formula Q�n�f non possa essere di grado maggiore di �n � � �efacilmente dimostrabile considerando il polinomio

��n�x �nY

j��

�x� xj�� �xj zeri di Ln�x

in quanto il suo integraleZ �

����n�xdx �

Z �

��

nYj��

�x� xj�dx � �

Page 114: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

�e una quantit�a strettamente maggiore di zero� mentre verr�a valutato ugualea zero dalla � ��� �Si osservi che �e nullo il valore del polinomio ��n�x intutti i nodi della formula di quadratura�

La positivit�a dei pesi A�n j �e importante perch�e permette d�assicurare

la convergenza� al crescere del numero dei nodi �n � �� delle formulegaussiane in virt�u del seguente

Teorema ����� Se f�x � C��a� b data una formula di quadratura dellaforma

Qn�f �nX

j��

A�n j f�xj

se esiste una costante k tale per cui

nXj��

���A�n j

��� � k� �n

alloralimn

Qn�f � I�f��

Osservazione ����� EssendoPn

j��A�n j � b � a e A

�n j � � per ogni j�

baster�a prendere� nel teorema ������ k � b� a�

Osservazione ����� Il teorema ����� non pu�o essere applicato nel caso delleformule di Newton�Cotes perch�e� per n � �� compaiono pesi A

�n j negativi�

L�esistenza delle formule gaussiane �e legata all�esistenza di una famiglia dipolinomi ortogonali sull�intervallo �non necessariamente �nito �a� b rispettoad una opportuna funzione peso �x� non negativa in �a� b �

De�nizione ����� Si dicono ortogonali in �a� b rispetto alla funzione peso �x� non negativa in �a� b i polinomi pn�x che veri�canoZ b

a

pi�xpj�x �xdx � ci�ij� ci � ���

Se ci � � per ogni i� i polinomi si dicono ortonormali�

Osservazione ����� Nel caso dei polinomi Ln�x di Legendre� �a� b � ���� � e �x � ��

Page 115: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� FORMULE GAUSSIANE �M�F�� ��

Vale il seguente

Teorema ����� Se esiste la famiglia di polinomi ortogonali in �a� b rispettoalla funzione peso non negativa �x allora� dato

I�f �

Z b

a

f�x �xdx � ���

esiste la formula di quadratura gaussiana

Q�n�f �nX

j��

A�n j f�xj � ���

dove

A�n j �

Z b

a

l�j �x �xdx

dove l�j �x �e il j�esimo polinomio di Lagrange di grado n sui nodi xk k ��� �� ���� n ed i nodi xk sono gli zeri interni all�intervallo �a� b del polinomioortogonale di grado n��

Osservazione ���� Si fa osservare che mentre la funzione integranda in���� comprende anche la funzione peso �x� nella ���� compare solo laf�x valutata nei nodi� Questo ci permette di rimuovere eventuali singolarit�anella funzione integranda� scaricandole sulla funzione peso �x�

����� Integrali impropri

Per concludere il capitolo faremo solo un breve cenno al problema del calcolodegli integrali della forma

I��f �

Z �

f�xdx� I��f �

Z �

f�xpxdx

�integrali impropri dove f�x e g�x � f�x px

sono funzioni integrabili rispetti�

vamente in R e in ��� � � Vedremo che in questi casi si possono ancora utilizzarele formule di Newton�Cotes �eventualmente aperte� anche se le formule diGauss risultano in generale pi�u e�cienti�

Page 116: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

Per calcolare l�integrale I� con pre�ssata precisone ��� si pu�o operare�sfruttando l�additivit�a dell�integrale� nel modo seguente� si pone

I��f �

Z �

f�xdx �

Z a

f�xdx �

Z �

a

f�xdx � Ia�f � I�f

e si sceglie a in modo che ����Z �

a

f�xdx

���� � ��

�tale a deve esistere per l�integrabilit�a di f�x e si calcola Ia�f con unaopportuna formula di Newton�Cotes Qn�f in modo che

jIa�f�Qn�fj � ��

dove �� e �� sono stati scelti in modo che �� � �� � ��Per calcolare l�integrale I� �e possibile utilizzare la formula del punto medio

composita in modo da evitare il calcolo della funzione integranda in x � ��singolarit�a� Se la funzione g�x �e su�cientemente regolare la formula delpunto medio converge all�integrale I�� anche se tale convergenza pu�o essereassai lenta�

Osservazione ����� Dall�osservazione ����� segue che basta costruire la fa�miglia di polinomi ortogonali in ��� � rispetto alla funzione pesi �x � �p

x

per ottenere una formula di quadratura gaussiana che permette di calcolarein modo esatto I� ogniqualvolta f�x �e un polinomio� Di fatto tale formulaesiste ed �e correlata ai polinomi di Chebyshev�

Per una pi�u esauriente trattazione dell�uso delle formule gaussiane sirimanda alla bibliogra�a speci�ca�

�� Riepilogo �M�F��

Sintetizzando quanto esposto nel presente capitolo abbiamo

�� Se l�intervallo �a� b �e �nito� e f�x �e regolare le formule di Newton�Cotescomposite funzionano in generale bene�

Page 117: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� RIEPILOGO �M�F�� ���

�� Si possono avere problemi coi tempi di calcolo se si richiede moltaaccuratezza �numero eccesivo di sottodivisioni dell�intervallo�

�� La convergenza� per le formule di Newton�Cotes� �e garantita al cresceredel numero di nodi �composite non al crescere dell�ordine �grado delpolinomio interpolante�

�� Volendo ottenere un risultato �con una pre�ssata precisione� �e meglioutilizzare formule adattive �integratori automatici� anche se in generalei programmi sono pi�u complessi�

� Se l�intervallo �a� b �e in�nito si possono ancora usare le Newton�Cotespur di spezzare opportunamente l�integrale �additivit�a�

�� Le formule gaussiane sono pi�u accurate� convergenti e permettono ditrattare� in alcuni casi� gli integrali singolari� �Solo cenni nel corso�

����� Esercizi

Esercizio ���� Dimostrare che il sistema ��� �e non singolare se i nodi xisono distinti�

Esercizio ���� Costruire la formula ����

Esercizio ���� Dimostrare la formula ���� Suggerimento� costruito ilpolinomio di Lagrange su un punto� integrando l�errore d�interpolazione������

Esercizio ��� Veri�care che la formula del punto medio �e gaussiana Gauss�Legendre con n � ��

Esercizio ���� Calcolare� se possibile� con errore relativo minore di ���� iseguenti integrali

I� �

Z �

e�x���dx� I� �

Z �

e�x���

xdx� I� �

Z �

cos x

x � �dx�

Page 118: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� CAPITOLO �� FORMULE DI QUADRATURA �M�FRONTINI�

Page 119: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Capitolo �

Equazioni di�erenziali ordinarie�M�Frontini�

Molti problemi della �sica� dell�ingegneria e delle discipline scienti�che ingenerale sono governati da equazioni di�erenziali ordinarie� come esempio sipensi alle equazioni

m��x� f �equazione della dinamica

Ri � Ldi

dt� vc � Vs �equazione circuito RLC

che permettono rispettivamente di studiare� date le opportune condizioniiniziale� il moto di un corpo di massa m soggetto ad una forza f � ed il varia�re� nel tempo� dell�intensit�a di corrente �i e del voltaggio del condensatore�vc in un circuito elettrico con una resistenza �R� un�induttanza �L� uncondensatore �di capacit�a C� con i � C dvc

dt ed un generatore �Vs�

Non sempre �e possibile determinare in �forma chiusa�� mediante i meto�di forniti nei corsi di Analisi Matematica� la soluzione di tali problemi� edanche quando ci�o �e possibile tale soluzione pu�o risultare tanto complicatada non essere facilmente utilizzabile� Va notato inoltre che spesso in molteapplicazioni non �e necessario conoscere la soluzione ovunque ma basta averneuna �buona stima� solo in alcuni punti� Per questi motivi i metodi numericirivestono un ruolo sempre pi�u importante nel calcolo scienti�co�

Prima di addentrarci nella presentazione dei metodi numerici richiamiamoalcuni risultati classici nella teoria delle equazioni di�erenziali ordinarie�

��

Page 120: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

Dato un problema ai valori iniziali �problema di Cauchy nella forma�y� � f�x� yy�a � y�

����

dove f�x� y esprime il legame fra una funzione e la sua derivata prima�mentre la condizione iniziale y�a � y� impone il passaggio della soluzioneper il punto �a� y�� vale il seguente

Teorema � �� Se la funzione f�x� y del problema ��� �e de�nita e conti�nua in a � x � b per ogni y ed inoltre veri�ca la condizione di Lipschitz

jf�x� u� f�x� vj � L ju� vj

a � x � b� �u� v ed L � costante�

allora esiste ed �e unica la soluzione del problema� cio�e

��y�xjy� � f�x� y con y�a � y���

Corollario � �� Condizione su�ciente per l�esistenza ed unicit�a della so�luzione di ��� �e che esista la derivata parziale di f rispetto ad y ed inoltre�

jfy�x� yj � L� a � x � b� �y��

Oltre a richiedere l�esistenza e l�unicit�a della soluzione del problema �����e importante che il problema ���� sia ben posto �nel senso di Hadamard ov�vero che� la soluzione� che esiste ed �e unica� dipenda con continuit�a dai dati�Il concetto di problema ben posto �e importante nelle applicazioni perch"e cipermette di risolvere un problema partendo anche da condizioni a�ette da er�rori ��non troppo grandi� ed ottenere comunque una soluzione �abbastanzavicina� alla soluzione vera del problema�

Esempio � �� Come esempio di problema ben posto si consideri l�equazionedi�erenziale

y� � y

che ammette come soluzione la famiglia di curve

y�x � cex� c � R�

Page 121: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���

Se si considera il corrispondente problema di Cauchy con la condizioneiniziale

y�� � �

si ottiene l�unica soluzioney�x � ex�

Il problema �e ben posto� infatti� se si considera il problema perturbato�y� � yy�� � � � �

si ottiene la soluzioney�x � �� � � ex

che� per � piccoli� poco si discosta� dalla soluzione y�x � ex si valuti l�errorerelativo��

Per trattare equazioni di�erenziali di ordine superiore al primo� ricordia�mo che un�equazione di ordine m

y�m � f�x� y� y�� y��� ���� y�m�� �

pu�o essere scritta sotto forma di sistema di m equazioni del primo ordine�����������

z���x � z��xz���x � z��x���z�m���x � zm�xz�m�x � f�x� z�� z�� ���� zm

dove si �e posto

z��x � y�x� z��x � y��x� z��x � y���x� ���� zm�x � y�m�� �x

per cui� date le opportune condizioni iniziali� siamo ricondotti ad un problemadi Cauchy �

z� � f�x� zz�a � z�

formalmente equivalente al sistema �����

Fatte salve le ipotesi dei teoremi precedenti� ci proponiamo ora di forniremetodi numerici che� partendo dal valore iniziale dato nel problema �����forniscano la soluzione in pre�ssati punti xi di un dato intervallo �a� b �

Page 122: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

�� Metodi ad un passo �M�F��

I metodi ad un passo sono metodi numerici che� �ssato un passo d�integrazio�ne h� determinano la soluzione nel punto x� � a� h� x� � a� �h���� partendodal valore noto in x� � a� mediante relazioni della forma

yn�� � yn � h��xn� yn� h� n � �� �� �� ���

dove con yn si �e indicata la stima numerica di y� xn� mentre ��xn� yn� h �euna funzione che dipender�a da f�xn� yn ed h� �I valori trovati al passo n� �dipendono solo dai valori ottenuti al passo precedente� metodi ad un passo�

Esistono due grandi famiglie di metodi ad un passo di cui parleremo nelseguito

�� i metodi di Taylor�

�� i metodi di Runge�Kutta�

����� Metodi di Taylor

I metodi di Taylor si basano sull�utilizzo dello sviluppo in serie di Taylor dellafunzione incognita y�x� pi�u precisamente dato il problema�

y� � f�x� yy�x� � y�

� ����

supposta y�x �su�cientemente regolare� si ha

y�x� � h � y�x� � y�x� � hy��x� �h�

�y���x� � ��� �

hk

k�y�k �� ����

� � �x�� x� � h � x� � x��

Arrestandoci� nella ����� alla derivata prima� tenuto conto della ����� siha

y�x� � y�x� � hy��x� � y�x� � hf�x�� y� ����

che� dopo n passi� fornisce

yn�� � yn � hf�xn� yn� ���

La ��� prende il nome di metodo di Eulero� in questo caso ��xn� yn� h �f�xn� yn�

Page 123: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ���

Dalle ���� e ��� con n � �� �e immediato ottenere l�errore di troncamentoche si commette ad ogni passo� precisamente si ha

Ee � y�x�� y� �h�

�y�����

Osservazione ���� Il metodo di Eulero gode di una chiara interpretazionegeometrica cfr� la �gura ���� il valore assunto dalla soluzione nel punto x�y�x� �e approssimato con il valore y� assunto� nello stesso punto� dallaretta tangente la linea integrale y�x nel punto x��

−1 −0.5 0 0.5 1 1.5 2−1

0

1

2

3

4

5

6

7

8

X0 X1

y(X1)

y1

Errore

Eulero

y(x)

0

^

>

Figura ����

Il metodo di Eulero �e esatto �errore di troncamento nullo se la soluzioney�x del problema ���� �e una retta� diremo quindi che il metodo di Eulero �eun metodo del primo ordine�

Volendo ottenere un metodo d�ordine superiore �maggiore accuratezzasar�a necessario introdurre un altro termine dello sviluppo ���� contenentey���x�� Dal problema ���� �e facile ricavare tale termine osservando che

y���x �d

dx�y��x �

d

dx�f�x� y�x �

� fx�x� y�x � fy�x� y�xy��x �

� fx�x� y�x � fy�x� y�xf�x� y�x�

Page 124: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

Si ottiene cos�i un metodo del secondo ordine ponendo

yn�� � yn � hf�xn� yn �h�

��fx�xn� yn � fy�xn� ynf�xn� yn

che comporta� veri�carlo per esercizio� un errore di troncamento

E� � y�x�� y� �h�

�y������

In generale un metodo di Taylor di ordine k si ottiene ponendo

��xn� yn� h � Tk�xn� yn� h

dove

Tk�xn� yn� h � f�xn� yn �h

�y���xn � ��� �

hk��

k�y�k �xn

e le derivate y�j �x �j � �� �� ���� k sono calcolate derivando successivamentey��x � f�x� y� Il metodo diviene quindi

yn�� � yn � hTk�xn� yn� h

con errore di troncamento

Ek � y�x�� y� �hk��

�k � ��y�k�� ���

Osservazione ���� Nei metodi di Taylor Eulero escluso si ha l�incon�veniente che devono essere calcolate preventivamente le derivate y�j �x de�rivando �simbolicamente� l�equazione y��x � f�x� y� Questo pu�o portare�per k elevati� ad espressioni di Tk�xn� yn� h molto complesse anche quandof�x� y �e semplice� rendendo il metodo poco e�ciente bench�e accurato�

Prima di procedere nella determinazione di altri metodi �e opportuno sof�fermarci sulla seguente domanda� cosa vuole dire risolvere un problema diCauchy� Se la variabile indipendente x rappresenta il tempo possiamo direche risolvere un problema di Cauchy vuole dire� essere in grado di prevederel�evolvere nel tempo di un fenomeno �rappresentato dalla y�x che all�istan�te a � x� si trovava nella con�gurazione y�� Sar�a quindi cruciale riuscire adottenere� con il metodo numerico scelto� una buona stima della soluzione delproblema� non solo vicino� ma anche �lontano� dalla condizione iniziale�

Page 125: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ��

Nella determinazione dei precedenti metodi abbiamo supposto d�eseguireun solo passo d�integrazione �dalla condizione iniziale in x� al punto x� �x� �h� Nella pratica essendo interessati a conoscere la soluzione in tutto unintervallo �a� b �dalla condizione iniziale in a � x� al punto b � xn � x��nhdovremo eseguire� �ssato h� n passi d�integrazione� E� chiaro che� essendo�a� b �ssato� se si riduce l�ampiezza del passo h si dovr�a aumentare il numeron di passi in modo che nh � b � a� Questa banale constatazione �e crucialeper far so�ermare il lettore sul fatto che� contrariamente a quanto avvieneper i metodi iterativi per il calcolo degli zeri �metodi di punto �sso� l�indiced�iterazione nei metodi per equazioni di�erenziali fa spostare da un �istante��n al �successivo� �n�� e non �convergere� verso un punto �sso�

Per quanto sopra detto introduciamo le seguenti de�nizioni

De�nizione ���� Errore di troncamento globale� �e la di�erenza frala soluzione y�xn�� soluzione �vera� del problema nel punto xn�� e yn��

soluzione �calcolata� dopo n � � passi� Con ovvio signi�cato dei simboliabbiamo

en�� � y�xn��� yn���

De�nizione ���� Errore di troncamento locale� �e la di�erenza fra lasoluzione y�xn�� soluzione �vera� del problema nel punto xn�� e un�� so�luzione �calcolata� partendo dalla condizione y�xn eseguendo un solo passod�integrazione� Per cui

�n�� � y�xn��� un�� ����

doveun�� � y�xn � h��xn� y�xn� h

�e la soluzione del problema di Cauchy�u� � f�x� uu�xn � y�xn

calcolata numericamente nel punto xn � h�

De�nizione ���� Errore di troncamento locale unitario� �e l�errore�hn�� di troncamento locale rapportato al passo d�integrazione h�

�hn�� ��n��

h�

Page 126: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

De�nizione ��� Ordine� un metodo di integrazione numerica si dice diordine p se l�errore di troncamento locale unitario va a zero come hp

�hn�� � O�hp�

De�nizione ���� Consistente� un metodo di integrazione numerica �econsistente se �e del primo ordine o� equivalentemente� se

limh�

Tk�xn� yn� h � f�xn� yn� ����

Osservazione ���� E� facile veri�care l�equivalenza fra la �� e il fattoche

�hn�� � O�h

dalla ��� si ha

�hn�� �y�xn��� y�xn

h� ��xn� y�xn� h ����

e passando al limite nella ���� essendo

limh�

�hn�� � �

si ha

limh�

y�xn��� y�xn

h� y��xn � f�xn� yn � ��xn� y�xn� ���

De�nizione ��� Convergenza� un metodo di integrazione numerica �econvergente se

limh�nh�x

yn � y�x

riducendo il passo h� tenendo �sso il punto d�arrivo x� la soluzione numericastima sempre meglio la soluzione analitica�

Una chiara interpretazione delle precedenti de�nizioni �e data in �gura����

Per quanto detto la sola consistenza non ci basta per garantire la bont�adei risultati ottenuti con il metodo numerico ma �e necessario garantire laconvergenza� D�altro canto mentre �e� in generale� relativamente sempliceveri�care la consistenza �basta considerare un passo d�integrazione assai pi�ucomplesso �e veri�care la convergenza �n�� passi d�integrazione���

Dimostriamo il seguente

Page 127: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ���

x

y(x)

y(x )

u(x )

y n+1

Consistenza

n+1

n+1

0

y(x)u(x)

y(x )-y = y(x )-u(x )+u(x )-yn+1 n+1 n+1 n+1 n+1 n+1

Errore Stabilita’

Figura ����

Teorema ���� Sotto le ipotesi del teorema d�esistenza ed unicit�a il metododi Eulero �e convergente�

Dimostrazione ���� Dalla ��� si ha

yn�� � yn � hf�xn� yn

dallo sviluppo in serie di Taylor segue

y�xn�� � y�xn � hf�xn� y�xn �h�

�y�����

Per l�errore al passo �n � ��esimo si avr�a

en�� � y�xn��� yn�� � en � h �f�xn� y�xn� f�xn� yn �h�

�y����

per le condizioni di esistenza ed unicit�a

jf�xn� y�xn� f�xn� ynj � L jy�xn� ynjper cui

en�� � en � hL jy�xn� ynj� h�

�y����

ovvero

jen��j � �� � hL jenj� h�

�jy����j �

Page 128: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

Posto n � �� �� ���� �n � �� si ha

je�j � �� � hL je�j� h�

�jy����j ���

je�j � �� � hL� je�j� �� � hLh�

�jy����j� h�

�jy����j

� � �jen��j � �� � hLn�� je�j�

nXj��

�� � hLjh�

�jy����j �

Supposto e� � � errore su condizione iniziale nullo� moltiplicando e divi�dendo l�ultima delle ��� per ��� � hL� � � si ottiene

jen��j � �� � hLn�� � �

hL

h�

�jy����j

ricordando che ehL � � � hL si ricordi lo sviluppo di Taylor ��� si ha

jen��j � e�n�� hL � �

hL

h�

�jy����j � e�n�� hL � �

�Ljy����j h � Ch �����

per cui l�errore globale va a zero come il passo d�integrazione h convergenza��

Il risultato appena dimostrato per il metodo di Eulero pu�o essere genera�lizzato a tutti i metodi ad un passo in virt�u del seguente

Teorema ���� Dato un metodo numerico di tipo Taylor ad un passo sela Tk�x� y� h �e �lipschitziana� ed il metodo �e consistente allora il metodo �econvergente��

Osservazione ��� La consistenza� se la Tk�x� y� h �e �lipschitziana�� �equindi condizione necessaria e su�ciente per la convergenza�

Osservazione ���� Essendo� nel metodo di Eulero ed in tutti i metodi adun passo di tipo Taylor

limh�

Tk�x� y� h � f�x� y

la �lipschitzianeit�a� �e garantita dal teorema di esistenza ed unicit�a� per cuisi ha garantita� in virt�u del teorema precedente� la convergenza�

Page 129: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ��

Osservazione ��� Nella ���� l�errore va a zero come h� si �e persa unapotenza di h rispetto all�errore locale di discretizzazione� Questo fatto nondeve stupire in quanto� esattamente come per le formule di quadratura com�poste� l�aver diviso l�intervallo d�integrazione in n parti nh � xn � x�� ciobbliga ad eseguire n passi d�integrazione� Questa osservazione giusti�ca an�che l�introduzione della de�nizione di errore di troncamento locale unitariola divisione per h ci fa perdere una potenza nell�ordine di convergenza�

Nella deduzione della ����� abbiamo trascurato gli errori d�arrotonda�mento che si commettono ad ogni singolo passo� Volendo tener conto anchedi questi si avrebbe� detto � l�errore d�arrotondamento�

jen��j � e�n�� hL � �

hL

�h�

�jy����j� �

������

��

h

�Ljy����j� �

hL

��e�n�� hL � �

�per cui� mentre l�errore di discretizzazione va a zero con h� l�errore d�arro�tondamento cresce come h��� Nasce quindi� come nel caso della derivazionenumerica� la necessit�a di scegliere� �ssata la precisione di macchina �� un hottimo� Derivando rispetto ad h il secondo membro della ������ ignorandole costanti� ed uguagliando a zero� si ha

d

dh

�h

�jy����j� �

h

��

�jy����j� �

h�� �

per cui� indicato con y��M � max�jy����j � si ha

hott �s

��

y��M� c�

� � �����

La ����� ci fornisce un legame fra h ed � che pu�o essere cos�i parafrasato�pi�u �e piccolo � �maggiore �e l�accuratezza nei calcoli pi�u si pu�o prendere hpiccolo �indipendentemente dalla costante c� purtroppo l�esponente �

�con�

trolla male la riduzione del passo in quanto� per avere h dell�ordine del ����� �enecessario prendere � dell�ordine del �����

Page 130: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

����� Metodi Runge�Kutta

I metodi di Rung�Kutta di ordine p sono metodi ad un passo in cui ��x� y� hviene determinata in modo d�ottenere� senza e�ettuare il calcolo delle derivated�ordine superiore della y�x� un ordine d�accuratezza p� Per fare questo siusa una combinazione di r �metodi di Runge�Kutta a r stadi derivate primedi y�x in opportuni punti dell�intervallo �x� x � h �

Osservazione ���� Il metodo d�Eulero �e un metodo di Runge�Kutta delprimo ordine ad uno stadio si usa solo la y��xn � f�xn� yn�

Vediamo come �e possibile costruire un metodo del secondo ordine a duestadi� Scriviamo ��xn� yn� h nella forma

��xn� yn� h � a�K� � a�K�

dove a� e a� sono opportuni coe�cienti e

K� � f�xn� yn

K� � f�xn � �h� yn � �hK�

con � � � � �� sono approssimazioni dei valori di y��x nei punti xn e xn��hottenute mediante f�xn� yn e f�xn � �h� yn � �hK� rispettivamente� Ilmetodo diviene quindi

yn�� � yn � h �a�f�xn� yn � a�f�xn � �h� yn � �hK� �

Una chiara interpretazione geometrica �e fornita in �gura �����Per determinare i parametri a�� a� e � imponiamo la massima accuratezza

possibile per il metodo� il che equivale ad eguagliare� nello sviluppo in seriedi Taylor di y�xn e di ��xn� yn� h in un intorno di xn� pi�u termini possibili�Arrestando lo sviluppo ai termini del secondo ordine si ha� per la ��xn� yn� h

��xn� yn� h � a�f�xn� yn � a�f�xn � �h� yn � �hK� �����

� a�f�xn� yn � a��f�xn� yn � �hfx�xn� yn �

��hf�xn� ynfy�xn� yn � O�h�

� �a� � a�f�xn� yn � �ha��fx�xn� yn �

�f�xn� ynfy�xn� yn � O�h��

Page 131: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ���

x

y(x)

0

y(x)u(x)

x x+hx+ μ h

k

y(x+h)

2=f(x+ μ h,y+ μ hk1

)

k 1 =f(x,y)

a1 k 1 +a2 k 2

y

yh

Figura ����

mentre per y�xn si ha

y�xn � h� y�xn

h� f�xn� yn �

h

��fx�xn� yn � f�xn� ynfy�xn� yn �O�h��

�����

Eguagliando i coe�cienti dei termini corrispondenti nei due sviluppi siottiene il seguente sistema non lineare�

a� � a� � ��a� � �

����

che� fra le in�nite soluzioni� ammette anche

�� metodo di Eulero modi�cato �a� � �a� � �� � �

che �e della forma

yn�� � yn � hf

�xn �

h

�� yn �

h

�f�xn� yn

��

Page 132: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

�� metodo di Heun �a� � �

a� � ��

� � �

che �e della forma

yn�� � yn �h

��f�xn� yn � f �xn � h� yn � hf�xn� yn �

Dalle ����� e ����� �e immediato osservare che i metodi risultano delsecondo ordine essendo l�errore di troncamento locale unitario O�h��

In �gura ���� �e riportata l�interpretazione geometrica del metodo diHeun�

x

y(x)

0 x

y n+1y(x)

f(x ,y )n n

n x n +h

f(x +h,y +hf(x +y ))n n n n

media

Figura ����

Osservazione ���� Nella ���� la condizione a� � a� � � implica la con�sistenza del metodo� per cui� sotto le consuete ipotesi di esistenza ed unicit�a�tutti i metodi ricavati dalla ���� risultano convergenti�

Quanto visto per due stadi pu�o essere esteso �non con facili conti� alcaso generale di r stadi� per cui otteniamo

��xn� yn� h �rX

i��

aiKi�x� y� h

Page 133: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ���

dove

K� � f�xn� yn �����

Ki � f�xn � �ih� yn � hi��Xj��

�ijKj� i � �� �� ���� r

che prendono il nome di metodi Runge�Kutta a r stadi espliciti� dove l�agget�tivo esplicito evidenzia il fatto che ogni Ki dipende solo dai Kj precedenti�Nella ����� i parametri ai� �i e �ij verranno scelti in modo d�ottenere unmetodo il pi�u accurato possibile�

Per i metodi espliciti� detto p l�ordine del metodo ed r il numero di stadi�vale la seguente regola� �al crescere del numero degli stadi non si guadagnasempre parimenti in accuratezza�� si ha infatti che

� esistono metodi con p � r solo per r � �� �� �� �

� per r � � �� � si hanno solo metodi di ordine p � r � �

� per r � �� si hanno solo metodi di ordine p � r � �

� per r � �� l�ordine �e p � r � ��

Per concludere ricordiamo il classico metodo di Runge�Kutta del quartoordine a quattro stadi

��xn� yn� h ��

��K� � �K� � �K� � K�

dove

K� � f�xn� yn

K� � f�xn �h

�� yn �

h

�K�

K� � f�xn �h

�� yn �

h

�K�

K� � f�xn � h� yn � hK��

Se nella ����� l�indice j della sommatoria pu�o andare �no ad r �ogniKi dipende da tutti gli altri si hanno i metodi di Runge�Kutta impliciti�se l�indice j della sommatoria pu�o andare �no ad i �ogni Ki dipende daiprecedenti Kj e da se stesso si hanno i metodi di Runge�Kutta semi�impliciti�per i quali si rimanda alla letteratura speci�ca�

Page 134: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

����� Sistemi del primo ordine

Quanto visto per un�equazione del primo ordine pu�o essere facilmente estesosia a sistemi del primo ordine che ad equazioni di grado superiore al primo�come di seguito illustrato per un�equazione del secondo ordine�

Dato il seguente problema del secondo ordine�y�� � f�x� y� y�y�a � y��a � �

abbiamo visto che �e possibile trasformarlo nel sistema equivalente�y� � zz� � f�x� y� z

con

�y�a � z�a � �

ed i metodi precedenti possono essere generalizzati al nuovo sistema utiliz�zando� con ovvio signi�cato dei simboli� la seguente notazione vettoriale

u� � f�x� u

ove

u �

�yz

�� u� �

�y�

z�

�� f �

�z

f�x� y� z

�con le condizioni

u�a � �� ove � �

��

��

Per cui� dato il metodo numerico

yn�� � yn � h��xn� yn� h

si potr�a scrivere� per il sistema precedente

un�� � un � h��xn� un� h�

Esempli�chiamo quanto sopra nel caso del metodo di Heun

yn�� � yn �h

��f�xnyn � f �xn��� yn � hf�xn� yn �

� yn �h

��K� � K�

Page 135: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI AD UN PASSO �M�F�� ��

ed all�equazioney��� � x�y�� � y�y��

posto �y � u�y� � u�y�� � u�

si ha� in forma vettoriale� �u�� � u�u�� � u�u�� � x�u� � u��u�

Un semplice programma in Matlab� per risolvere il problema precedente�pu�o essere il seguente

function uprime�uprim�x u�

uprime��u����u����u���x���u���u�������

return

end

������ definizioni iniziali

h����

x����

uold���

while x � xfin

k��uprim�x uold��

x�x�h�

utemp�uold�hk��

k��uprim�x utemp��

unew�uold�h�k��k�����

uold�unew�

����continua fino a che x � xfin

end�����fine while

Come si pu�o osservare dall�attenta lettura del listato del programma� pertrattare� con il metodo di Heun� altri sistemi di�erenziali basta modi�careopportunamente solo la function uprim�

Page 136: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

����� Stabilit�a numerica

Per i metodi precedenti abbiamo garantito un certo ordine di accuratezza ela convergenza per h � �� Quando si deve integrare numericamente� in unintervallo pre�ssato �a� b � un�equazione di�erenziale� il passo h deve essere�ovviamente� preso diverso da zero� Nasce spontanea la domanda� per ognipasso h la soluzione numerica sar�a una buona stima di quella reale � O equi�valentemente� gli errori introdotti in ogni singolo passo d�integrazione nonfanno �esplodere� la soluzione numerica � Questi interrogativi permettonodi introdurre il concetto di stabilit�a numerica� Diamo le seguenti

De�nizione ���� La soluzione di un problema di Cauchy �e detta inerente�mente stabile se

jy�xn��j � C jy�xnj � C � �� �n�De�nizione ���� Un metodo numerico si dice assolutamente stabile se

jyn��j � c jynj � c � �� �n�dove yn �e la soluzione numerica nel punto xn�

Per studiare la stabilit�a numerica considereremo �per semplicit�a d�espo�sizione l�equazione test �

y� � �yy�� � y�

� Re���� �����

che ammette come soluzione

y�x � y�e�x �����

che� essendo Re�� � �� �e inerentemente stabile�Se consideriamo il metodo di Eulero� otteniamo la soluzione

yn�� � yn � h�yn � �� � h� yn

e quindiyn�� � �� � h�n�� y�� ����

La soluzione ���� si comporta� per n��� come la soluzione ����� see solo se

j� � h�j � �� �����

Page 137: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI A PI �U PASSI �M�F�� ���

La ����� rappresenta la regione di stabilit�a assoluta del metodo di Eulero�se non stiamo all�interno di questa regione �cerchio del piano complesso concentro in ���� � e raggio unitario la soluzione numerica crescer�a in moduloallontanandosi sempre pi�u dalla soluzione vera y�x �che tende a zero pern���

Analogamente si potrebbe calcolare la regione di stabilit�a assoluta per imetodi di Heun e di Runge�Kutta del quarto ordine� per i quali si rimandaalla bibliogra�a speci�ca�

Per tutti i metodi ad un passo visti esiste quindi una condizione sulpasso d�integrazione h che� se violata� produce una soluzione inattendibi�le� Questo fatto si traduce dicendo che i metodi Runge�Kutta e di Taylorsono condizionatamente stabili�

�� Metodi a pi�u passi �M�F��

I metodi a p passi sono metodi che approssimano la soluzione in un puntoutilizzando i valori trovati in p passi precedenti� Si possono scrivere nellaforma

yn�p � �p��Xj��

jyn�j � h

pXj��

�jy�n�j � �����

� �p��Xj��

jyn�j � h

pXj��

�jf�xn�j� yn�j

dove i �p � � parametri j e �j individuano il singolo metodo� Se �p � �il metodo si dir�a esplicito� mentre se �p �� � il metodo sar�a detto implicito�Per determinare i parametri j e �j si possono seguire di�erenti approcci

�� utilizzando le formule di derivazione numerica� l�equazione di�erenzialefornisce un legame fra la funzione y�x e la sua derivata prima�

�� utilizzando le formule di quadratura� il problema ���� �e equivalente alproblema Z x��h

x�

y��xdx �

Z x��h

x�

f�x� ydx

ovvero

y�x� � h � y�x� �

Z x��h

x�

f�x� ydx� �����

Page 138: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

in cui compare il calcolo di un integrale�

�� utilizzando il concetto di ordine d�accuratezza� imponendo che la for�mula sia esatta ogniqualvolta la soluzione della ���� �e un polinomio digrado al pi�u �p�

Esempio ���� dalle formula di derivazione� consideriamo la formula

y��x �y�x � h� y�x� h

�h� O�h� �����

si ottiene il metodoyn�� � yn � �hy�n�� �����

per cui � � ��� � � �� �� � �� �� � �� �� � � che �e un LMM Linear Multystep Method a due passi� Il passaggio dalla ���� alla ���� si giusti�cafacilmente considerando� dovendo integrare la ��� nell�intervallo �x�h� x�h � il seguente schema�

yn�

y�x�k ���� yn��

�y�x

���� yn���

y�x�k ����

Esempio ���� dalle formule di quadratura� consideriamo la formula delpunto medio Z x��h

x��hf�xdx � �hf�x� �

h�

�f ����

applichiamola� tenuto conto del cambio d�intervallo e dello schema �����alla ���� ed otteniamo ancora

yn�� � yn � �hy�n��

che �e lo stesso metodo ottenuto nell�esempio precedente�

Esempio ���� imponendo l�ordine d�accuratezza� imponiamo che il metodonumerico a � passi

yn�� � � ��yn � �yn�� � h���y

�n � ��y

�n�� � ��y

�n��

������

sia esatto errore di troncamento nullo quando la soluzione �e un polinomiodi grado minore od uguale a �� Ricordando che i monomi �� x ed x� veri�canorispettivamente i problemi di Cauchy�

y� � �y��h � �

�y� � �y��h � �h �

�y� � �xy��h � h�

Page 139: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI A PI �U PASSI �M�F�� ��

le soluzioni ottenute con il metodo ���� nell�intervallo ��h� h si confrontilo schema seguente

�h�n���� �

�n��

���� h�n��

portano al seguente sistema lineare�� � � �� � � � h � �h � � ����h � � � � � h ��� � �� � ��h� � � ��h

� � � � � � h ���h�� � �� � � � �h���

Imponendo al metodo di essere esplicito �� � � si ha�� � �� � �

� � � � �� � ��

� � �� � ���

avendo ancora un grado di libert�a � equazioni e � parametri possiamo porre� � �� ottenendo cos�i � � ��� � � �� �� � �� �� � �� �� � � che �e ancorail metodo �����

Dall�ultimo esempio considerato si ottiene che il metodo del punto medio

yn�� � yn � �hy�n��

�e un metodo LMM a � passi esplicito e del secondo ordine�

Osservazione ���� Che l�errore di troncamento locale unitario sia O�h�lo si evince dalla ����� che poi il metodo non sia esatto per polinomi digrado � �e facilmente veri�cabile sul problema di Cauchy�

y� � �x�

y��h � �h� �

I metodi della forma ����� devono appoggiarsi� nei primi p� � passi� adun metodo di integrazione ad un passo �Taylor o Runge�Kutta che forniscale prime p� � valutazioni necessarie ad innescare il metodo� Se il metodo �eesplicito proceder�a poi automaticamente nelle iterazioni sucessive mentre� se�e implicito� bisogner�a risolvere ad ogni passo un�equazione implicita� A frontedi questi due inconvegnenti si ha il vantaggio che ad ogni singolo passo si ha

Page 140: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

un metodo di ordine p � p � � se implicito con una sola nuova valutazionedella f�xn� yn�

L�equazione implicita da risolvere ad ogni passo si presenta nella forma�cfr� �����

yn�p � h�pf�xn�p� yn�p � gn �����

dove

gn � �p��Xj��

jyn�j � h

p��Xj��

�jf�xn�j� yn�j

non dipende da yn�p� Per risolvere la ����� si pu�o utilizzare un metodo dipunto �sso nella forma

y�i�� n�p � h�pf�xn�p� y

�i n�p � gn� i � �� �� ��� �����

partendo da un punto iniziale y�� n�p �per esempio y

�� n�p � yn�p��� Ricordando

la condizione di convergenza per i metodi del primo ordine �cfr� teorema����� si ha il seguente

Teorema ���� Se f�x� y �e lipschitziana

kf�x� y�� f�x� yk � L ky� � yk � �a � x � b� y� y� � R

ed inoltreh���p��L � � ����

allora la successione ���� converge��

Osservazione ���� La condizione ���� a�erma che� pur di prendere hsu�cientemente piccolo� l�equazione implicita ��� pu�o essere risolta con ilmetodo di punto �sso �����

Come per i metodi ad un passo possiamo dare le de�nizioni di consistenzaed ordine per i LMM�

De�nizione ���� Un LMM della forma ���� �e consistente se l�errrore ditroncamento locale unitario va a zero con h� ovvero�

limh�

�hn � limh�

�y�xn � h �

Pp��j�� jy�xn � jh

h�

pXj��

�jy��xn � jh

� ��

Page 141: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI A PI �U PASSI �M�F�� ���

De�nizione ���� Un LMM della forma ���� �e di ordine r se

�hn � O�hr�

Osservazione ���� E� interessante osservare che mentre un LMM di or�dine r integra esattamente a meno degli errori di arrotondamento un qua�lunque problema di Cauchy la cui soluzione �e un polinomio di grado minoreod ugual ad r� questo non �e pi�u vero� in generale� se si usa un metodo diRunge�Kutta di ordine r� Si veda l�esercizio ������

Per quanto riguarda la convergenza nei metodi LMM la consistenza �datala lischipzianit�a di f�x� y �e solo condizione necessaria ma non pi�u su�cien�te� Infatti� mentre nei metodi ad un passo l�equazione di�erenziale del primoordine viene sostituita con un� equazione alle di�erenze anch� essa del primoordine� la cui unica soluzione �e un� approssimazione della soluzione del pro�blema di�erenziale� nei LMM l�equazione alle di�erenze� di ordine p � �� pu�oammettere �no a p soluzioni distinte� di cui una sola �soluzione principaleapprossima correttamente la soluzione del problema di�erenziale� mentre lerimanenti �soluzioni parassite introducono solo dell�errore� E� chiaro quin�di che� per avere convergenza� al tendere di h a zero� le soluzioni parassitedovranno essere inin�uenti�

Ricordando che data l�equazione alle di�erenze �a coe�cienti costanti

yn�p � �p��Xj��

jyn�j �����

�ottenuta� guarda caso� dalla ����� per h � � ammette la soluzione generale

ym �

p��Xj��

�j�rjm

dove le rj sono le p radici distinte del polinomio di grado p� associato alla����� �polinomio caratteristico�

�p �

p��Xj��

j�j � �

e le �j sono p costanti arbitrarie �determinabili una volta dati i primi p valoriy�� y�� ���� yp��� possiamo dare la seguente

Page 142: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

De�nizione ���� Un LMM della forma ���� �e zero stabile se

jrjj � �� �j �����

e se jrjj � � allora rj �e semplice� La ���� �e detta anche root condition ocondizione delle radici�

Si pu�o dimostrare il seguente

Teorema ���� Un metodo LMM �e convergente se e solo se �e consistente ezero stabile��

Dipendendo le radici del polinomio caratteristico dai coe�cienti j que�sti potranno essere scelti in modo opportuno� non solo per garantire laconsistenza e l�ordine ma anche la zero stabilit�a�

����� Metodi di Adams

I metodi di Adams�Bashforth e Adams�Moulton sono della forma

yn�p � yn�p�� � h

p��Xj��

�jy�n�j� �Adams� Bashforth

yn�p � yn�p�� � h

pXj��

�jy�n�j� �Adams�Moulton

dove i coe�cienti �j dipendono dalla particolare formula di quadratura uti�lizzata per approssimare� nell�identit�a

yn�p � yn�p�� �

Z xn�p

xn�p��

y��xdx

l�integrale a secondo membro� Se nella formula di quadratura si usa ilnodo xn�p �si �e costruito il polinomio interpolante su tutti i p � � nodixn� xn��� ���� xn�p si hanno i metodi di Adams�Moulton che sono impliciti�se invece nella formula di quadratura non si usa il nodo xn�p �si �e costruito ilpolinomio interpolante solo sui p nodi xn� xn��� ���� xn�p�� si hanno i metodidi Adams�Bashforth che sono espliciti�

Per l�ordine di accuratezza si pu�o dimostrare che� a parit�a di passi p� gliAdams�Bashforth hanno ordine p� mentre gli Adams�Moulton hanno ordinep � ��

Page 143: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI A PI �U PASSI �M�F�� ���

Questi metodi fanno tutti parte della famiglia dei metodi LMM con lacaratteristica di essere sicuramente zero stabili in quanto il loro polinomiocaratteristico� essendo p�� � �� e j � � per j � �� �� ���� p� �� �e

�p � �p�� � �

e quindi ha radici tutte nulle tranne una uguale ad � �root condition veri��cata�

Fra i metodi di Adams ricordiamo�

�� Adams�Bashforth

�a yn�� � yn � hf�xn� yn� metodo di Eulero �p � � � ordine ��

�b yn�� � yn�� � h�

��f�xn��� yn��� f�xn� yn � �p � � � ordine ��

�� Adams�Moulton

�a yn�� � yn� h�

�f�xn��� yn�� � f�xn� yn � metodo dei Trapezi �p �� � ordine �

�b yn�� � yn���h��

� f�xn��� yn�� � �f�xn��� yn��� f�xn� yn � �p �� � ordine �

����� Stabilit�a numerica

Per i metodi LMM oltre che di stabilit�a assoluta �come per i metodi ad unpasso si pu�o parlare anche di stabilit�a relativa� in quanto� anche quandola soluzione �e inerentemente instabile non �e garantito che �h la soluzionenumerica� che dipende da tutte le radici del polinomio caratteristico� �segua�bene la soluzione analitica� Analizzeremo brevemente questi concetti su dueclassici metodi� il metodo dei trapezi ed il metodo del punto medio�

Stabilit�a assoluta

Come per i metodi ad un passo ci proponiamo di studiare� �ssato il passo hd�integrazione� come si comporta la soluzione numerica calcolata con i LMM�Ci limiteremo� per motivi di tempo� ad analizzare la formula dei Trapezi�per una trattazione di carattere generale si rimanda ai testi in bibliogra�a�

Page 144: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

Consideriamo la solita equazione test ������ il metodo dei Trapezi ad essaapplicato produce

yn�� � yn �h

���yn�� � �yn

che� risolta rispetto a yn��� fornisce

yn�� �� � h�

�� h��

yn

ed� essendo Re�� � �� si ha ������ � h��

�� h��

����� � �� �����

per cui la soluzione numerica �e assolutamente stabile per ogni scelta del passoh�

Osservazione ��� La regione di stabilit�a del metodo dei Trapezi �e quinditutto il semipiano negativo del piano complesso� Questa maggiore �libert�a�nella scelta del passo h �e stata pagata con la necessit�a di dover risolvere� adogni passo� un�equazione implicita� In generale i metodi impliciti sono �pi�ustabili� degli espliciti�

Nella ����� la quantit�a��h�

��h��

�e l�unica radice del polinomio associato al

metodo dei Trapezi applicato all�equazione test ������ in generale� dette �ile p radici �supposte distinte del polinomio associato al generico LMM ap�plicato all�equazione test ������ la condizione di assoluta stabilit�a si traducein

j�ij � �� i � �� �� ���� p�

Stabilit�a relativa

Se nel problema test ����� abbiamo Re�� � � la soluzione tender�a� pern��� ad in�nito� Per i LMM si pu�o quindi introdurre il concetto di stabilit�arelativa� nel senso che� �e ragionevole richiedere che la soluzione numerica �chedovr�a anch�essa crescere si mantenga �vicino� alla soluzione vera �errorerelativo limitato� Poich�e la soluzione generale dell�equazione alle di�erenzedipende da tutte le p radici �i dovremo avere

j�ij � j��j � i � �� �� ���� p

Page 145: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���� METODI A PI �U PASSI �M�F�� ��

dove con �� si �e indicata la radice principale e con �i� i � �� �� ���� p leradici parassite le radici parassite devono essere controllate dalla radiceprincipale�

Vediamo di chiarire le ultime a�ermazioni fatte studiando la stabilit�a delmetodo del punto medio ������ Dal metodo

yn�� � yn � �hy�n��

e dal problema test ����� si ottiene l�equazione alle di�erenze

yn�� � �h�yn�� � yn � �

il cui polinomio caratteristico �e

p�r � r� � �h�r � � �����

se � � � allora si hanno due radici distinte� � � r� � � e r� � ��� come �efacile veri�care osservando che

p��� � �p�� � �h� � �� p�� � ���

La soluzione generale dell�equazione alle di�erenze

yn � Arn� � Brn�

�dove A e B sono costanti da determinarsi� si comporta� per n��� comeyn � Arn� in quanto il termine in r� tende a zero� Essendo r� � �� la soluzio�ne cresce in modulo per cui il metodo non �e assolutamente stabile� Se venisseutilizzato per risolvere un problema di Cauchy con soluzione inerentementestabile produrrebbe� al crescere di n� una soluzione numerica inattendibile�

Se nel problema test � � � �soluzione inerentemente instabile le radicidel polinomio ����� divengono �� � r� � � e r� � �� come �e facile veri�careosservando che

p��� � �p�� � �h� � �� p�� � ���

La soluzione generale dell�equazione alle di�erenze si comporta� in questocaso come yn � Brn� � poich�e� per n � �� il termine in r� tende a zero�Il metodo �e quindi relativamente stabile e fornir�a una soluzione numericacon errore relativo limitato �la radice parassita r� al crescere di n in�uenzasempre meno la soluzione�

Page 146: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

Osservazione ���� Garantita la stabilit�a relativa bisogner�a scegliere il pas�so d�integrazione h su�cientemente piccolo per garantire l�accuratezza�

Osservazione ��� Nei metodi ad un passo si parla solo di stabilit�a asso�luta in quanto non ha senso parlare di stabilit�a relativa si osservi al riguardoche il polinomio caratteristico ha una sola radice�

� Riepilogo �M�F��

Sintetizzando quanto esposto nel presente capitolo abbiamo

�� I metodi ad un passo �Taylor e Runge�Kutta� possono essere uti�lizzati sfruttando la sola condizione iniziale fornita nel problema diCauchy� mentre� per utilizzare i metodi LMM a p passi� �e necessariopreventivamente ottenere le prime p� � valutazioni�

�� Per aumentare l�accuratezza nei metodi di Taylor e Runge�Kutta �enecessario aumentare il numero di valutazioni funzionali ad ogni passo�mentre� per i LMM a p passi basta sempre una sola valutazione puravendo ordine p �espliciti o p � � �impliciti�

�� Il concetto di consistenza permette di legare l�errore di troncamentolocale alla convergenza e precisamente

�a la consistenza �e condizione necessaria e su�ciente per la conver�genza nei metodi a un passo�

�b la consistenza �condizione necessaria � la zero stabilit�a �rootcondition garantiscono la convergenza nei metodi a p passi�

�� Il concetto di zero stabilit�a si introduce solo per i metodi a p passi inquanto� con questi metodi� si sostituisce al problema continuo� dipen�dente da una sola condizione iniziale� un problema discreto dipenden�te da p condizioni �per �ssare le p costanti arbitrarie della soluzionegenerale dell�equazione alle di�erenze associata al metodo numerico�

� Per lo stesso motivo solo per i LMM si pu�o parlare di stabilit�a relativa�

Page 147: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

��� RIEPILOGO �M�F�� ���

�� Oltre all�accuratezza �e cruciale il concetto di assoluta stabilit�a che pu�ovincolare il passo h d�integrazione �regione di assoluta stabilit�a tantoda rendere praticamente inutilizzabile un metodo �se h �e troppo piccolonon si procede nell�integrazione�

�� I metodi impliciti� pur essendo pi�u laboriosi �equazione da risolvere adogni passo� hanno regioni di stabilit�a pi�u ampie �permettono passi hpi�u grandi�

����� Esercizi

Esercizio ���� Dato il problema di Cauchy�y� � y � x� � �x� �y�� � �

la cui soluzione �e y�x � x� � �� calcolare con passo h � ��� la soluzione in��� � con

�� il metodo di Heun�

�� il metodo dei Trapezi esplicitando yn�� in funzione di yn�

Giusti�care i risultati ottenuti�

Esercizio ���� Integrare nell�intervallo ��� � il seguente problema di Cau�chy �

y� � �� yy�� � �

come passi d�integrazione h � ���� e h � ��� utilizzando

�� la formula del punto medio�

�� la formula dei trapezi�

Giusti�care i risultati ottenuti�

Esercizio ���� Dato il metodo numerico

yn�� ��

�yn�� � �yn�� �

�yn � �hf�xn��� yn���

Page 148: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

���CAPITOLO �� EQUAZIONI DIFFERENZIALI ORDINARIE �M�FRONTINI�

�� �e consistente �

�� E� convergente�

Giusti�care le risposte date�

Esercizio ��� Scrivere una function in MATLAB richiamabile con

�xn� yn���adabas��xn xn� h yn yn� f�

che esegua un passo del metodo d�integrazione di Adams�Bashforth diordine ��

Esercizio ���� Integrare nell�intervallo ��� �� � utilizzando diversi passi d�in�tegrazione� il sistema �

y� � x�tx� � �y�t

�y�� � �x�� � �

mediante

�� il metodo di Eulero�

�� il metodo di Heun�

Dire qual��e la linea i cui punti hanno coordinate P �t � �x�t� y�t etracciarne il gra�co�

Page 149: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

Bibliogra a

�� K�E�Atkinson� An Introduction to Numerical Analysis� John Wiley #Sons� New York� ���

�� J�C�Butcher� The Numerical Analysis od Ordinary Di�erential Equa�tions� Runge�Kutta and General Linear Methods� John Wiley # Sons�Chichester� ���

�� S�D�Conte e C�de Boor� Elementary Numerical Analysis� McGraw�HillKogakusha� Ltd�� Tokio� ���

�� G�Dahlquist e A�Bjorck� Numerical Methods� Prentice�Hall� EnglewoodCli�s� NJ� ����

� F�Fontanella e A�Pasquali� Calcolo Numerico� Metodi e algoritmi�Vol����� Pitagora Editrice� Bologna� ���

�� G�Gambolati� Metodi Numerici� Cortina� Padova� ��

�� L�Gotusso� Calcolo Numerico� Clup� Milano� ���

�� G�H�Golub e C��Van Loan� Matrix Computation� The John HopkinsPress� Baltimore� ��

� J�D�Lambert� Numerical Methods for Ordinary Di�erential Systems�The Initial Value Problem� John Wiley # Sons� Chichester� ��

��� A�M�Ostrowski� Solution of equations and systems of equations�Academic Press� New York� ���

��� J�R�Rice� The Approximation of functions� Vol I�II� Addison�Wesley�New York� ��

��

Page 150: INDICE Esercizi - Altervista · Decomp osizione LU Decomp osizione di Cholesky Analisi dell errore MF Ric hiami sulle norme di v ettore. INDICE Stima dell errore Meto di iterativi

� � BIBLIOGRAFIA

��� J�Stoer e R�Bulish� Introduction to Numerical Analysis� Springer�Verlag�New York� ���

��� G�Strang� Linear algebra and its applications� Academic Press� NewYork� ���

��� A�Stroud� Numerical Quadrature and Solutio of Ordinary Di�erentialEquations� Springer�Verlag� New York� ���