Laboratorio Matematico-Informatico: Teoria · a.a 2010-2011, docente: Francesca Mazzia 1.1 Calcolo...

Capitolo 1

Laboratorio

Matematico-Informatico: Teoria

Questo capitolo riassume le lezioni di teoria del corso di Laboratorio Matematico Informaticoa.a 2010-2011, docente: Francesca Mazzia

1.1 Calcolo Scientifico

Il calcolo scientifico rappresenta l’insieme degli strumenti, delle tecniche e delle teorie ne-cessarie per risolvere con il computer problemi della scienza e della tecnica. Il suo scopoprincipale e sviluppare metodi efficienti per calcolare approssimazioni di quantita che sonodifficili o impossibili da ottenere per mezzi analitici.

Il calcolo scientifico nasce negli anni ’40 con la modellizzazione di problemi balistici e diarmi nucleari nella Seconda Guerra Mondiale; si sviluppa negli anni ’70 e ’80 per risolvereproblemi industriali, es. progettazione di aerei; si consolida negli anni ’90 su problemi diInformatica, Fisica, Chimica, Ingegneria, Economia, ....

Una buona parte del calcolo scientifico richiede la risoluzione di problemi matematici alcalcolatore che puo essere eseguita in modo efficiente solo se:

• si conosce lo strumento computazionale da usare.

• si conosce il problema da risolvere

• si costruisce un algoritmo che risolve il problema con una accuratezza desiderata edentro i limiti delle risorse (tempo, memoria, ...) disponibili.

L’ ambiente computazionale utilizzato per risolvere problemi matematici al calcolatoree costituito da:

• hardware: PC o supercomputer (computer vettoriali, computer paralleli);

• sistemi operativi e linguaggi: UNIX, LINUX, WINDOWS, ... C, C++, Java, For-tran95, ...

2 CAPITOLO 1. LABORATORIO MATEMATICO-INFORMATICO: TEORIA

• strumenti per il data management: costruzione di database contenenti tutte le infor-mazioni rilevanti per un particolare progetto applicativo;

• strumenti per la visualizzazione: rappresentazione grafica finalizzata ad una velocecomprensione del calcolo.

Il linguaggio di programmazione per la comunita scientifica e stato per molti anni ilFortran (introdotto negli anni ’50). Ora anche altri linguaggi, specialmente il C e il C++ eJava, sono usati per il Calcolo Scientifico. La disponibilita di un linguaggio che puo essereusato su una molteplicita di macchine determina la portabilita del software da macchina amacchina con pochissimi cambiamenti nelle performance del codice.

Ambienti che sono molto utilizzati per risolvere problemi del calcolo scientifico sono iProblem Solving Environments (PSE). Un PSE e un ambiente di alto livello che contienetutto il software necessario per risolvere una determinata classe di problemi, fra cui:

• metodi avanzati per la soluzione;

• selezione automatica o semiautomatica dei metodi di soluzione;

• modi per incorporare facilmente nuovi metodi di soluzione.

Gli utenti possono usarlo senza avere una conoscenza specialistica del computer. Alcunidei PSE piu usati per il calcolo scientifico sono: Matlab; Maple; Mathematica; Scilab;Octave; R; Python; Maxima.

In questo corso utilizzeremo il Matlab e il Maple. Il MATLAB (MATrix LABorato-ry) http://www.mathworks.com/ e un PSE per il calcolo scientifico commerciale ad alteprestazioni, integra analisi numerica, calcolo con matrici, grafica.

Il MAPLE http://www.maplesoft.com/ e un PSE per il calcolo scientifico commercia-le ad alte prestazioni. Fu sviluppato per la prima volta nel 1981 dal Symbolic Computa-tion Group all’Universita di Waterloo, Canada ed e particolarmente indicato per il calcolosimbolico.

Costruire un algoritmo efficiente per il calcolo scientifico puo essere una sfida complessae difficile. Lo scopo di questo corso e farvi conoscere metodi e strumenti di base da poter uti-lizzare per risolvere problemi piu complessi che potreste incontrare in futuro. In particolareci soffermeremo su alcuni problemi di particolare importanza quando si lavora con i numerial calcolatore fra cui: la rappresentazione dei numeri di macchina, il condizionamento di unproblema e la stabilita degli algoritmi.

1.2 Il processo di risoluzione numerica

Per risolvere un problema del mondo reale al calcolatore si eseguono i seguenti passi:

• Problema reale → Modello Matematico: questo passo spesso richiede una semplifica-zione del modello, per esempio si suppongono trascurabili alcune grandezze fisiche

1.3. SORGENTI DI ERRORI 3

• Modello Matematico → Metodo Numerico: questo passo spesso genera delle appros-simazione nel metodo risolutivo (es. un procedimento infinito approssimato medianteun procedimento finito: errore di troncamento);

• Metodo Numerico → programma: in questo passo si introducono gli errori di arro-tondamento: i dati numerici elaborati e i risultati delle operazioni eseguite vengonoarrotondati (si opera con aritmetica finita, errore di round-off)

Il processo di risoluzione numerica si occupa dello sviluppo di metodi numerici che ten-gano conto della natura del problema, delle risorse hardware e software, dello sviluppo di unalgoritmo o un set di istruzioni che descrivono come risolvere un problema con un calcolato-re. Questo richiede l’esame dei diversi aspetti che riguardano la scrittura, la comprensione,la valutazione di algoritmi e la loro implementazione in un determinato ambiente di calcolo;la formulazione di problemi test; lo sviluppo della fase di testing e della misura dell’effi-cienza. L’ efficienza corrisponde alla ottimizzazione della complessita di tempo e di spazioe richiede uno studio accurato dell’algoritmo e della fase di implementazione attraversola ricerca di tecniche che consentono l’effettiva minimizzazione dei calcoli e una correttapianificazione delle aree di memoria da utilizzare.

Alcuni criteri per un buon codice sono:

• affidabilita;

• robustezza;

• portabilita;

• leggibilita;

• buona documentazione;

• ampia fase di testing.

1.3 Sorgenti di errori

Nel processo di risoluzione numerica e necessario essere capaci di controllare le diversesorgenti di errori in modo che non interferiscano con i risultati calcolati. Di solito gli errorisi propagano dalla loro sorgente a quantita che sono calcolate dopo, alcune volte con unfattore di amplificazione considerevole. Fra le diverse sorgenti di errori abbiamo:

• Errori nei dati di input: sono errori sistematici che dipendono dallo strumen-to di misura utilizzato per rilevare i dati di input o errori casuali che dipendonodall’ambiente di sperimentazione.

• Errori di arrotondamento: π non e rappresentato esattamente nel nostro compu-ter. Il prodotto esatto di due numeri con s cifre contiene 2s o 2s− 1 cifre e per essererappresentato nel nostro calcolatore deve essere arrotondato. Gli errori di arroton-damento sono legati alla rappresentazione dei dati ed esecuzione delle operazioni inaritmetica finita.


• Errori di troncamento: questi errori si verificano quando un processo infinito vienetroncato, per esempio una serie infinita viene troncata dopo un numero finito i termini.

• Semplificazione del Modello: passando dal problema reale al modello matematico,per esempio, si suppongono trascurabili alcune grandezze fisiche.

• Errori Umani: In tutti i lavori di tipo numerico bisogna aspettarci errori umanicome: errori in calcoli manuali. Bisogna sempre considerare che anche i libri di testopossono contenere errori in tabelle e formule. Ci possono essere errori nei programmi,errori di battitura o scambio di operatori e spesso errori di macchina.

Bisogna essere abili nel verificare la correttezza dei risultati. Ci sono errori controlla-bili e incontrollabili, Gli errori nei dati di input e le semplificazioni del modello possonoessere considerati incontrollabili. Gli errori di troncamento sono controllabili. Gli errori diarrotondamento sono controllabili fino ad un certo punto.

Se A e una quantita che vogliamo calcolare e A e un’approssimazione di A, allora l’errorecommesso e la differenza fra i due valori:

errore = A− A;

L’errore assoluto e il valore assoluto dell’errore:

errore assoluto = |A− A|;e l’errore relativo si ottiene normalizzando l’errore assoluto con il valore esatto, se A 6= 0:

errore relativo =|A− A||A| .

L’errore relativo e piu significativo dell’errore assoluto. E ragionevole chiedere che l’er-rore relativo sia minore di un valore prefissato. Se conosciamo una maggiorazione dell’erroreassoluto, cioe:

|A− A| < tol.

possiamo fare una stima del valore esatto:

A− tol ≤ A ≤ A + tol.

Se conosciamo una maggiorazione dell’errore relativo, cioe:

|A− A||A| < tol.

possiamo fare una stima del valore esatto:

A

1 + tol≤ A ≤ A

1− tol.

Se riteniamo accettabili approssimazioni in cui

|A− A||A| < 0.001,

1.3. SORGENTI DI ERRORI 5

A A Errore Assoluto Errore Relativo

1 0.99 0.01 0.011 1.01 0.01 0.01-1.5 -1.2 0.3 0.2100 99.99 0.01 0.0001100 99 1 0.01

Tabella 1.1: Esempi di errore assoluto e relativo.

allora siano A = 123457 e A = 123500, calcoliamo l’errore relativo:

43

123457= 0.00034,

e poiche l’errore e minore di 0.001, l’approssimazione e accettabile.

Siano invece A = 341.5 e A = 300, l’errore relativo e:

41.5

341.5= 0.121,

e quindi l’approssimazione non e accettabile.

Se due quantita sono approssimativamente uguali, useremo la notazione ≈ per indicarequesta relazione.

Questa e una notazione ambigua. E vero che 0.99 ≈ 1? Forse si. E vero che 0.8 ≈ 1?Forse no. Sia h un parametro reale che tende a zero tale che limh→0 Ah = A allora,

Ah ≈ A

per ogni h “sufficientemente piccolo”.

Sia n un parametro intero che tende all’infinito tale che limn→∞ An = A allora,

An ≈ A

per ogni n “sufficientemente grande”.

Un modo per scrivere la derivata prima di una funzione e:

f ′(x) = limh→0

f(x + h)− f(x)

h

Possiamo quindi concludere che per h sufficientemente piccolo

f(x + h)− f(x)

