Complementi di Matematicae Calcolo Numerico

A.A. 2014-2015

Laboratorio 2 - 9/3/2015

Quando si risolve un problema reale con l’ausilo di un computer si

avra inevitabilmente a che fare con degli errori. Tali errori possono

essere di varia natura a seconda della causa che ne e l’origine.

• Errori nei dati del problema dovuti ad esempio alla im-

precisione dello strumento di misura (microscopio,bilancia,ter-

mometro, barometro, buretta, densimetro etc.) utilizzato per

misurare una data quantita che si suppone nota per il problema

da risolvere.

• Errori di arrotondamento dovuti al fatto che per quan-

to grande e potente un computer ha pur sempre dimensione

finita e quindi e impossibile rappresentare su di esso l’infinita

dei numeri reali, bisogna accontentarsi di un sottoinsieme di

numeri con un numero finito di cifre.

• Errori di troncamento dovuti al fatto che molti metodi

numerici richiederebbero l’esecuzione di una sequenza infini-

ta di operazioni di calcolo che inevitabilmente devono essere

“troncati” dopo un numero finito di passaggi. Esempi: quan-

do si approssima con una somma parziale la somma di una

serie, quando si approssina un integrale con una formula di

quadratura, quando si risolve un sistema linare con un metodo

iterativo etc.

Definizione Se x e una quantita esatta ed x e una sua approssi-

mazione

|x− x|, |x− x||x|

sono detti rispettivamente errore assoluto ed errore relativo.

L’errore assoluto non tiene conto dell’ordine di grandezza della

quantita misurata. Diverso e un errore di un grammo nella misura

del dosaggio di una confezione di preparato medicinale o nella

misura della quantita di cemento necessara a costruire un intero

palazzo.

Viceversa l’errore relativo fornisce indicazioni sulla qualita dell’ap-

prossimazione in quanto non e influenzato dall’ordine di grandezza

della quantita misurata anche se non da informazioni sulla effettiva

distanza tra la quantita esatta e quella approssimata.

La conoscenza di entrambi gli errori fornisce indicazioni precise

sulla bonta di una approssimazione.

Esercizio 1: Si calcolino l’errore assoluto e quello relativo ottenu-

ti approssimando i valori del fattoriale dei primi 10 numeri natu-

rali con l’approssimazione di Stirling. ovvero per n = 1, . . . , 10

S(n) =√2 π n

(

ne

)n

. Visualizzare sulle colonne di una tabella: i

valori di n, n! S(n), |n!− S(n)|, |n!−S(n)||n!| .

2

Soluzione: Le seguenti righe di codice

n=[1:10]’;

Sn= sqrt(2*pi*n).*(n./exp(1)).^n;

fn=factorial(n);

Ea=abs(fn-Sn);

Er=Ea./abs(fn);

disp(’ n n! S(n) Ea Er’)

[n fn Sn Ea Er]

produrranno il risultato

n n! S(n) Ea Er

ans =

1.0000e+00 1.0000e+00 9.2214e-01 7.7863e-02 7.7863e-02

2.0000e+00 2.0000e+00 1.9190e+00 8.0996e-02 4.0498e-02

3.0000e+00 6.0000e+00 5.8362e+00 1.6379e-01 2.7298e-02

4.0000e+00 2.4000e+01 2.3506e+01 4.9382e-01 2.0576e-02

5.0000e+00 1.2000e+02 1.1802e+02 1.9808e+00 1.6507e-02

6.0000e+00 7.2000e+02 7.1008e+02 9.9218e+00 1.3780e-02

7.0000e+00 5.0400e+03 4.9804e+03 5.9604e+01 1.1826e-02

8.0000e+00 4.0320e+04 3.9902e+04 4.1760e+02 1.0357e-02

9.0000e+00 3.6288e+05 3.5954e+05 3.3431e+03 9.2128e-03

1.0000e+01 3.6288e+06 3.5987e+06 3.0104e+04 8.2960e-03

Si noti come l’errore assoluto aumenta all’aumentare di n mentre

l’errore relativo diminuisce.

3

