Complementi di Matematicae Calcolo Numerico

A.A. 2015-2016

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

avrà inevitabilmente a che fare con degli errori. Tali errori possono

essere di varia natura a seconda della causa che ne è 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 quantità 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 è impossibile rappresentare su di esso l’infinità

dei numeri reali, bisogna accontentarsi di un sottoinsieme di

numeri con un numero finito di cifre.

• Errori di troncamento dovuti al fatto che molti metodinumerici 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 è una quantità esatta ed x̃ è 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

quantità misurata. Diverso è un errore di un grammo nella misu-

ra del dosaggio di una confezione di preparato medicinale o nella

misura della quantità di cemento necessara a costruire un intero

palazzo.

Viceversa l’errore relativo fornisce indicazioni sulla qualità dell’ap-

prossimazione in quanto non è influenzato dall’ordine di grandezza

della quantità misurata anche se non da informazioni sulla effettiva

distanza tra la quantità esatta e quella approssimata.

La conoscenza di entrambi gli errori fornisce indicazioni precise

sulla bontà di una approssimazione.

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

nuti approssimando i valori del fattoriale dei primi 10 numeri na-

turali 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’infinità dei numeri reali si

utilizza un sottoinsieme di dimensione finita F ⊂ R detto insiemedei 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, β è detto base (tipicamente

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

0 β − 1) e p è detto esponente.Matlab utilizza la rappresentazione in doppia precisione, con la

quale è 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 più grande ed il più 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

Sarà Inf anche il risultato della valutazione di una espressione che

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

Matlab continuerà i calcoli.

Nota: Matlab riesce a rappresentare numeri più 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 mol-

to grossolano, infatti questo è solo un possibile formato del sistema

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

presentazione interna. Quest’ultima è data appunto dal numero

floating point “che più si avvicina” ad numero reale.

Ciò 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 è 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 dipenderà 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) è tale che

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

ǫM = β1−t è detto epsilon macchina, dove β è 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 è predefinita la variabile eps contenente il valore del-

l’epsilon macchina ǫM .

>> eps

ans =

2.2204e-16

8

ǫM non va confuso con il più piccolo numero rappresentabile in

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

può variare al più 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 può anche essere definito come il più 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/2altrimenti STOP

(2)

Soluzione: Come ogni linguaggio di programmazione Matlab pos-

siede istruzioni per effettuare cicli. Nel codice sotto riportato ve-

diamo un esempio di utilizzo dell’istruzione while, che permette

di ripetere ciclicamente un certo insieme di istruzioni fintanto che

una data condizione è 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 associatività. Alcune proprietà delle opera-

zioni 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 identità (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 realtà

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” è y2, mentre

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

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

In aritmetica esatta è 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)n

e 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 più semplice comando Matlab per disegnare un grafico è :

plot(x,y)

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

dimensione. Il comando plot(x,y) rappresenterà 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 funzio-

ni 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

Sarà 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 chesviluppando 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 di-

segni nei due casi il grafico del polinomio utilizzando il comando

plot.

Si osservi che, in aritmetica esatta, per l’identità (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 possibilità di scegliere il tipo di linea,

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

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 già 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 possibilità:

• possiamo disegnare il grafico di f e aggiungere successivamenteil 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 continuacon le seguenti modalità:

• 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 è lunga e/o complessa e la

funzione deve essere valutata in istanti successivi per diversi valori

delle variabili da cui dipende, è 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 è possibile de-

finire 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’

è 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: perchè il comando

>> plot(x,f)

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

È possibile definire funzioni che dipendono da più variabili o para-

metri

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

Attenzione: y = f(2, 10) è 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 − axcon 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 finestrautilizzando 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 kdisponendoli in 3 finestre distinte e corredando le figure di titoli,

ed etichette sugli assi.

21