h≈ f ′(x)

L’uguaglianza approssimata verifica le proprieta transitiva, simmetrica e riflessiva:

A ≈ B,B ≈ C → A ≈ CA ≈ B → B ≈ A

A ≈ A


Un’altra notazione e la notazione dell”’O grande”, conosciuta come ordine asintotico.Supponiamo di avere un valore y e una famiglia di valori che lo approssimano yh, Se esisteuna costante C > 0, indipendente da h, tale che:

|y − yh| ≤ C|β(h)|,

per h sufficientemente piccolo, allora diciamo che:

y = yh + O(β(h)) per h→ 0,

cioe y−yh e dell’ordine di β(h), β(h) e una funzione del parametro h tale che limh→0 βh = 0.Ci concentriamo sul modo in cui l’errore dipende dal parametro h e ignoriamo dettagli menoimportanti come il valore di C.

L’utilizzo e analogo se abbiamo una successione xn che approssima x per valori di ngrandi:

|x− xn| ≤ C|β(n)|, x = xn + O(β(n))

Teorema 1.3.1 Teorema di Taylor Sia f(x) una funzione avente n+1 derivate continuesu [a, b] per qualche n ≥ 0, e siano x,x0 ∈ [a, b]. Allora

f(x) = pn(x) + Rn(x)

con

pn(x) =n

∑

k=0

(x− x0)k

k!f (k)(x0)

e

Rn(x) =1

n!

∫ x

x0

(x− t)nf (n+1)(t)dt.

Inoltre esiste un punto ξx tra x e x0 tale che:

Rn(x) =(x− x0)

n+1

(n + 1)!f (n+1)(ξx)

Ritorniamo all’approssimazione della derivata prima di una funzione e vediamo quale el’errore che commettiamo se la approssimiamo con il rapporto incrementale. Sappiamo che:

f ′(x0) ≈f(x0 + h)− f(x0)

h= f ′

h(x0)

per h sufficiantemente piccolo. Vogliamo calcolare come f ′h(x0) si avvicina a f ′(x0). Usiamo

il Teorema di Taylor, con n = 1, x = x0 + h per esprimere f(x0 + h):

f(x0 + h) = f(x0) + hf ′(x0) +h2

2f ′′(ξh)

quindi

f ′h(x0) =

f(x0 + h)− f(x0)

h

=f(x0) + hf ′(x0) + h2

2 f ′′(ξh)− f(x0)

h= f ′(x0) + h

2f ′′(ξh) = f ′(x0) + O(h).

1.4. RAPPRESENTAZIONE DEI NUMERI 7

La quantita τ(h) =h

2f ′′(ξ) si chiama errore di troncamento o errore di discretizzazione

e dipende da h. Possiamo dire che l’errore va a zero come O(h)

f ′(x0) ≈f(x0 + h)− f(x0)

h= f ′

h(x0)

Questa formula per approssimare la derivata si chiama metodo alle differenze in avanti.

Consideriamo adesso il rapporto:

f(x + h)− f(x− h)

2h

Utilizzando il Teorema di Taylor abbiamo:

f(x + h)− f(x− h) =

= f(x) + hf ′(x) +h2

2f ′′(x) +

h3

6f ′′′(x) +

h4

24f (4)(ξ1)

−f(x) + hf ′(x)− h2

2f ′′(x) +

h3

6f ′′′(x)− h4

24f (4)(ξ2) =

= 2hf ′(x) +h3

3f ′′′(x) +

h4

24(f (4)(ξ1)− f (4)(ξ2))

In definitiva abbiamo

f(x + h)− f(x− h)

2h=

= f ′(x) +h2

6f ′′′(x) +

h3

48(f (4)(ξ1)− f (4)(ξ2)) =

= f ′(x) + O(h2)

Questa tecnica per approssimare la derivata prima si chiama metodo delle differenze centrali.In questo caso l’errore va a zero come O(h2)

1.4 Rappresentazione dei numeri

L’utilizzo in modo corretto del calcolatore per fare calcoli di tipo scientifico, richiede laconoscenza di come sono rappresentati i numeri e degli errori che derivano da questa rap-presentazione che sono chiamati errori di arrotondamento. L’uso dei numeri reali richiedeuna attenzione particolare, essendo questi infiniti, mentre il calcolatore ci da la possibilitadi rappresentarne solo un numero finito.

La nostra notazione per rappresentare i numeri e una notazione posizionale a base 10.Cio significa che se scriviamo 123 intendiamo esprimere il numero:

1 · 102 + 2 · 101 + 3 · 100.


Cioe la cifra che occupa la posizione piu a destra deve essere moltiplicata per 100 (=1), quella che occupa la posizione centrale per 10, ecc. Questa convenzione e ormai uni-versalmente adottata, al contrario di altre “convenzioni” che quotidianamente usiamo perla comunicazione di nostre idee ed esperienze come, ad esempio, le lingue, gli alfabeti, ledirezioni di scrittura, la musica, ecc.

Proprio per la sua universale adozione, la notazione posizionale a base 10 puo apparire amolti come l’unica possibile. Cio non e vero. Anzi, dal punto di vista storico risulta che talerappresentazione si e affermata solo da pochi secoli. In precedenza erano state usate tantealtre rappresentazioni, alcune in base 10 ma non posizionali (ad esempio quella romana),altre posizionali ma non a base 10 (ad esempio la sumerica). Torniamo alla rappresentazioneattuale.

Un qualunque numero lo possiamo esprimere nella forma:

±αpαp−1 . . . α0.α−1α−2 . . . α−q . . . (1.1)

sottointendendo:

±αp · 10p + . . . + α0 · 100 + α−1 · 10−1 + . . . + α−q · 10−q . . .

I coefficienti αi sono uno dei seguenti simboli: 0, 1, 2, ..., 9. Nulla vieta di usare al postodi 10 un altro numero, diciamo β. In tale caso la stringa (intesa come successione ordinatadi simboli) descritta dalla (1.1) assumerebbe un altro significato. Infatti, supposto che isimboli α siano compresi tra 0 ed β − 1, essa sarebbe equivalente a:

±αp · βp + . . . + α0 · β0 + α−1 · β−1 + . . . + α−q · β−q . . .

Si ottiene cosı la rappresentazione posizionale di un numero in base β.

Nel campo scientifico si preferisce usare una notazione teoricamente equivalente alla(1.1), cioe:

±d0.d1d2 . . . dt . . . · βq (1.2)

con d0 6= 0. Quindi il numero 3.57 · 102 e espresso in notazione esponenziale normalizzata,mentre i numeri 73.54 ·102 o 0.57 ·10−3 non lo sono. Con questa rappresentazione si possonorappresentare tutti i numeri eccetto lo zero. Essa non e univoca. Ad esempio in base 10 idue numeri:

1.200.....1.199.....

rappresentano lo stesso numero. La parte d0.d1d2 . . . dt . . . viene chiamata mantissa e noi laindicheremo con m. Il numero q viene chiamato esponente, mentre β e la base. Ogni numeroreale non nullo sara quindi indicato in questa rappresentazione, detta rappresentazionenormalizzata, con:

±m · βq

1.4. RAPPRESENTAZIONE DEI NUMERI 9

Il numero zero sara rappresentato con 0. Ovviamente nella rappresentazione di unnumero intero tutti i valori di saranno nulli da un certo indice (maggiore di uno) in poi. Eovvio che:

1 ≤ m < β. (1.3)

Una stessa stringa rappresenta numeri diversi a seconda della base usata. Ad esempiola stringa 123 in base quattro rappresenta il numero:

1 · 42 + 2 · 4 + 3 = ventisette

cioe il numero che in base 10 indichiamo con 27. Nel caso in cui β fosse piu grande di dieci, ilnumero di simboli che usualmente si usano per indicare le cifre tra 0 ed β−1, cioe 0,1,2,. . .,9non bastano piu. Si fa ricorso alle lettere maiuscole dell’alfabeto per integrare i simbolimancanti. Ad esempio se β = 16, i simboli usati sono 0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F.Intendendo con questo che il numero dieci viene rappresentato con A, l’undici con B, ecc..Il numero sedici sara ovviamente indicato con 10.

Esiste una convenienza nell’uso di una base al posto di un’altra? Diamo uno sguardoalla tabella seguente in cui il numero (non la stringa!) centoventitre e rappresentato indiverse basi:

base rappresentazione

16 7B10 1235 4434 13232 1111011

E evidente che lo stesso numero richiede un numero di cifre crescente al diminuire dellabase. Cio puo risultare scomodo. Vi e pero certamente un vantaggio nell’usare le basi piupiccole. Infatti i calcoli risultano semplificati. La tavola pitagorica, che in base dieci eformata da un quadrato di 100 elementi, nella base 2 risulterebbe formata da 4 elementi.

La base due e ancora piu interessante perche in tale base tutti i numeri sono rappresenta-ti da successioni dei due simboli 0 ed 1 che possono essere messe in corrispondenza con ognisistema che abbia solo due stati (ad esempio acceso-spento, magnetizzato-smagnetizzatoecc.). Difatti e per questo motivo che essa risulta indispensabile per la codifica di infor-mazioni nei calcolatori. Ogni informazione, per essere memorizzata nel calcolatore, vieneopportunamente codificata nel linguaggio comprensibile alla macchina che e appunto ba-sato sui simboli del sistema di numerazione binario, 0 ed 1. Abbiamo gia avuto occasionedi dire che questi due simboli rappresentano lo stato fisico dei circuiti di cui e costituita lamemoria, che possono essere spenti (=0) o accesi (=1).

Ovviamente la rappresentazione (1.2) necessita un numero infinito di cifre e quindisarebbe impossibile memorizzarla in un calcolatore. Nei calcolatori ad ogni numero vieneriservato uno spazio fisso di memoria, per cui sia il numero di cifre della mantissa che


il numero di cifre dell’esponente non possono superare dei limiti prefigurati. I numeri dimacchina, detti anche numeri floating point, sono del tipo:

±d0.d1d2 . . . dt · βq (1.4)

con d0 6= 0, 0 ≤ di ≤ β − 1 e qmin ≤ q ≤ qmax, piu, naturalmente, lo zero. Denotiamocon realmin il piu piccolo numero positivo e con realmax il piu grande numero positivorappresentabile.