Rappresentazione dei numeri reali:

Per rappresentare su un calcolatore l’infinita dei numeri reali si

utilizza un sottoinsieme di dimensione finita F ⊂ R detto insieme

dei numeri floating-point (virgola mobile) o numeri macchina, ge-

neralmente rappresentati in forma esponenziale normalizzata con

x = ±(0.a1 a2 . . . at) βp a1 6= 0 (1)

dove si intende

x = ±(a1β−1 + a2β

−2 + · · · + atβ−t) βp

e β, t, p sono numeri naturali positivi, β e detto base (tipicamente

β = 2), t e il numero massimo di cifre memorizzabili (comprese tra

0 β − 1) e p e detto esponente.

Matlab utilizza la rappresentazione in doppia precisione, con la

quale e possibile rappresentare ogni numero reale x che si trova

approssimativamente negli intervalli:

−1.88× 10308 ≤ x ≤ −2.23× 10−308

2.23× 10−308 ≤ x ≤ 1.80× 10308

Le variabili predefinte realmax, realmin di Matlab contengono

il piu grande ed il piu piccolo numero macchina positivo rappre-

sentabile in doppia precisione

>> realmax

ans= 1.7977e+308

>> realmin

ans = 2.2251e-308

4

Att.ne!

Nell’espressione 1.7977e+308, la e indica l’esponente nella poten-

za di 10 per cui si intende moltiplicato il numero 1.7977 e non ha

nulla a che vedere con il numero di Nepero e = 2, 7182818284590 . . . !

1.7977e+308 = 1.7977 ∗ 10308

Nota: Ad un numero maggiore di realmax Matlab associa il

valore Inf o “infinito”. Quindi ad esempio:

>> 2*realmax

ans =

Inf

Sara Inf anche il risultato della valutazione di una espressione che

contenga una divisione per 0 e senza generare un segnale di errore

Matlab continuera i calcoli.

Nota: Matlab riesce a rappresentare numeri piu piccoli di realmin

cambiando la rappresentazione ed utilizzando parte dei bit normal-

mente utilizzati per la mantissa per l’esponente. Questo comporta

una predita di cifre significative rispetto agli altri numeri macchi-

na se tale perdita supera una certo livello allora Matlab restituisce

semplicemente 0.

>> realmin/10

ans = 2.2251e-3090

>> realmin/1e+16

ans = 0

5

Att.ne! a non confondere la rappresentazione interna di un nu-

mero floating point in base 2 con il formato con cui tale numero

viene visualizzato a video da Matlab in rappresentazione decimale.

Numeri reali come 1/3 = 0.3 1/e = 0.367879441171442 . . .

vengono quindi approssimati con numeri macchina che hanno un

numero finito di cifre significative, digitando

>> 1/3

ans =

0.3333

>> 1/ exp(1)

ans =

0.3679

ne vediamo una rappresentazione.

Usare solo 4 cifre decimali per rappresentare i numeri sarebbe molto

grossolano, infatti questo e solo un possibile formato del sistema

per visualizzare un numero reale che non coincide con la sua rap-

presentazione interna. Quest’ultima e data appunto dal numero

floating point “che piu si avvicina” ad numero reale.

Cio si evidenzia cambiando il formato di visualizzazione.

6

Altri possibili formati sono:

>> format long

>> 1/3

ans= 0.3333333333333

>> 1/exp(1)

ans = 0.367879441171442

>> format short e

>> 1/3

ans = 3.3333e-01

>> 1/exp(1)

ans = 3.6788e-01

>> format long e

>> 1/3

ans = 3.333333333333333e-01

>> 1/exp(1)

ans = 3.678794411714423e-01

Esercizio 2 (Esistenza di un realmax e di un realmin.)

Partendo da 1 moltiplicare o dividere per 1e50

7

Aritmetica finita ed errori

Una conseguenza dell’approssimazione dei numeri reali con numeri

macchina e la costante presenza di errori di arrotondamento che con

l’esecuzione di successive operazioni aritmetiche si possono propa-

gare ed accumulare, per questo motivo in aritmetica finita il risul-