Se per esempio t = 1, β = 10, qmin = −2, qmax = 2 allora:

+1.0 · 10−2 e il piu piccolo numero positivo rappresentabile in notazione esponenzialenormalizzata. I numeri successivi sono 1.1 · 10−2, 1.2 · 10−2, 1.3 · 10−2, . . . 2.0 · 10−2, . . . , Ilpiu grande numero rappresentabile in questa notazione e 9.9 · 102. Quanti numeri possiamorappresentare in totale? Ogni numero e individuato dal segno, dalle due cifre d0, d1 e dall’e-sponente e. Ci sono 9 possibilita per d0 ∈ 1, 2, 3, . . . , 9, 10 possibilta per d1 ∈ 0, 1, , 2, . . . , 9e 5 possibilita per q ∈ −2,−1, 0, 1, 2. Quindi in totale abbiamo 900 numeri piu lo zero.

La quantita che devo aggiungere ad 1.0 per avere il successivo numero di macchina eǫM = 0.1, tale numero si chiama ǫ di macchina.

Un numero reale rappresentato tramite la (1.2) deve essere prima espresso in cifre binarie0 ed 1 per poter essere memorizzato in un calcolatore.

1.5 Errore assoluto e relativo

Consideriamo il numero reale positivo:

x = d0.d1d2 . . . dt . . . · βq

Se consideriamo i numeri di macchina con t cifre per la mantissa, x sara compreso fra:

x1 = d0.d1d2 . . . dt · βq

e

x2 = x1 + 0.0 . . . 01 · βq = x1 + β−t · βq

graficamente si ha:

d0.d1d2 . . . dt . . . · βq

| | |x1 x2

Vi sono due modi diversi per approssimare questo numero mediante i numeri floatingpoint, in cui solo t cifre della mantissa sono rappresentabili:

1. il troncamento;

2. l’arrotondamento.

1.5. ERRORE ASSOLUTO E RELATIVO 11

x troncato arrotondatot=2 t=2

5.672 5.67 5.67-5.672 -5.67 -5.675.677 5.67 5.68-5.677 -5.67 -5.68

Tabella 1.2: Esempi di troncamento e arrotondamento

Nel primo caso il numero x viene sempre approssimato con x1 trascurando tutte le cifresuccessive a t; nel secondo caso il numero viene approssimato con x1 se x < (x1 + x2)/2,cioe se x si trova nella prima meta dell’intervallo [x1x2], altrimenti viene approssimato conx2. Se x = (x1 + x2)/2 allora x viene approssimato con x1 se il numero intero d0 . . . dt

e pari, altrimenti viene approssimato con x2. Vedi la Tabella 1.2 con alcuni esempi diarrotondamento e troncamento.

Denotiamo con fl(x) l’approssimazione di x. Se viene usato il troncamento, qualunquesia la posizione di x nell’intervallo, x viene approssimato con il numero di macchina allasua sinistra e quindi l’errore assoluto, cioe la distanza in valore assoluto di fl(x) da x, sarasempre minore dell’ampiezza dell’intervallo [x1, x2]:

|fl(x)− x| < β−t · βq

Nel secondo caso invece l’errore assoluto sara minore della meta dell’ampiezza dell’intervallo:

|fl(x)− x| < 1

2β−t · βq

L’errore relativo di un numero x 6= 0, approssimato con fl(x) e definito da:

∣

∣

∣

∣

fl(x)− x

x

∣

∣

∣

∣

Teorema 1.5.1 Se x 6= 0 e usiamo t cifre per rappresentare la mantissa, allora:

∣

∣

∣

∣

fl(x)− x

x

∣

∣

∣

∣

≤ u (1.5)

dove:

u =

{

β−t in caso di troncamento12β−t in caso di arrotondamento

Il numero u e detto unita di arrotondamento.

Dimostrazione. Sia x = d0.d1d2 . . . dt . . . · βq, allora |x| ≥ βq e, nel caso di troncamento,si ha:

|fl(x)− x||x| ≤ βq−t

βq≤ β−t,


da cui segue la tesi.

Si procede analogamente nel caso di arrotondamento e nel caso di numeri negativi.

Come conseguenza del Teorema 1.5.1 si ha che l’errore relativo, piuttosto che l’erroreassoluto, e legato alle cifre esatte di un numero x. Se |x| < realmin, allora si ha un erroredetto errore di “underflow”, oppure fl(x) si pone uguale a 0 e quindi:

∣

∣

∣

∣

x− fl(x)

x

∣

∣

∣

∣

= 1,

cioe l’errore relativo e molto piu grande di u. Se |x| > realmax il numero non e rappresen-tabile e si genera un errore detto di “overflow”.

Poniamo:

ǫ =fl(x)− x

x

Sappiamo che |ǫ| ≤ u. Dalla precedente si ricava:

fl(x) = x(1 + ǫ), (1.6)

cioe per ogni x tale che realmin ≤ |x| ≤ realmax esiste un numero ǫ in valore assolutominore od uguale ad u tale che la sua approssimazione mediante un numero di macchinapuo essere espressa dalla (1.6). Se vale la (1.5) possiamo dire che il numero e rappresentatocon t cifre esatte.

1.6 Standard IEEE

Descriviamo piu in dettaglio come viene rappresentato un numero usando lo standard floa-ting point IEEE, che costituisce un insieme di regole definito dall’istituto degli ingegnerielettrici e elettronici per la rappresentazione ed l’elaborazione dei numeri floating point neicomputer. Lo standard IEEE specifica esattamente cosa sono i numeri floating point e comesono rappresentati nell’hardware e ha 4 scopi principali:

1. rendere l’aritmetica floating point il piu accurata possibile;

2. produrre risultati sensati in situazioni eccezionali;

3. standardizzare le operazioni floating point fra i computer;

4. Dare al programmatore un controllo sulla manipolazione delle eccezioni;

Il punto (1) si ottiene in due modi. Lo standard specifica esattamente come un nume-ro floating point dovrebbe essere rappresentato nell’hardware e richiede che le operazioni(addizione,radice,...) siano il piu accurate possibili. Per il punto (2) lo standard produceInf (Infinito) per indicare che un risultato e piu grande del massimo numero floating pointrappresentabile e NaN (Not a Number) per indicare che un risultato non e definito. Primadello standard la maggior parte dei computer avrebbero bloccato un programma in talicircostanze. Il punto (3) ha alcune conseguenze:

1.6. STANDARD IEEE 13

(i) I numeri floating point possono essere trasferiti fra due computer che utilizzano lostandard IEEE in binario senza perdere precisione;

(ii) I dettagli dell’aritmetica floating point dovrebbero essere conosciuti dal programma-tore;

(iii) Lo stesso programma su computer diversi dovrebbe produrre lo stesso risultato, anchese questo non e vero in pratica.

Ricordiamo che l’unita base di informazione immagazzinata da un computer e il bit, unavariabile il cui valore e o 0 o 1. I bit sono organizzati in gruppi di 8 chiamati byte. L’unitapiu comune e la parola di 4 byte (32 bit). Per accuratezze piu alte e raddoppiata a 8 byteo 64 bit. Vi sono 2 possibili valori per un bit, 28 = 256 possibili valori per un byte e 232

differenti parole di 4 byte.

I due tipi di numeri rappresentati sono interi (fixed point) e reali (floating point). Unnumero reale ha tipo float in c e real in Fortran e Matlab. Nella maggior parte dei compi-latori c e in Matlab un float ha per default 8 byte invece di 4. Il formato IEEE sostituiscela base 10 con la base 2.

L’aritmetica con gli interi e molto semplice. Vi sono 232 ≈ 4 · 109 interi a 32 bit checoprono l’intervallo da −2 · 109 a 2 · 109. Addizione, sottrazione e moltiplicazione sonofatte esattamente se la risposta e compresa nell’intervallo. La maggior parte dei computerdanno risultati imprevedibili se il risultato e fuori dal range (overflow). Lo svantaggiodell’aritmetica con gli interi e che non possono essere rappresentate le frazioni e l’intervallodei numeri e piccolo.

Quando una parola a 32 bit e interpretata come un numero floating point il primo bite il bit del segno, s = + o s = −. I successivi 8 bit formano l’esponente e i rimanenti23 bit determinano la mantissa. Vi sono 2 possibili segni, 256 esponenti (che variano da0 a 255) e 223 ≈ 8.4 milioni di possibili mantisse. Per avere gli esponenti negativi si usauna convenzione: lo zero e nella posizione 127, dopo ci sono i numeri positivi e prima inumeri negativi. Quindi in memoria viene rappresentato q∗ = q + 127. Inoltre il primobit della mantissa, che rappresenta d0, e sempre uguale a 1, quindi non vi e necessitadi memorizzarlo esplicitamente. Nei bit assegnati alla mantissa viene memorizzato m∗ em = 1.m∗. Un numero floating point positivo ha quindi il valore x = +(1.m∗)22

q∗−127 e lanotazione (1.m∗)2 indica che 1.m∗ e interpretata in base 2.

Per esempio il numero 2.752 · 103 = 2572 puo essere scritto:

2752 = 211 + 29 + 27 + 26 == 211(1 + 2−2 + 2−4 + 2−5) == 211(1 + (0.01)2 + (0.0001)2 + (0.00001)2)) == 211(1.01011)2

allora la rappresentazione di questo numero avrebbe segno + esponente q∗ = q + 127 =138 = (10001010)2 e m∗ = (010110 . . . 0)2.

Il caso q∗ = 0 ( che corrisponde a 2−127) e il caso q∗ = 255 (che corrisponde a 2128 hannouna interpretazione differente e complessa che rende la IEE diversa dagli altri standard.

Se q∗ = 0 il valore del numero memorizzato e x = ±(0.m∗)22−126. Questo e chiamato

underflow graduale (l’underflow e la situazione in cui il risultato di una operazione e diversa


da zero ma e piu vicina a zero di qualsiasi numero floating point). I numeri corrispondentivengono chiamati denormalizzati. L’underflow graduale ha la conseguenza che due numerifloating point sono uguali se e solo se sottraendone uno dall’altro si ha esattamente zero.

Se escludiamo i numeri denormali allora il piu piccolo numero positivo floating point ea = realmin = q−126, il successivo numero positivo b ha q∗ = 1 e m∗ = 0.0 · · · 01, la distanzafra a e b e 223 volte piu piccola, della distanza fra a e zero. Senza l’underflow graduale vi eun gap fra zero ed il numero floating point piu vicino.

L’altro caso e q∗ = 255 che ha due sottocasi, Inf se m∗ = 0 e NaN se m∗ 6= 0. Siail c che il Matlab stampano Inf o NaN quando si stampa una variabile floating point checontiene questo risultato. Il computer produce Inf se il risultato di una operazione e piugrande del piu grande numero rappresentabile.

Operazioni invalide che producono come risultato NaN sono: Inf/Inf , 0/0, Inf−Inf .Operazioni con Inf hanno il significato usuale: Inf + finito = Inf , Inf/Inf = NaN ,Finito/Inf = 0, Inf − Inf = NaN .

E chiaro che l’accuratezza delle operazioni Floating Point e determinata dalla grandez-za degli errori di arrotondamento. Questo errore di arrotondamento e determinato dalladistanza di un numero dal numero floating point piu vicino. Eccetto che per i numeri denor-malizzati, i numeri floating point differiscono di un bit, l’ultimo bit di m∗. Cioe 2−23 ≈ 10−7

in singola precisione. Questo e l’errore relativo e non l’errore assoluto.

L’aritmetica IEEE in doppia precisione usa 8 byte (64 bit). Vi e un bit del segno, 11 bitper l’esponente e 52 bit per la mantissa. Quindi in doppia precisione ǫ = 2−52 ≈ 2.22·10−16 .Cioe l’aritmetica in doppia precisione ha circa 15 digit decimali di accuratezza, invece di 7.Vi sono 211 possibili esponenti che variano da 1023 a -1022. Il piu grande numero e circa10307 ≈ 21023

Riassumendo:

IEEE - singola precisione. Occupa 4 byte = 32 bits

←− 8 −→ ←− 23 −→ —

s q m

segno esponente mantissa

Esso rappresenta (−1)s · 2q−127 · (1.m) Si noti che d0 = 1 nella mantissa nondeve essere esplicitamente memorizzato, quindi m = d1d2 . . . dt . Range di numeripositivi normalizzati da 2−126 a 2127×1.11 · · · 1 ≈ 2128 (da 1.2×10−38 a 3.4×1038)

1.7. OPERAZIONI CON I NUMERI DI MACCHINA 15

IEEE - doppia precisione. Occupa 8 byte = 64 bits

←− 11 −→ ←− 52 −→ —

s q m

segno esponente mantissa

Esso rappresenta (−1)s · 2q−1023 · (1.m)Range di numeri positivi normalizzati da 2−1022 a 21023 × 1.11 · · · 1 ≈ 21024 (da2.2× 10−308 a 1.8× 10308)

IEEE - Eccezioni aritmetiche e risultati di default

tipo di esempio risultatoeccezione di default

operazione 0/0 0×∞ NaNinvalide

√−1 (Not a number)

Overflow ±∞

Divisione Numero finito non nullo/0 ±∞per zero

Underflow numeridenormali

1.7 Operazioni con i numeri di macchina

Sia o ∈ {+,−,×, /} e siano x e y due numeri floating point. E improbabile che l’esattovalore di xoy sia un numero floating point.

Esempio t = 3, β = 10

1.111 · 103 × 1.111 · 102 = 1.234321 · 103

che richiede piu di tre cifre decimali.

Il computer dovrebbe eseguire le operazioni aritmetiche di base in modo che il risultatofinale sia il risultato esatto arrotondato al piu vicino numero floating point. Il modello delleoperazioni aritmetiche richiesto dall’IEEE e il seguente:

fl(xoy) = (xoy)(1 + ǫ), |ǫ| ≤ u

ofl(xoy) = (xoy)/(1 − ǫ), |ǫ| ≤ u


L’aritmetica floating point IEEE richiede questo. Supponiamo di voler eseguire l’addizionetra due numeri x = α0.α1α2 . . . αt ·βp e y = d0.d1d2 . . . dt ·βq, in cui p > q. Dobbiamo trasfo-mare i due numeri in modo da avere i due esponenti uguali a p. Per fare cio bisogna che l’areadi memoria in cui si eseguono le operazioni abbia uno spazio riservato alle mantisse superiorea quello usuale. Nel nostro caso si pone y = (d0.d1d2 . . . dt ·βq−p) ·βp = .0 . . . 0d0d1 . . . dt ·βp

e poi si sommano le mantisse. Si procede quindi al troncamento o all’arrotondamento.In generale quindi il risultato di una operazione puo considerarsi come il risultato esattodell’operazione tra x ed y seguito dall’arrotondamento o dal troncamento.

Le operazioni con i numeri di macchina (aritmetica floating point), non godono di tuttele proprieta di cui godono le corrispondenti operazioni con i numeri reali. In generale nonvale piu la proprieta associativa. Consideriamo ad esempio il seguente caso in cui β = 10 et = 2. Eseguiamo la somma 5.24 · 10−2 + 4.04 · 10−2 + 1.21 · 10−1 in due modi diversi:

1) 5.24 · 10−2 + (4.04 · 10−2 + 1.21 · 10−1) = 5.24 · 10−2 + (0.404 + 1.21)10−1 = 5.24 · 10−2 +1.61 · 10−1 = (0.524 + 1.61)10−1 = 2.13 · 10−1;

2) (5.24 · 10−2 + 4.04 · 10−2) + 1.21 · 10−1 = 9.28 · 10−2 + 1.21 · 10−1 = (0.928 + 1.21)10−1 =2.14 · 10−1.

Supponiamo ora di avere due numeri reali x ed y e di volere fare la somma x + y. Nelcalcolatore x sara rappresentato da fl(x) = x(1 + ǫx) e y da fl(y) = y(1 + ǫy). Inoltrefl(fl(x) + fl(y)) = (fl(x) + fl(y))(1 + ǫ). Quindi fl(fl(x) + fl(y)) = (x(1 + ǫx) + y(1 +ǫy))(1 + ǫ) con |ǫ|, |ǫx|, |ǫy| ≤ u

Calcoliamo l’errore relativo trascurando i termini contenenti ǫǫx ed ǫǫy,

|(x + y)− (fl(x) + fl(y))(1 + ǫ)||x + y|

=|x + y − (x + xǫx + y + yǫy)(1 + ǫ)|

|x + y| ≈

≈ |ǫ|+ |ǫx||x||x + y| + |ǫy|

|y||x + y| .

Se x + y non e piccolo l’errore relativo e dello stesso ordine degli errori relativi |ǫ|, |ǫx|ed |ǫy|.Se x + y e molto piccolo l’errore relativo e grande.

Dal precedente risultato risulta evidente che se il denominatore non e piccolo, comeavviene certamente nel caso in cui i due numeri siano dello stesso segno, l’errore relativo edello stesso ordine degli errori relativi ǫ, ǫ1 ed ǫ2 e quindi in valore assoluto minore di u.Cio non avviene nel caso in cui il denominatore risulti molto piccolo. Cio risultera evidentenel seguente esempio.

Supponiamo infatti di voler effettuare la differenza tra i due numeri reali x1 = 0.147554326ed x2 = 0.147251742 ed operare con t = 5 ed β = 10. Usando l’arrotondamento sarafl(x1) = 1.47554 ·10−1 e fl(x2) = 1.47252 ·10−1 . Si ha fl(fl(x1)−fl(x2)) = 3.02000 ·10−4 ,mentre la vera differenza e 3.02584 ·10−4 . Le ultime tre cifre della mantissa risultano errate.Cio e dovuto al fatto che dopo aver eseguito la differenza fl(x1) − fl(x2) = 0.000302, larappresentazione normalizzata 3.02000 · 10−4 ha introdotto tre zeri alla fine della mantissa.Ovviamente questi tre zeri non sono esatti. L’errore relativo risultante e:

1.7. OPERAZIONI CON I NUMERI DI MACCHINA 17

3.02584 − 3.02000

3.02584≈ 10−3

che e piuttosto elevato. Quindi bisogna fare particolare attenzione a questo fenomeno, cheviene chiamato cancellazione numerica. Essa si presenta quando si esegue la sottrazione tranumeri molto vicini.

Allo stesso modo si puo analizzare il prodotto:

fl(fl(x) ∗ fl(y)) = fl(x(1 + ex) ∗ y(1 + ey)) = x(1 + ex)y(1 + ey)(1 + e∗)

l’errore relativo e:

|x(1 + ex)y(1 + ey)(1 + e∗)− xy|/|xy|

cioe

|(1 + ex)(1 + ey)(1 + e∗)− 1|

semplificando ed eliminando i termini che contengono i prodotti di pi‘uerrori si ottiene:

errore relativo nel prodotto ≈ |ex|+ |ey|+ |e∗|

Spesso e comodo vedere la (2.7) da un punto di vista leggermente diverso. Consideriamoper esempio il caso in cui op sia l’addizione. Sara:

fl(x + y) = (x + ǫx) + (y + ǫy),

cioe il risultato dell’operazione di macchina puo essere visto come l’operazione esatta sudati perturbati (x+δx) e (y+δy). Ovviamente nel caso considerato sara δx = ǫx e δy = ǫy.Nel caso in cui op sia la moltiplicazione si ha:

fl(xy) = xy(1 + ǫ) = (x + δx)(y + δy) con δx = 0, δy = ǫy.

La tecnica di considerare la perturbazione nei dati che porterebbe allo stesso risultatofinale se le operazioni fossero eseguite con precisione infinita e molto utile e viene chiamataanalisi degli errori all’indietro (backward).

Esercizio 1.7.1 Sia t = 3 e supponiamo che le operazioni siano eseguite con t = 6 cifredecimali, eseguire le seguenti operazioni aritmetiche dopo over trasformato i numeri nellarappresentazione floating point:

1.235 + 1.527 · 10−4,

.5 · 103 − 2.5,

1.7435 · 10−5 − 1.7448 · 10−5.

Calcolare l’errore relativo e l’errore assoluto commesso nell’eseguire i calcoli.