tato calcolato dipendera dalla sequenza con cui vengono eseguite

le varie operazioni.

L’errore relativo che si commette approssimando un numero reale

x con il numero floating point fl(x) e tale che

|fl(x)− x||x| ≤ C ǫM

ǫM = β1−t e detto epsilon macchina, dove β e la base e t

il numero di cifre significative dell’insieme dei numeri macchina

considerato. Nel caso della doppia precisione β = 2 e t = 53.

In Matlab e predefinita la variabile eps contenente il valore del-

l’epsilon macchina ǫM .

>> eps

ans =

2.2204e-16

8

ǫM non va confuso con il piu piccolo numero rappresentabile in

un sistema floating point, esso definisce invece una stima di quanto

puo variare al piu l’errore relativo quando si approssima un numero

reale con un numero macchina e quindi fornisce una stima del

numero di cifre significative corrette nella rappresentazione in un

determinato sistema floating point di un numero reale.

ǫM puo anche essere definito come il piu piccolo numero macchina

positivo tale che fl(1 + x) > 1.

Esercizio 3 Scrivere uno script-file che calcoli ǫM attraverso la

seguente ricorrenza:

e0 = 1/2, k = 0

se ek + 1 > 1 → ek+1 = ek/2

altrimenti STOP

(2)

Soluzione: Come ogni linguaggio di programmazione Matlab possiede

istruzioni per effettuare cicli. Nel codice sotto riportato vediamo un

esempio di utilizzo dell’istruzione while, che permette di ripetere

ciclicamente un certo insieme di istruzioni fintanto che una data

condizione e verificata.

epsmacchina.m

e=1/2; k=0;

while (e+1)>1

e=e/2; k=k+1;

end

e=2*e, k

9

Esempio di non associativita. Alcune proprieta delle oper-

azioni in aritmetica esatta possono non valere in aritmetica finita

in virgola mobile (floating point). Ad esempio:

>> a=1.0e+308;

>> b=1.1e+308;

>> c=-1.001e+308;

>> (a+b)+c

ans =

Inf

>> a+(b+c)

ans =

1.0990e+308

Esempio di cancellazione numerica. In aritmetica esatta,

usando la nota identita (a+b)(a−b) = a2−b2, si ottiene facilmente

√x2 + 1− x =

1√x2 + 1 + x

∀x ∈ R . (3)

Calcolando con Matlab:

>> x=77777777;

>> y1=sqrt(x^2+1)-x

y1 =

0

>> y2=1/(sqrt(x^2+1)+x)

y2 =

10

6.4286e-09

>> err=abs((y1-y2)/y2)

err =

1

In vista della uguaglianza (3), in aritmetica esatta i valori y1 e y2

dovrebbero essere uguali, e l’errore relativo err nullo. In realta

si osserva che i risultati ottenuti (y1 e y2) sono assai diversi. Il

risultato finale dipende fortemente da come viene effettivamente

calcolata la funzione (errore di arrotondamento, dovuto all’artit-

metica finita del calcolatore). Il risultato “corretto” e y2, mentre

il risultato dato da y1 e soggetto a un fenomeno di cancellazione.

Esercizio 4 (Esempio di errore dovuto all’aritmetica finita.)

In aritmetica esatta e ben noto il seguente limite

limn→∞

(

1 +1

n

)n

= e ,

dove e = 2, 718.. rappresenta il numero di Nepero.

Si scriva uno script file Matlab che una volta letto da tastiera un

valore per n calcoli en =(

1 + 1n

)ne l’errore assoluto ean = |en−e|

commesso. Per approssimare e si esegua ripetutamente tale script

file per n = 102, 104, 108, 1012, 1014, 1016. In aritmetica esatta, i

valori ean dovrebbero tendere a zero in quanto le en sopra calcolate

tendono a e. Cosa accade invece?

11

Grafici-2D

Il piu semplice comando Matlab per disegnare un grafico e :

plot(x,y)

dove x = (x1, . . . , xn) e y = (y1, . . . , yn) sono 2 vettori di ugual

dimensione. Il comando plot(x,y) rappresentera in una finestra