1.8 Propagazione degli errori e condizionamento

Un problema numerico e una chiara e non ambigua descrizione della connessione funzionalefra i dati di input, variabili indipendenti del problema, e i dati di output, cioe i risultatidesiderati. I dati di output devono essere determinati in modo univoco e devono dipendere inmodo continuo dai dati di input. Un algoritmo e una descrizione completa di operazioni bendefinite attraverso cui ogni insieme di dati di input ammissibile e trasformato nell’insiemedei dati di output. Le operazioni sono operazioni aritmetiche e logiche che un computerpuo eseguire.

La propagazione degli errori e di grande interesse per pianificare e analizzare esperimentiscientifici. Anche gli errori di arrotondamento di ogni passo di un calcolo si propagano sulrisultato finale. Per molti algoritmi si puo fare un’analisi degli errori di arrotondamento incui si mostra che il risultato finale e il risultato esatto di un problema in cui i dati di inputsono stati perturbati di poco.

Consideriamo, per esempio, una funzione y = f(x) e supponiamo di volerla calcolarein un punto particolare x il quale e stato misurato con un errore ∆x. Questo errore siripercuote sulla variabile y con un errore ∆y che, nel caso in cui f ammetta derivata primain un intorno di x si puo stimare facilmente. Si ha:

y + ∆y = f(x + ∆x) ≃ f(x) + f ′(x)∆x

da cui:∆y

y≃ x

f ′(x)

f(x)

∆x

x.

Si deduce che un errore relativo |∆x|/|x| sulla variabile indipendente produce un errorerelativo sulla variabile dipendente che risultera notevolmente amplificato se la quantita:

K(x, f) =

∣

∣

∣

∣

xf ′(x)

f(x)

∣

∣

∣

∣

e molto grande. Questa quantita e detta indice di condizionamento o numero di condizionedel problema.

Definizione 1.8.1 Dato un problema numerico y = f(x) ∈ Rm, x ∈ Rn, sia x fissato e as-sumiamo che x, f(x) siano diversi da zero. La sensibilita di y rispetto a piccoli cambiamentidi x si misura utilizzando il numero di condizione relativo:

K(f, x) = limǫ→0

sup‖h‖ǫ

{‖f(x + h)− f(x)‖‖f(x‖ /

‖h‖‖x‖

}

dove ‖ · ‖ e una norma vettoriale,per esempio

‖x‖p = (|x1|p + |x2|p + . . . |xn|p)(1/p), p = 1, 2

‖x‖∞ = max1≤i≤n|xi|

1.8. PROPAGAZIONE DEGLI ERRORI E CONDIZIONAMENTO 19

Il numero di condizione relativo misura la massima amplificazione di una perturbazione suidati di input

‖y − y‖ ≤ Kǫ‖y‖+ O(ǫ2)

quindi ci dobbiamo aspettar di avere circa s = log10(K) cifre significative in meno nellasoluzione.

Se il numero di condizione e molto grande il problema si dice mal condizionato, se e digrandezza moderata si dice ben condizionato. Cosa significa molto grande? Dipende dal-l’accuratezza che richiediamo in output quando risolviamo un problema e dalla grandezzadelle peturbazioni. Se ǫ = 10−16 e richiediamo 8 cifre significative corrette nella soluzio-ne, allora K deve essere minore di 108 altrimenti il problema non e risolvibile nei limitidell’accuratezza desiderata.

Se otteniamo una soluzione poco accurata per un problema mal condizionato non pos-siamo farci niente. Se fai una domanda sciocca, ottieni una risposta sciocca.

Esempio 1.8.1 Sia f(x) = log x. Si ha K(x, f) = 1/ log(x) da cui risulta che il calcolo dellogaritmo e un problema mal condizionato se x e vicino ad 1.

Esempio 1.8.2 Consideriamo il problema di calcolare le radici di un polinomio. In questocaso i dati del problema sono i coefficienti del polinomio e il risultato finale le sue radici.Se consideriamo il seguente polinomio:

(x− 1)4 = x4 − 4x3 + 6x2 + 4x + 1

sappiamo che la radice e 1 con molteplicita 4. Perturbiamo il termine noto, ottenendo ilpolinomio:

(x− 1)4 = x4 − 4x3 + 6x2 + 4x + 1− 10−8

il quale ha le radici:

x1 = 1.01x2 = 0.99

x3 = 1 + i0.01x4 = 1− i0.01

quindi una pertubazione relativa di 1e-8 sul termine noto produce una perturbazione relativadi 1e-2 sul risultato. Possiamo concludere che il problema e mal condizionato.

Quando il risultato di un problema non si ottiene soltanto calcolando il valore di unafunzione regolare in un punto, lo studio della propagazione degli errori e del condizionamentodel problema puo essere molto laboriosa.

Studiare il condizionamento del problema e molto importante nel caso in cui si vogliarisolvere il problema usando il calcolatore. In questo caso infatti una perturbazione sui dativiene sempre provocata quando gli stessi vengono memorizzati.


1.9 Analisi della stabilita degli algoritmi

Di solito la poca accuratezza del risultato dipende dal mal condizionamento del problema, avolte, pero, dipende dall’algoritmo. Ci sono algoritmi che eseguiti in aritmetica di macchinafanno propagare in maniera molto elevata gli errori di arrotondamento di ogni istruzione.

Definizione 1.9.1 Un algoritmo si dice instabile se introduce grandi errori quando e ap-plicato a un problema ben condizionato.

Consideriamo un algoritmo con dati di input x1, x2, . . . , xr, x = (x1, x2, . . . , xr)T e dati

di output y1, y2, . . . , yn, y = (y1, y2, . . . , yn)T e siano y1, y2, . . . , yn, y = (y1, y2, . . . , yn)T

l’output generato dall’algoritmo eseguito sul calcolatore. Possiamo fare due diversi tipidi analisi per studiare la stabilia di un algoritmo l’analisi forward (in avanti) e l’analisibackward (all’indietro).

• L’analisi Forward consiste nel trovare una maggiorazione dell’errore che si produ-ce sull’output andando ad analizzare gli errori che si introducono nei singoli passidell’algoritmo. Si determina, quindi una maggiorazione della quantita

‖y − y‖‖y‖ ≤ CF u, ‖y‖ 6= 0

con u l’unita di arrotondamento e CF una costante non molto grande.

• L’analisi Backward consiste nel trovare x che, quando l’algoritmo viene eseguito in arit-metica esatta, genera y come risultato. Quindi si trova una maggiorazione dell’errorerelativo:

‖x− x‖‖x‖ ≤ CBu, ‖x‖ 6= 0

con u l’unita di arrotondamento e CB una costante non molto grande.

Come sono legate l’analisi forward e quella backward? Se K e il numero di condizionedel problema si ha che

‖y − y‖‖y‖ ≤ K

‖x− x‖‖x‖ ≤ KCBu

ne segue che un algoritmo che e backward stabile e anche forward stabile. Ci sono peroalcuni problemi per cui non si puo fare l’analisi backward. Un esempio e il prodotto esternofra vettori:

A = xyT =

x1y1, . . . x1yn

x2y1, . . . x2yn

... . . ....

xry1 . . . xryn

.

Quando eseguiamo il prodotto esterno al calcolatore abbiamo un errore per ogni elemento

fl(xiyj) = xiyj(1 + δij) con |δij | < u

1.9. ANALISI DELLA STABILITA DEGLI ALGORITMI 21

La matrice A ottenuta utilizzando l’aritmetica di macchina avra elementi

Aij = xiyj + xiyjδij

dunqueA = A + ∆

con∆ij = xiyjδij

non e possibile trovare delle perturbazioni sui dati di input che generano A. Quando l’in-sieme dei dati di input e piu grande dell’insieme dei dati di output e molto probabile chel’analisi backward non si possa effettuare. In tal caso si esegue un’analisi mista.

Capitolo 2

Laboratorio Matematico

Informatico: Matlab

Questo capitolo riassume le lezioni di laboratorio del corso di Laboratorio MatematicoInformatico a.a 2010-2011, docente: Francesca Mazzia

2.1 Introduzione al Matlab

MATLAB (MATrix LABoratory): PSE per il scientifico ad alte prestazioni e la visualiz-zazione, integra analisi numerica, calcolo con matrici, grafica. E dotato di un particolarelinguaggio di programmazione di tipo interpretativo che consente di sviluppare le proprieapplicazione. Ogni istruzione viene interpretata ed eseguita immediatamente.Entrati in ambiente Matlab, compare il prompt:

>>

il che significa che il calcolatore e pronto a ricevere le istruzioni Matlab ed ad eseguirle. Sevogliamo eseguire la somma 3+2 basta scrivere:

>> 3+2

e premere il tasto invio. Otteniamo immediatamente la risposta:

ans =

5

>>

e il calcolatore e di nuovo pronto a ricevere un nuovo comando. Il risultato dell’operazioneeseguita e memorizzato nella variabile ans. Le variabili nei linguaggi di programmazionesono dei contenitori di valori. Se vogliamo assegnare il risultato dell’operazione alla variabiledi nome var basta scrivere:

23

24 CAPITOLO 2. LABORATORIO MATEMATICO INFORMATICO: MATLAB

>> var = 3+2

var =

5

Per evitare la visualizzazione del risultato dell’istruzione eseguita si scrive il simbolo ; allafine della riga, per esempio:

>> var = 3+2;

Questa istruzione si chiama istruzione di assegnazione. Possiamo assegnare a una va-riabile il valore di una espressione piu complicata. Per far questo dobbiamo sapere a qualisimboli corrispondono le principali operazioni aritmetico-logiche. In Matlab questi sono:

+ addizione− sottrazione∗ moltiplicazione/ divisione^ elevamento a potenza& and logico| or logico~ not logico== uguale~ = diverso< minore<= minore uguale> maggiore>= maggiore uguale

Le espressioni aritmetiche e logiche seguono le classiche regole di precedenza, che epossibile modificare utilizzando le parentesi tonde.

Per avere informazioni su tutte le operazioni che possono essere eseguite in Matlab sipuo digitare sul pulsante help e quindi su Help Window, compare una finestra con tutte leinformazioni sul Matlab divise per capitoli. Si puo usare anche il comando help. Eseguimoi seguenti comandi:

>> help

>> help help

>> help(’=’)

>> help(’*’)

Possiamo ora eseguire:

>> var = 5*3/(3+2)^2

var =

0.6000

2.1. INTRODUZIONE AL MATLAB 25

E anche possibile calcolare i valori delle principali funzioni. Ad esempio sono predefinite lefunzioni in una variabile:

abs(x) valore assoluto di x,log(x) logaritmo naturale,exp(x) esponenziale con base e,sin(x) seno,cos(x) coseno,tan(x) tangente.

E possibile quindi eseguire espressioni del tipo:

>> var = (5* 3/(3+cos(5)))^log(2) * 100

var =

286.6104

Se eseguiamo il comando

>> format short e

e rieseguimo l’operazione

>> var = (5* 3/(3+cos(5)))^log(2) * 100

var =

2.8661e+02

La variabile var = 2.8661e+002 e ora una variabile reale espressa in notazione esponenziale,essa e equivalente ad 2.8661 · 102.

Anche per queste funzioni e possibile usare il comando help. Si provi ad eseguire:

>> help cos

>> help tan

>> help log

Oltre alle variabili, sono di particolare utilita anche le costanti. Le costanti sono deicontenitori di valori che non cambiano nel corso dell’esecuzione di un algoritmo. In Matlabvi sono delle costanti predefinite, come, ad esempio, la costante pi che contiene il valore diπ:

>> pi

ans =

3.1442


Le variabili possono essere scalari oppure vettori o matrici. I vettori e le matrici sonovariabili che contengono un certo numero di componenti. Ad esempio dopo l’istruzione:

>> A = [ 3, 5, 6];

la variabile A sara un vettore di tre elementi con A(1)=3, A(2)=5, A(3)=6. Gli elementidi un vettore si indirizzano specificando l’indice tra parentesi tonde. Se vogliamo sapere lasomma degli elementi di A possiamo eseguire:

>> A(1) + A(2) + A(3)

ans =

14

Dopo l’istruzione:

>> A = [ 3, 5, 7

8, 9, 10

2, 3, 5 ] ;

la variabile A sara una matrice con tre righe e tre colonne con elementi A(1,1)=3, A(1,2)=5,. . ., A(2,1)=8, . . ., A(3,3)=5. Ogni elemento di una matrice viene individuato da una coppiadi indici, il primo indice e detto indice di riga, il secondo indice di colonna.

Il Matlab consente anche di memorizzare una sequenza di istruzioni in un file.

I nomi dei file Matlab, detti anche M-file, devono avere l’estensione .m, cioe devono esseredel tipo nomefile.m. Essi possono essere scritti utilizzando un editore di testi richiamabilecon il comando:

>> edit

oppure selezionando con il mouse file e new se il file e nuovo o open se il file e stato scrittoprecedentemente.

Gli M-file possono essere di due tipi script e function. Gli M-file di tipo script contengonosemplicemente una sequenza di istruzioni Matlab. Essi utilizzano tutte le variabili giadefinite, dette variabili globali, le quali possono essere elencate con il comando:

>> who

oppure:


>> whos

Gli M-file di tipo function, invece, definiscono una nuova funzione Matlab. Essi accettanodei dati di input e restituiscono dei dati di output come risultato della loro elaborazione.La prima istruzione in un M-file di tipo function deve essere:

function [ o1 , o2 , . , on ] = NomeFunction ( i1 , i2 , . , m)

dove o1,o2 ,., on sono le variabili di output e i1 , i2 ,., im sono le variabili di input.

Quando abbiamo terminato di lavorare con il Matlab possiamo uscire con il comando:

>> quit

Quando si esce dal Matlab tutte le variabili definite vanno perdute. E possibile tuttaviasalvarne il contenuto in un file mediante il comando:

>> save NomeFile

Se nome-file non e specificato il contenuto e salvato nel file matlab.mat. Se nome-file eseguito da un elenco di variabili vengono memorizzate solo le variabili considerate. Quandosi entra in Matlab e possibile caricare le variabili precedentemente memorizzate in un filemediante il comando:

>> load NomeFile

Anche in questo caso se nome-file non e specificato si leggono le variabili del file matlab.mat.

Eseguiamo ora la seguente istruzione

>> 2+4

sul video compare

ans =

6

Se ora eseguiamo i comandi who e whos otteniamo

>> who

Your variables are:


ans

>> whos

Name Size Bytes Class

ans 1x1 8 double array

Grand total is 1 elements using 8 bytes

In Matlab possiamo usare anche variabili complesse

>> y= sqrt(-1)

y =

0 + 1.0000i

Se rieseguimo whos

>> whos

Name Size Bytes Class

ans 1x1 8 double array

x 1x1 8 double array

y 1x1 16 double array (complex)

Grand total is 3 elements using 32 bytes

Il comando format permette di cambiare il modo in cui il numero viene visualizzato,per esempio possiamo visualizzare la costante pi

>> pi

ans =

3.1416

>> format long

>> pi

ans =

3.14159265358979

>> format long e


>> pi

ans =

3.141592653589793e+000

>> format short e

>> pi

ans =

3.1416e+000

Una costante importante e la costante eps che contiene la precisione di macchina:

>> eps

ans =

2.2204e-016

Altre costanti collegate alla rappresentazione dei numeri sono:

>> realmin

ans =

2.2251e-308

>> realmax

ans =

1.7977e+308

Il Matlab utilizza l’aritmetica IEEE e permette di rappresentare infinito e NaN:

>> 1/0

Warning: Divide by zero.

ans =

Inf

>> x=1/0



x =

Inf

Ci sono alcune function matlab che permettono di verificare se un numero e finito o see uguale ad NaN. Eseguima le seguenti operazioni:

>> isfinite(x)

ans =

0

>> 0/0


ans =

NaN

>> y=0/0


y =

NaN

>> isnan(y)

ans =

1

Il Matlab consente di scrivere i file di tipo script e di tipo function utilizzando lasequenza, la selezione e la ripetizione. Per la selezione si usa la parola chiave if. Eseguiamo:

>> help if

IF IF statement condition.

The general form of the IF statement is

IF expression

statements

ELSEIF expression


statements

ELSE

statements

END

The statements are executed if the real part of the expression

has all non-zero elements. The ELSE and ELSEIF parts are optional.

Zero or more ELSEIF parts can be used as well as nested IF’s.

The expression is usually of the form expr rop expr where

rop is ==, <, >, <=, >=, or ~=.

Example

if I == J

A(I,J) = 2;

elseif abs(I-J) == 1

A(I,J) = -1;

else

A(I,J) = 0;

end

See also RELOP, ELSE, ELSEIF, END, FOR, WHILE, SWITCH.

Per i cicli di ripetizione possiamo utilizzare l’istruzione for e l’istruzione while

>> help for

FOR Repeat statements a specific number of times.

The general form of a FOR statement is:

FOR variable = expr, statement, ..., statement END

The columns of the expression are stored one at a time in

the variable and then the following statements, up to the

END, are executed. The expression is often of the form X:Y,

in which case its columns are simply scalars. Some examples

(assume N has already been assigned a value).

FOR I = 1:N,

FOR J = 1:N,

A(I,J) = 1/(I+J-1);

END

END

FOR S = 1.0: -0.1: 0.0, END steps S with increments of -0.1


FOR E = EYE(N), ... END sets E to the unit N-vectors.

Long loops are more memory efficient when the colon expression appears

in the FOR statement since the index vector is never created.

The BREAK statement can be used to terminate the loop prematurely.

See also IF, WHILE, SWITCH, BREAK, END.

>> help while

WHILE Repeat statements an indefinite number of times.

The general form of a WHILE statement is:

WHILE expression

statements

END

The statements are executed while the real part of the expression

has all non-zero elements. The expression is usually the result of

expr rop expr where rop is ==, <, >, <=, >=, or ~=.

The BREAK statement can be used to terminate the loop prematurely.

For example (assuming A already defined):

E = 0*A; F = E + eye(size(E)); N = 1;

while norm(E+F-E,1) > 0,

E = E + F;

F = A*F/N;

N = N + 1;

end

See also FOR, IF, SWITCH, BREAK, END.

Come esempio per l’analisi degli errori di arrotondamento scriviamo il seguente file ditipo script dal nome scriptrad.m . I commenti nel file vengono preceduti dal %

% Scr i p t f i l e% esempio e r r o r i d i arrotodamento% eseguiamo per n v o l t e l a rad i c e quadrata d i x e% per n v o l t e i l quadrato de l r i s u l t a t o . Prendiamo x=100

n = input ( ’ n = ’ )% l ’ i s t r u z i on e input v i s u a l i z z a l a s t r i n g a% e l e g g e i l dato d i input da t a s t i e r a .


xt = 100;x = 100;

for i =1:nx = sqrt ( x ) ;

end

for i =1:nx = x . ˆ 2 ;

end

%v i sua l i z z i amo l ’ e r rore r e l a t i v oe r r = abs (x−xt )/abs ( xt )

Adesso eseguiamo il file con diversi valori di n con la seguente istruzione

>> scriptrad

n = 10

n =

10

err =

6.3665e-014

>> scriptrad

n = 20

n =

20

err =

1.1555e-010

>> scriptrad

n = 30

n =


30

err =

2.2990e-007

>> scriptrad

n = 40

n =

40

err =

1.8975e-004

>> scriptrad

n = 60

n =

60

err =

9.9000e-001

>> xt

xt =

100

>> x

x =

1

Un’altro esempio e il calcolo delle seguenti funzioni equivalenti


x ∗ (√

x2 + 1− x)

x ∗√

x2 + 1− x2

x/(√

x2 + 1 + x)

con diversi valori di x.

Scriviamo il seguente file di tipo script dal nome scriptlimitefx.m

%%r i s u l t a t i d i v e r s i ca l co lando f un z i on i e q u i v a l e n t i%echo on% ca l c o l o d e l l e s e gu en t i f un z i on i e q u i v a l e n t i% 1) x ∗( s q r t ( xˆ2+1)−x% 2) x∗ s q r t ( xˆ2+1)−xˆ2% 3) x /( s q r t ( xˆ2+1)+x )pause% facciamo i l g r a f i c o d e l l e f un z i on i% ne l l ’ i n t e r v a l l o [10ˆ0 ,10ˆ9]% per f a r e i l g r a f i c o f i s s i amo a l c un i% punt i s u l l ’ asse d e l l e x e ca l co l i amo% l e t r e f un z i on i in q u e s t i punt ipause% con l ’ i s t r u z i on e log space ca l co l i amo% 500 punt i in s c a l a l o g a r i tm i c ax=logspace ( 0 , 9 , 500 ) ;pause% ca l c o l o l e f un z i on i% su t u t t o i l v e t t o r e d e l l e xr1=x . ∗ ( sqrt ( x.ˆ2+1)−x ) ;r2=x .∗ sqrt ( x.ˆ2+1)−x . ˆ 2 ;r3=x . / ( sqrt ( x.ˆ2+1)+x ) ;pause%fac c i o i l g r a f i c o in s c a l a l o g a r i tm i c a% su l l ’ asse xsemilogx (x , r1 , ’ r . ’ , x , r2 , ’b . ’ , x , r3 , ’ g− ’ ) ;