grafica una linea che collega i punti di coordinate (xi, yi), i =

1, . . . , n

Esempio

Disegnare il grafico della funzione f (x) = 2 sin(x) cos(x) + 2x

nell’intervallo [0, π/2]

>> x=linspace(0, pi/2,10);

>> y=2*sin(x).*cos(x)+2*x;

>> plot(x,y)

Basta un solo comando per valutare f in un vettore di punti gra-

zie alle operazioni componente per componente; inoltre le fun-

zioni matematiche (quali sin e cos) applicate a vettori eseguono

l’operazione componente per componente.

Se cambio il vettore delle ascisse x devo ricalcolare il vettore y

contenente i valori assunti dalla funzione f nei nuovi punti prima

di fare il grafico:

>> x=linspace(0,pi);

>> y=2*sin(x).*cos(x)+2*x;

>> plot(x,y)

12

Sara utile, per non dover manualmente ripetere la successione di

istruzioni ogni volta che si cambia un parametro, memorizzare il

lavoro in uno script-file.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.60

0.5

1

1.5

2

2.5

3

3.5

0 0.5 1 1.5 2 2.5 3 3.50

1

2

3

4

5

6

7

13

Esercizio 5

Scrivere uno script-file Matlab che, letti da tastiera gli estremi a

e b di un intervallo, disegni il grafico della funzione sin(x) per x ∈[a, b]. A tale scopo si generi un vettore z di 50 punti equispaziati

in [a, b], si valuti la funzione nei punti di z e si disegni il grafico

usando la funzione plot.

Esercizio 6 (Esempio grafico di cancellazione numerica.)

Per disegnare il grafico del polinomio p(x) = (x−1)6 si osservi che

sviluppando la potenza del binomio si ha la seguente uguaglianza

(x− 1)6 = x6 − 6x5 + 15x4 − 20x3 + 15x2 − 6x + 1 . (4)

Si utilizzino le due espressioni sopra citate per calcolare p(x) nei

punti x=[0.995:0.00001:1.005]; in due finestre distinte si dis-

egni nei due casi il grafico del polinomio utilizzando il comando

plot.

Si osservi che, in aritmetica esatta, per l’identita (4) i due grafici

disegnati dovrebbero essere identici. Invece, a causa degli effetti di

cancellazione nel calcolo del termine di destra il grafico presenta

forti oscillazioni artificiali.

14

Il comando plot prevede la possibilita di scegliere il tipo di linea,

il suo colore, ecc (vedi help plot). La sintassi generale e:

plot(x1, y1,’specifiche linea1’)

Per esempio per una linea rossa tratteggiata

plot(x,y,’r--’)

Il comando doc LineSpec mostra una lista di tutte le opzioni

disponibili per specificare lo stile della linea.

L’istruzione plot crea una nuova finestra grafica solo se non ci

sono finestre grafiche gia aperte, altrimenti utilizza l’ultima finestra

creata, e sovrascrive il nuovo grafico a quello creato in precedenza.

Pertanto se vogliamo visualizzare contemporaneamente i grafici di

due funzioni

f (x) = sin(x) + x, g(x) = x2 + cos(x) x ∈ [0, π]

abbiamo due possibilita:

• possiamo disegnare il grafico di f e aggiungere successivamente

il grafico di g specificando a Matlab di non cancellare il primo

grafico tramite il comando hold on

>> x=linspace(0,pi);

>> fx=sin(x)+x;

>> plot(x,fx,’r’)

>> hold on

>> gx=x.^2+cos(x);

>> plot(x,gx,’b’)

15

• possiamo disegnare i due grafici con un solo comando plot

>> x=linspace(0,pi);

>> fx=sin(x)+x;

>> gx=x.^2+cos(x);

>> plot(x,fx,’r’,x,gx,’b’)

0 0.5 1 1.5 2 2.5 3 3.50

1

2

3

4

5

6

7

8

9

Il comando figure genera una nuova finestra, la numera come

indicato e ci permette di disegnare i 2 grafici in 2 diverse finestre:

x=linspace(0,pi);