La funzione x=logspace(d1,d2,N) genera un vettore x con N elementi con spaziaturalogaritmica da 10d1 a 10d2.L’analoga funzione x=linspace(d1,d2,N) genera un vettore x con N elementi con spaziaturalineare da d1 a d2.L’istruzione semilogy genera grafici in scala logaritmica. L’istruzione plot genera graficiin scala lineare.Il punto prima di un operatore aritmetico permette di eseguire l’operazione sui vettorielemento per elemento.

Eseguiamo il file

>> scriptlimitefx


Algoritmi per la conversione di base

Vogliamo costruire una funzione Matlab per trasformare un numero x intero espresso inbase 10 in una base β > 1. Desideriamo quindi trovare le cifre α0, α1, . . . , αp tali che:

x = αpβp + αp−1β

p−1 + . . . + α0β0

Se consideriamo la seguente uguaglianza:

x = α0 + β(α1β0 + . . . + αpβ

p−1)

si noti che il valore di α0 si puo ottenere come resto della divisione di x per β. Infatti:

x = α0 + βx1

con x1 = α1β0 + . . . + αpβ

p−1, intero. Possiamo ripetere lo stesso procedimento su x1 pertrovare α1 ( α1 sara il resto della divisione di x1 per β).

La seguente e una funzione Matlab che esegue la conversione di base.

function a l f a = ConvInter i (x , beta )% dat i d i input% x = numero in t e ro in base 10% be ta = nuova base% da t i d i ou tput% a l f a = v e t t o r e contenente i c o e f f i c i e n t i in base N%

i =1;while ( x ˜= 0) & ( i <= 100)

a l f a ( i ) = rem(x , beta ) ;x = f ix ( x/beta ) ;

i=i +1;endi f ( x ˜=0) , error ( ’ r appr e s en taz ione troppo lunga ’ ) , end

Il simbolo % individua una linea di commento, cioe una linea in cui scriviamo informazioniche possono essere utili per capire il funzionamento del programma. Tali linee non vengonoeseguite dal Matlab, ma sono stampate se si esegue il comando help ConvInteri.

La funzione rem(x,beta) restituisce il valore del resto della divisione di x per beta, mentrela funzione fix(y) ha come valore la parte intera di y.

Il vettore alfa di output conterra i valori:

alfa (1) = α0, alfa (2) = α1, . . .

e importante inserire un controllo sul numero massimo di iterate perche un numero potrebbeavere una rappresentazione troppo lunga. Cio viene fatto mediante la condizione i <= 100nel ciclo while.

Esercizio 2.1.1 Sia x, 0 < x < 1, un numero in base 10, scrivere una funzione Matlab checalcoli i coefficienti α−i del numero in una base β > 1. Inserire un controllo sul numero


massimo di coefficienti consentito perche un numero potrebbe avere una rappresentazioneinfinita nella base β prescelta.

Suggerimento: essendo 0 < x < 1, si ha che:

x = α−1β−1 + α−2β

−2 + . . . + α−pβ−p

quindi, moltiplicando x per beta, si ha:

x ∗ β = α−1 + α−2β−1 + . . . + α−pβ

−p+1

e α−1 risulta essere la parte intera di x ∗ beta. Il procedimento si puo iterare applicandoloa:

x1 = α−2β−1 + . . . + α−pβ

−p+1.

Esercizio 2.1.2 Scrivere una funzione Matlab per convertire un numero reale in base 10in una base β > 1.

Algoritmo per determinare la base usata da un calcolatore

Per determinare la base usata da un calcolatore e utile considerare che fl(n + u) > nper 1 ≤ n < β e fl(β + u) = β possiamo quindi utilizzare la seguente procedura

function [ base ] = CalcolaBasebase = 1while ( base + eps ) > base

base = base + 1 ;end

La costante eps contiene il valore della precisione di macchina.

Algoritmi per il calcolo del valore di un polinomio.

p(x) = a1xN + a2x

N−1 + . . . + aN+1

Come possiamo valutarlo al variare di x?

Un algoritmo standard, che implementiamo nella function matlab pol1, e:

% Algoritmo per i l c a l c o l o de l v a l o r e d i% un pol inomio in un punto o in un insieme d i punt i% da t i d i input% p : v e t t o r e con i c o e f f i c i e n t i d e l pol inomio% p (1) x ˆ(n)+P(2) x ˆ(n−1) + . . . + p (n) x+p(n+1)% x : v e t t o r e in cu i v a l u t e r e i l pol inomio% da t i d i ou tput ;% y : v a l o r e de l pol inomio in xfunction y=pol1 (p , x )n = length (p)−1; % i s t r u z i on e l eng t h c a l c o l a l a lunghe z za de l v e t t o r e


y = p(n+1);for i=n:−1:1

y = y+p( i )∗x . ˆ ( n−i +1);end

Contiamo le operazioni aritmetiche:

addizioni : Nmoltiplicazioni : 1 + 2 + 3 + · · ·N = N(N+1)

2

Ogni termine ajxN−j e stato calcolato indipendentemente dagli altri termini.

Possiamo modificare l’algoritmo, che implementiamo nella function matlab pol2, calco-lando ricorsivamente xj = x ∗ xj−1.

L’algoritmo diventa:

% Algoritmo per i l c a l c o l o de l v a l o r e d i% un pol inomio in un punto o in un insieme d i punt i% da t i d i input% p : v e t t o r e con i c o e f f i c i e n t i d e l pol inomio% p (1) x ˆ(n)+P(2) x ˆ(n−1) + . . . + p (n) x+p(n+1)% x : v e t t o r e in cu i v a l u t e r e i l pol inomio% da t i d i ou tput ;% y : v a l o r e de l pol inomio in xfunction y=pol2 (p , x )n = length (p)−1; % i s t r u z i on e l eng t h c a l c o l a l a lunghe z za de l v e t t o r ey = p(n+1);xp = 1 ;for i=n:−1:1

xp = xp .∗ x ;y = y+p( i )∗xp ;

end

Le operazioni aritmetiche sono:addizioni: Nmoltiplicazioni : N + N − 1 = 2N − 1Il costo e molto inferiore rispetto al primo algoritmo. Esempio: N=20primo algoritmo: 210 moltiplicazionisecondo algoritmo: 39 moltiplicazioni

Un algoritmo ancora piu efficiente e la regola di Ruffini-Horner, che esegue le moltipli-cazioni in modo innestato

ESEMPI:N = 2 : p(x) = a2 + x(a1 + a0x)N = 3 : p(x) = a3 + x(a2 + x(a1 + xa0))N = 4 : p(x) = a4 + x(a3 + x(a2 + x(a1 + xa0)))


Il numero di operazioni e, rispettivamente, 2, 3 e 4 moltiplicazioni. Il secondo algoritmone richiedeva 3,5 e 7.

In generale:p(x) = aN+1 + x(aN + · · ·+ x(a2 + a1x)) · · ·)

L’algoritmo, che implementiamo nella function matlab horner, e :

% Algoritmo d i horner per i l c a l c o l o de l v a l o r e d i% un pol inomio in un punto o in un insieme d i punt i% da t i d i input% p : v e t t o r e con i c o e f f i c i e n t i d e l pol inomio% p (1) x ˆ(n)+P(2) x ˆ(n−1) + . . . + p (n) x+p(n+1)% x : v e t t o r e in cu i v a l u t e r e i l pol inomio% da t i d i ou tput ;% y : v a l o r e de l pol inomio in xfunction y=horner (p , x )m = length (p ) ;% i l grado de l pol inomio e ’ n=m−1y = p ( 1 ) ;for i =2:m

y = y .∗ x + p( i ) ;end

Le operazioni aritmetiche sono:addizioni: Nmoltiplicazioni : N

L’analisi degli algoritmi dal punto di vista del costo computazionale va sotto il nome diCOMPLESSITA DI CALCOLO.

L’algoritmo di Horner non solo e meno costoso ma anche piu stabile.

Scriviamo adesso il seguente file di tipo script dal nome scriptpol.m

% s c r i p t f i l e per con f ron tare va r i a l g o r i tm i per% i l c a l c o l o de i v a l o r i d e l pol inomio ( x−1)ˆ6%p = [ 1 −6 15 −20 15 −6 1 ] ;x = linspace ( 0 . 995 , 1 . 005 , 1000 ) ;

% Algoritmo 1w1 = pol1 (p , x ) ;

% Algoritmo2w2 = pol2 (p , x ) ;


% Algoritmo d i hornery = horner (p , x ) ;

% espre s s i onez = (x−1) .ˆ6 ;

plot (x , y , ’ . r ’ , x , z , ’−g ’ ,x ,w1 , ’ . b ’ , x , w2 , ’ . y ’ )

Eseguiamolo con l’istruzione

>> scriptpol

Approssimazione della derivata prima.

Sappiamo che per definizione la derivata prima di una funzione f(x) e data da

f ′(x) = limh→0

f(x + h)− f(x)

h

e utilizzando il polinomio di Taylor possiamo dire che:

f(x + h)− f(x)

h=

f(x) + hf ′(x) + h2

2 f ′′(ξ)− f(x)

h

con x ≤ ξ ≤ x + h.