fx=sin(x)+x;

gx=x.^2 +cos(x);

figure(1)

plot(x,fx)

figure(2)

plot(x,gx)

16

I comandi title, xlabel,ylabel,legend ci permettono di com-

pletare il nostro grafico con un titolo, delle etichette sugli assi e una

legenda:

>> x=linspace(0,pi);

>> fx=sin(x)+x;

>> gx=x.^2 +cos(x);

>> plot(x,fx,’r’, x,gx,’b’)

>> title(’Grafici di funzioni’)

>> xlabel(’Asse x’)

>> ylabel(’Asse y’)

>> legend(’f’,’g’)

0 0.5 1 1.5 2 2.5 3 3.50

1

2

3

4

5

6

7

8

9Grafici di funzioni

Asse x

Ass

e y

fg

Il comando grid on inserisce una griglia sul grafico.

17

Esercizio 7

Disegnare il grafico delle seguenti funzioni:

• f (x) = 2 log(x+2)√x+1

x ∈ [1, 2] linea blu tratteggiata

• g(x) = x2+2x+1x2+1

x ∈ [0, 5] linea rossa continua

con le seguenti modalita:

• disegnare due grafici distinti nella stessa figura

• disporre i due grafici in due finestre distinte.

18

Funzioni Simboliche (inline)

Assegnata una funzione del tipo f (x) = (sin(x) + x)2 vogliamo

valutare i valori assunti da f per diversi valori di x.

Quando l’espressione della funzione e lunga e/o complessa e la fun-

zione deve essere valutata in istanti successivi per diversi valori delle

variabili da cui dipende, e utile poter definire la funzione una volta

per tutte senza dover riscrivere la sua espressione ogni volta che la

si vuole valutare in punti differenti. In Matlab e possibile definire

una funzione in modo simbolico. Il comando inline definisce una

funzione in linea, ovvero direttamente nello spazio di lavoro, senza

ricorrere a file esterni. Per esempio, definiamo:

>> f=inline(’(sin(x)+x).^2’,’x’)

dove il primo elemento in input

’(sin(x)+x).^2’

e una stringa che definisce l’espressione matematica che deve essere

trasformata in funzione, mentre i successivi elementi in input (nel

nostro caso solo ’x’) sono l’elenco ordinato delle variabili da cui

deve dipendere la funzione.

Attenzione: ricordarsi di utilizzare operazioni con i punti se si

vuole che la funzione operi sui vettori!

19

Ad una funzione cosı definita non sono associati dei valori numerici

(verificare con whos f). Per associarle valori numerici scriviamo,

per esempio

>> x=0:0.01:2*pi;

>> y=f(x);

Tali valori numerici vengono conservati nel vettore y (verificare con

whos y). Possiamo poi usarli, per esempio, per disegnare il grafico

di f con il comando

>> plot(x,y)

Attenzione: perche il comando

>> plot(x,f)

non funziona? Si potrebbe corregere in plot(x,f(x))?

E possibile definire funzioni che dipendono da piu variabili o parametri

>> f=inline(’2*x+a’,’a’,’x’);

Attenzione: y = f (2, 10) e diverso da z = f (10, 2).

20

Esercizio 8. Dopo averle definite con il comando inline, fare

un grafico delle seguenti funzioni:

a) f (x) = 2 log(x+2)√x+1

in [1, 2], linea nera tratteggiata;

b) f (x) = (x2 − 1)e1

x2−1 in [1, 10], linea rossa continua.

Esercizio 9. Dopo aver definito la funzione f (a, x) = e − ax

con il comando inline come funzione di due variabili, disegnare i

grafici di f in [−1, 1] per a = 1 e per a = −1 nella stessa finestra

utilizzando 2 colori diversi e corredandola di legenda;

Esercizio 10. Dopo aver definito la funzione f (k, x) = cos(kπx)

con il comando inline come funzione di due variabili k = 1, 2, 3

e x ∈ [−1, 1] disegnare i grafici di f per i differenti valori di k

disponendoli in 3 finestre distinte e corredando le figure di titoli,

ed etichette sugli assi.

21