Da qui deriva che

f(x + h)− f(x)

h= f ′(x) +

h

2f ′′(ξ)

La quantita τ(h) =h

2f ′′(ξ) si chiama errore di troncamento o errore di discretizzazione

e dipende da h. Possiamo dire che l’errore va a zero come O(h) cioe

∣

∣

∣

∣

τ(h)

h

∣

∣

∣

∣

e limitato per h→ 0.

Consideriamo il rapporto

f(x + h)− f(x− h)

2h


Utilizzando lo sviluppo in serie di Taylor abbiamo:

f(x + h)− f(x− h) =

= f(x) + hf ′(x) +h2

2f ′′(x) +

h3

6f ′′′(x) +

h4

24f (4)(ξ1)

−f(x) + hf ′(x)− h2

2f ′′(x) +

h3

6f ′′′(x)− h4

24f (4)(ξ2) =

= 2hf ′(x) +h3

3f ′′′(x) +

h4

24(f (4)(ξ1)− f (4)(ξ2))

In definitiva abbiamo

f(x + h)− f(x− h)

2h=

= f ′(x) +h2

6f ′′′(x) +

h3

24(f (4)(ξ1)− f (4)(ξ2)) =

= f ′(x) + O(h2)

Scriviamo un file di tipo script per eseguire lo studio dell’errore nel calcolo della derivataprima. Chiamiamo tale file scriptderivata.m:

%% stud i o d e l l ’ e r rore ne l c a l c o l o d e l l a de r i v a t a prima d i exp (1)% con l e d i f f e r e n z e in avan t i e l e d i f f e r e n z e c e n t r a l i%

h = logspace (−16 ,0 ,1000);

% approssimazione d e l l a de r i v a t a de l primo ordine

d f 1 = (exp(1+h)−exp ( 1 ) ) . / h ;

% approssimazione d e l l a de r i v a t a de l secondo ordine

d f 2 = (exp(1+h)−exp(1−h ) ) . / ( 2∗ h ) ;

% gra f i c o d e l l ’ e r rore

loglog (h , abs (exp(1)− d f 1 ) , ’ . b ’ ,h , abs (exp(1)− d f 2 ) , ’ . r ’ ) ;

L’istruzione loglog genera grafici in scala logaritmica su entrambi gli assi del sistemadi riferimento.

Eseguiamo il file

>> s c r i p t d e r i v a t a


Scriviamo adesso una function per approssimare la derivata prima dando in outputanche la stima dell’errore:

% func t i on per approssimare l a de r i v a t a prima% usando i l rapporto incrementa le% func t i on [ dfh , e r r r e l s ]= de r i v a t a f ( f , x , h , f l a g e s )% Dati d i input :% f : f unz ione% x : v a l o r e in cu i v a l u t a r e l a de r i v a t a% h : passo da u t i l i z z a r e% f l a g e s : opz iona le , = 1 : stampa i l g r a f i c o% Dati d i output% dfh : v a l o r e de l rapporto incrementa le% e r r r e l s : e r rore r e l a t i v o s t imatofunction [ dfh , e r r r e l s ]= d e r i v a t a f e r r A ( f , x , h , f l a g e s )

dfh = ( feval ( f , x+h)− feval ( f , x ) ) . / h ;dfh2 =( feval ( f , x+(h/2))− feval ( f , x ) ) . / ( h /2 ) ;e r r r e l s = abs ( dfh−dfh2 )./(1+ dfh2 ) ;

i f nargin == 4i f f l a g e s

loglog (h , e r r r e l s , ’ . ’ )end

end

questa funzione la possiamo eseguire tramite le seguenti istruzioni:

>> h = logspace (−56 ,0 ,500)>> [ dh , e r r ]= d e r i v a t a f e r r A ( ’ exp ’ ,1 , h , 1 ) ;

E necessario definire dei criteri di arresto che interrompano il succedersi delle iterazioninel momento in cui e stata raggiunta una stima buona dello zero.

Si puo arrestare il procedimento quando l’errore assoluto e minore di atol

err < atol

e/o l’errore relativo e minore di rtol

err

|dh| < rtol

con atol,rtol tolleranze specificate dall’utente. Il primo criterio e utile se |dh| ≈ 0.

Poiche non conosciamo il valore di dh, e preferibile usare il seguente criterio misto

err

atol + rtol|dh| < 1

Adesso, poiche abbiamo a disposizione la stima dell’errore, diamo in input le tolleranzeassolute e relative e calcoliamoci una approssimazione della derivata con un errore mistoche soddisfa le tolleranze di input. La function diventa:


% func t i on per approssimare l a de r i v a t a prima% usando i l rapporto incrementa le% func t i on [ dfh , e r r r e l s , h]= d e r i v a t a f t o l A ( f , x , a to l , r t o l , f l a g e s )% Dati d i input :% f : f unz ione% x : v e t t o r e in cu i v a l u t a r e l a de r i v a t a% ato l , r t o l : t o l l e r a n z e% f l a g e s : opz iona le , f l a g e s =1 stampiamo i l g r a f i c o% d e l l a de r i v a t a e i l g r a f i c o d e l l ’ e r rore s t imato% in ogni punto% Dati d i output% dfh : v a l o r e de l rapporto incrementa le% e r r r e l s = errore r e l a t i v o misto s t imato% h = passo u t i l i z z a t ofunction [ dfh2 , e r r r e l s , h]= d e r i v a t a f t o l ( f , x , ato l , r t o l , f l a g e s )h = 1 ;

dfh = ( feval ( f , x+h)− feval ( f , x ) ) . / h ;dfh2 =(feval ( f , x+(h/2))− feval ( f , x ) ) . / ( h /2 ) ;e r r r e l s = abs ( dfh−dfh2 ) . / ( a t o l+r t o l ∗abs ( dfh2 ) ) ;

while max( e r r r e l s ) > 1h=h/2 ;dfh = dfh2 ;dfh2 = ( feval ( f , x+(h/2))− feval ( f , x ) ) . / ( h /2 ) ;e r r r e l s = abs ( dfh−dfh2 ) . / ( a t o l+r t o l ∗abs ( dfh2 ) ) ;

end


figure (1 )plot (x , dfh2 , ’−. ’ )figure (2 )plot (x , e r r r e l s , ’ o ’ )

endend

Possiamo ottimizzare l’algoritmo utilizzando un passo diverso per i diversi valori di x,otteniamo la seguente function:

% func t i on per c a l c o l a r e i l rapporto incrementa le%%func t i on [ dfh2 , e r r r e l s , h]= d e r i v a t a f t o l h v ( f , x , a to l , r t o l , f l a g e s )%% DATI DI INPUT% f : funz ione d i cu i v a l u t a r e i l rap . inc .


% x : punto in cu i v a l u t a r e l a funz ione% a t o l : t o l l e r a n z a a s so l u t a% r t o l : t o l l e r a n z a r e l a t i v a% f l a g e s : opz i ona l e%%% DATI DI OUTPUT% dfh2 : rapporto incrementa le% e r r r e l s : e r rore r e l a t i v o s t imato% h : pa s s i u t i l i z z a t i%function [ dfh2 , e r r r e l s , h]= d e r i v a t a f t o l h v ( f , x , ato l , r t o l , f l a g e s )

h = ones ( s ize ( x ) ) ;index = 1 : length ( x ) ;

dfh2 ( index ) = ( feval ( f , x ( index)+h( index ))− feval ( f , x ( index ) ) ) . / h( index ) ;dfh = zeros ( s ize ( x ) ) ;e r r r e l s = zeros ( s ize ( x ) ) ;

while ˜isempty ( index )dfh ( index ) = dfh2 ( index ) ;dfh2 ( index ) =(feval ( f , x ( index )+(h( index )/2))− feval ( f , x ( index ) ) ) . / ( h( index ) / 2 )e r r r e l s ( index ) = abs ( dfh ( index)−dfh2 ( index ) ) . / ( a t o l+r t o l ∗abs ( dfh2 ( index ) ) ) ;index = find ( e r r r e l s > 1 ) ;h( index ) = h( index )/2 ;

end


figure (1 )plot (x , dfh2 , ’ ∗ ’ )figure (2 )plot (x , e r r r e l s , ’ ∗ ’ )

endend

close a l lx = linspace ( 1 , 2 , 2 0 ) ;[ dh , err , h]= d e r i v a t a f t o l ( @sin , x , 1 e−4,1e−5 ,1) ;hold on[ dhv , errv , hv]= d e r i v a t a f t o l h v ( ’ s i n ’ , x , 1 e−4,1e−5 ,1) ;

La @ prima del nome di una funzione indica l’handle, o identificativo della funzione,se si utilizza l’handle come dato di input non e’ necessario utilizzare la function feval pervalutare la funzione e si puo semplificare il codice nel seguente modo:


function [ dfh2 , e r r r e l s , h]= d e r i v a t a f t o l h v h and l e ( f , x , ato l , r t o l , f l a g e s )

h = ones ( s ize ( x ) ) ;index = 1 : length ( x ) ;

dfh2 ( index ) = ( f ( x ( index)+h( index))− f ( x ( index ) ) ) . / h( index ) ;dfh = zeros ( s ize ( x ) ) ;e r r r e l s = zeros ( s ize ( x ) ) ;

while ˜isempty ( index )dfh ( index ) = dfh2 ( index ) ;dfh2 ( index ) =( f ( x ( index )+(h( index )/2))− f ( x ( index ) ) ) . / ( h( index ) / 2 ) ;e r r r e l s ( index ) = abs ( dfh ( index)−dfh2 ( index ) ) . / ( a t o l+r t o l ∗abs ( dfh2 ( index ) ) ) ;index = find ( e r r r e l s > 1 ) ;h( index ) = h( index )/2 ;

end


figure (1 )plot (x , dfh2 , ’ ∗ ’ )figure (2 )plot (x , e r r r e l s , ’ ∗ ’ )

endend

Laboratorio Matematico-Informatico: Teoria · a.a 2010-2011, docente: Francesca Mazzia 1.1 Calcolo...

Documents

Transcript of Laboratorio Matematico-Informatico: Teoria · a.a 2010-2011, docente: Francesca Mazzia 1.1 Calcolo...