Introduzionecalvino.polito.it/~puppo/didattica/dispensa.pdf · Introduzione Questa dispensa ... una...

Introduzione

Questa dispensa contiene alcuni testi d’esame per i corsi di calcolo numerico per gli studenti diTelecomunicazioni, Elettronica e Aerospaziale.

Nella prima parte, sono stati inseriti i testi degli esami, come sono stati forniti agli studenti il giornodell’esame.

Nella seconda parte, gli esercizi sono stati suddivisi per argomento, e risolti. Nella soluzione compareuna descrizione del ragionamento seguito, insieme ai comandi Matlab che permettono di ottenere irisultati richiesti. Inoltre, vengono forniti anche i risultati ottenuti.

Per utilizzare al meglio questa dispensa si consiglia di procedere nel modo seguente:

• Per prima cosa, studiate la teoria. Nel frattempo, consultate anche le dispense power pointche mostrano come vanno usate le functions di Matlab per realizzare gli schemi che vengonodiscussi nella parte di teoria. Alla fine di ogni argomento, fate gli esercizi che compaiono nel sitodel materiale didattico, in formato PDF. Notate che le nozioni di teoria vi servono per sapereche tipo di risultati dovreste ottenere e per rendervi conto di errori nell’impostazione dei vostriprogrammi.

• Finito lo studio preliminare, siete pronti per cominciare ad affrontare i problemi d’esame. Lavoratesu un tema d’esame per volta, senza consultare le soluzioni.

• Quando avete finito di risolvere un tema d’esame, consultate finalmente le soluzioni, per vederese i risultati combaciano, o se sono stati aggiunti dei commenti.

• Se vi siete incastrati irrimediabilmente su un esercizio e non riuscite a venirne a capo, consultatela soluzione di quell’esercizio, ma ricordatevi anche di ripassare l’argomento trattato, percheprobabilmente vi siete “un po’ affrettati” su uno dei punti precedenti.

Buon lavoro!

Gabriella Puppo

1

Temi d’esame

INGEGNERIA DELLE TELECOMUNICAZIONI - CALCOLO NUMERICO

Esame del 22 novembre 2002 - Prova di laboratorio

Turno 1

1. Scrivere una function che calcoli le matrici quadrate A[n] n × n , con n ≥ 3, tali che:

(A[n])i,j =

{

2 se j = i,−1 se j = i + 2 o j = i − 2,0 altrimenti.

1 ≤ i, j ≤ n

a) Stampare solo la matrice corrispondente a n = 7.

b) Per ogni valore di n assegnato, calcolare l’autovalore minimo λn di A[n], utilizzando lafunction eig per il calcolo degli autovalori.

c) Stampare un grafico dell’autovalore minimo λn in funzione di n, per n = 3, . . . , 100.

[Riportare i comandi Matlab usati per ottenere tutti i risultati.]

2. Considerare le due funzioni:

f(x) =

∫ x

0

(cos(5s)2) ds y(x) =1.1

3x.

Disegnare un grafico delle due funzioni sull’intervallo [0, 2], e dire se esistono punti di intersezionefra le due curve.


3. Calcolare il polinomio di grado 3, P 3(x) che interpola la funzione:

f(x) = e−x2

sui punti x1 = −2, x2 = 1, x3 = 2, e x4 = 3.2.

• Usare il polinomio P 3 per stimare la derivata di f(x) nel punto x = 3/2;

• Calcolare l’errore fra la stima di f ′ ed il risultato esatto.


2


Esame del 22 novembre 2002 - Prova di laboratorio

Turno 2

1. Considerare il problema differenziale:

{

y′ = f(y)y(0) = y0

,

dove:

y =

(

y1

y2

)

y0 =

(

10

)

e f(y) =

(

−2 −1a 3

)

∗

(

y1

y2

)

.

Sono assegnati tre valori per il parametro a, cioe a = 1; a = 2; a = 3.

• Disegnare su uno stesso grafico la soluzione y1, corrispondente ad ognuno dei valori di aassegnati: otterro quindi 3 funzioni sulla stessa figura.

• Trovare per quale fra i 3 valori di a y1(2) e piu vicino al numero −14.

2. Considerare la funzione:f(x) = sin(2x)e−x/2.

Interpolare f con una funzione spline cubica, Sn(x), sull’intervallo [−1, 3] con una griglia equi-spaziata di n + 1 punti. Studiare i due casi n = 5 e n = 10:

• disegnare il grafico di Sn e di f su una griglia grafica di 201 punti, sull’intervallo [−1, 3];

• calcolare l’errore fra Sn e f sulla griglia grafica in norma infinito.

3. Calcolare:∫ 7/3

0

√

|x − 1| dx, Risultato esatto: I =2

3+

2

3

√

(

4

3

)3

.

Usare la function quadl, con valore di tolleranza 10−8. Calcolare l’integrale in due modi:

• su tutto l’intervallo;

• spezzando l’intervallo in modo da avere una funzione integranda regolare.

Quale metodo da il risultato piu accurato? Quale metodo richiede meno operazioni? N.B. quadl

in output puo anche fornire il numero di operazioni effettuate.


3


Esame del 13 febbraio 2003 - Prova di laboratorio

Turno 1

1. E’ data la funzione:

y(x) =

{ √

(x) x ≤ 1(x − 2)2 x > 1

Disegnare un grafico della funzione assegnata sull’intervallo [0, 2], e calcolare l’integrale dellafunzione fra 0 e 2, usando le functions di Matlab.

2. Considerare la funzione:f(x) = | cos(2x)| ex.

Costruire una griglia {xi} equispaziata con 11 punti sull’intervallo [3, 5].

• Interpolare la funzione f sulla griglia {xi} con un polinomio che interpoli i dati esattamente.

• Interpolare la funzione f sulla griglia {xi} con un polinomio lineare a tratti.

Disegnare un grafico che contenga f ed i due interpolanti. Inoltre valutare l’errore in normainfinito per ognuno dei due problemi di interpolazione assegnati. Come si potrebbe fare perottenere degli interpolanti piu accurati per questa funzione?

3. Costruire una function che calcoli la matrice A n per n cosı definita:

A(i, j) =

{

2 se i ≤ j1 altrimenti.

con n assegnato in input.

Per n = 100, risolvere il sistema lineare:

A2x = b, con b(i) ≡ 1, i = 1, . . . , 100.

Cercate, se riuscite, di risolvere il sistema lineare in maniera efficiente, senza calcolare A2, e conuna sola fattorizzazione LU .


4

CORSI DI LAUREA TLC, ELN, AER - CALCOLO NUMERICO

Esame dell’ 1 settembre 2003 - Prova di laboratorio

1. Costruire la matrice quadrata 10 per 10 definita dalla relazione:

(A)i,j ={

0 se j > i,1 altrimenti

1 ≤ i, j ≤ 10

e stamparla. Costruire inoltre il vettore colonna x, definito da xi = 1 per i = 1, · · · , 10.

Stampare la matrice i cui elementi sono (A10)i,j e la matrice i cui elementi sono (Ai,j)10. Quali

comandi permettono di ottenere queste due matrici?

Per n = 1, · · · , 10, calcolare il vettore bn = An ∗ x. Disegnare un grafico di ||bn|| in funzione di n.


2. Considerare le due funzioni:

f(x) =

∫ x

0

(10− s)(log s)2 ds g(x) = x2.

Disegnare un grafico delle due funzioni sull’intervallo [1, 11], e stimare la posizione dei punti diintersezione fra le due curve, nell’intervallo [1, 11], usando le informazioni fornite dal grafico.


3. Considerare il problema di Cauchy:

{

y′ = t cos(y2 + t)y(1) = 0

Risolvere il problema differenziale fino a t = 10, usando la ode-function ode45 di Matlab. De-terminare il passo di integrazione minimo ed il passo di integrazione massimo utilizzati nel corsodell’integrazione.


5

CORSO DI LAUREA AER - CALCOLO NUMERICO

Esame dell’ 8 settembre 2003 - Prova di laboratorio

1. Costruire la matrice 20 per 10 definita dalla relazione:

(A)i,j ={

2 se i = j,1 altrimenti

1 ≤ i ≤ 20; 1 ≤ j ≤ 10.

ed il vettore b ∈ R20 tale che bi = (−1)i+1. Stampare A e b. Considerare ora il sistemasovradeterminato

Ax = b

e risolverlo nel senso dei minimi quadrati nei due modi seguenti.

(a) Risolvere il problema con il metodo QR. Sia x1 la soluzione ottenuta in questo modo, e siar1 = ||Ax1 − b|| il residuo corrispondente.

(b) Risolvere il problema con il metodo delle equazioni normali. Siano rispettivamente x2 e r2

la soluzione ed il corrispondente residuo ottenuti con questa tecnica.

Calcolare infine le quantita ∆x = ||x1 − x2|| e ∆r = |r1 − r2|. Stampare e commentare i risultatiottenuti.


2. Considerare l’integrale:

I =

∫ 2

−1

xex dx,

il cui valore esatto e I = e2 + 2e−1. Per k = 1, · · · , 10, stimare il valore dell’integrale, usando lafunction quad di Matlab, con tolleranza tol = 10−k. Stampare l’errore fra integrale stimato evalore esatto, per ogni valore di k assegnato. Per quale valore di k si ottiene il risultato migliore?Quanto vale l’errore minimo?


3. Considerare il problema di Cauchy:

{

y′ = t cos(y2 + t)y(1) = 0

Risolvere il problema differenziale fino a t = 10, usando la ode-function ode45 di Matlab.Determinare i valori massimi e minimi della soluzione nell’intervallo [1, 10].


6

Capitolo 1

Algebra lineare

1.1 Esercizio 1, 22/11/02 Turno 1

Scrivere una function che calcoli le matrici quadrate A[n] n × n, con n ≥ 3, tali che:

(A[n])i,j =

{

2 se j = i,−1 se j = i + 2 o j = i − 2,0 altrimenti.

1 ≤ i, j ≤ n

a) Stampare solo la matrice corrispondente a n = 7.

b) Per ogni valore di n assegnato, calcolare l’autovalore minimo λn di A[n], utilizzando la functioneig per il calcolo degli autovalori.

c) Stampare un grafico dell’autovalore minimo λn in funzione di n, per n = 3, . . . , 100.

Soluzione

Questo e il testo della function con la quale creo le matrici richieste dall’esercizio:

function a=tridiag(n)

if n<3; return; end;

a=2*eye(n);

for k=3:n

a(k-2,k) = -1;

a(k,k-2) = -1;

end

La matrice corrispondente a n = 7 e data da:

>>tridiag(7)

ans =

2 0 -1 0 0 0 0

0 2 0 -1 0 0 0

-1 0 2 0 -1 0 0

0 -1 0 2 0 -1 0

0 0 -1 0 2 0 -1

0 0 0 -1 0 2 0

0 0 0 0 -1 0 2

Per costruire il grafico richiesto, uso i seguenti comandi:

nmax = 100;

for n=3:nmax

7

a = tridiag(n);

v = eig(a);

vmin(n) = min(v);

end

plot(3:nmax,vmin(3:nmax)

Il grafico ottenuto si trova in Fig. 1.1.

0 10 20 30 40 50 60 70 80 90 1000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 1.1: Esercizio 1.1: Autovalore minimo di A[n] in funzione di n.

1.2 Esercizio 3, 13/2/2003

Costruire una function che calcoli la matrice A n per n cosı definita:

A(i, j) =

{

2 se i ≤ j1 altrimenti.

con n assegnato in input.Per n = 100, risolvere il sistema lineare:

A2x = b, con b(i) = 1, i = 1, . . . , 100.

Cercate, se riuscite, di risolvere il sistema lineare in maniera efficiente, senza calcolare A2, e con unasola fattorizzazione LU .Soluzione

La function seguente calcola la matrice richiesta:

function a=mat(n)

a=ones(n);

for i=1:n

for j=i:n

a(i,j) = 2;

end

end

Per calcolare il sistema in modo efficiente, osservo che: A2x = b ⇒ A(Ax) = b, quindi invece del sistemaassegnato posso risolvere i due sistemi:

Ax1 = b quindi Ax = x1.

8

Poiche entrambi i sistemi hanno la stessa matrice dei coefficienti, posso utilizzare due volte la stessafattorizzazione LU:

n=100;

a=mat(n);

b=ones(n,1);

[l,u]=lu(a);

x1=u\(l\b);

x=u\(l\x1);

Attenzione:l’uso delle parentesi qui e essenziale, perche

u\(l\b)

e diverso da

u\l\b

.

1.3 Esercizio 1, 1/9/03

Costruire la matrice quadrata 10 per 10 definita dalla relazione:

(A)i,j ={

0 se j > i,1 altrimenti

1 ≤ i, j ≤ 10

e stamparla. Costruire inoltre il vettore colonna x, definito da xi = 1 per i = 1, · · · , 10.Stampare la matrice i cui elementi sono (A10)i,j e la matrice i cui elementi sono (Ai,j)

10. Qualicomandi permettono di ottenere queste due matrici?

Per n = 1, · · · , 10, calcolare il vettore bn = An ∗ x. Disegnare un grafico di ||bn|| in funzione di n.Soluzione

Questo e il segmento di programma che costruisce la matrice richiesta ed il vettore x:

% Costruisce la matrice

k=10; a=zeros(k);

for i=1:k-1

for j=i+1:k

a(i,j) = 1;

end

end

a

% Costruisce il vettore x

x=ones(k,1);

La matrice A e data da:

a =

0 1 1 1 1 1 1 1 1 1

0 0 1 1 1 1 1 1 1 1

0 0 0 1 1 1 1 1 1 1

0 0 0 0 1 1 1 1 1 1

0 0 0 0 0 1 1 1 1 1

0 0 0 0 0 0 1 1 1 1

0 0 0 0 0 0 0 1 1 1

0 0 0 0 0 0 0 0 1 1

0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0

9

La matrice i cui elementi sono (A10)i,j si ottiene con il comando

a^10

ed e data da:

>> a^10

ans =

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

La matrice i cui elementi sono (Ai,j)10 si ottiene con il comando

a.^10

e coincide con la matrice A iniziale.L’ultima parte dell’esercizio si risolve nel modo seguente:

for n=1:k

b(n)=norm((a^n)*x);

end

plot(b)

Il grafico di ||bn|| compare in figura 1.2.

1 2 3 4 5 6 7 8 9 100

50

100

150

Figura 1.2: Esercizio 1.3: ||bn|| in funzione di n.

1.4 Esercizio 1, 8/9/03

Costruire la matrice 20 per 10 definita dalla relazione:

(A)i,j ={

2 se i = j,1 altrimenti

1 ≤ i ≤ 20; 1 ≤ j ≤ 10.

10

ed il vettore b ∈ R20 tale che bi = (−1)i+1. Stampare A e b. Considerare ora il sistema sovradeterminato

Ax = b

e risolverlo nel senso dei minimi quadrati nei due modi seguenti.

1. Risolvere il problema con il metodo QR. Sia x1 la soluzione ottenuta in questo modo, e siar1 = ||Ax1 − b|| il residuo corrispondente.

2. Risolvere il problema con il metodo delle equazioni normali. Siano rispettivamente x2 e r2 lasoluzione ed il corrispondente residuo ottenuti con questa tecnica.

Calcolare infine le quantita ∆x = ||x1 − x2|| e ∆r = |r1 − r2|. Stampare e commentare i risultatiottenuti.Soluzione

Cominciamo a costruire le matrici richieste:

a=ones(20,10);

for i=1:10

a(i,i) = 2;

end

b=ones(20,1);

for i=2:2:20

b(i)=-1;

end

Tralascio la stampa delle matrici ottenute, per brevita.Risolvo ora nel senso dei minimi quadrati. Ricordiamo innanzitutto che l’operatore backslash di

Matlab, nel caso dei sistemi sovradeterminati, usa il metodo QR. Quindi per la parte a) dell’eserciziosono sufficienti le seguenti istruzioni:

% Risolve il problema con il metodo QR

x1=a\b;

r1=norm(a*x1-b);

Per risolvere il problema con il metodo delle equazioni normali e invece necessario formare prima lamatrice AT A e poi risolvere il sistema quadrato AT Ax = AT b:

% Risolve il problema con il metodo delle equazioni normali

aa=a’*a;

bb=a’*b;

x2 = aa\bb;

r2=norm(a*x2-b);

Infine, stampa i due risultati richiesti:

dx=norm(x1-x2)

dr=abs(r1-r2)

I valori ottenuti sono:

dx =

2.935275148808155e-15

dr =

4.440892098500626e-16

I due valori ∆x e ∆r sono molto piccoli, perche il metodo QR ed il metodo delle equazioni normalisono due modi equivalenti per risolvere il sistema lineare nel senso dei minimi quadrati. La soluzioneed il residuo ottenuti sono quindi approssimativamente gli stessi, se l’errore di tipo floating point epiccolo. In questo caso, l’errore numerico e piccolo, perche la matrice A e ben condizionata: infatticond(a) = 14.8660. Se la matrice A fosse mal condizionata, le due soluzioni x1 ed x2 ed i relativiresidui potrebbero essere anche molto diversi fra loro.

11

Capitolo 2

Approssimazione di funzioni


Calcolare il polinomio di grado 3, P 3(x) che interpola la funzione:

f(x) = e−x2

nei punti x1 = −2, x2 = 1, x3 = 2, e x4 = 3.2.

• Usare il polinomio P 3 per stimare la derivata di f(x) nel punto x = 3/2;

• Calcolare l’errore fra la stima di f ′ ed il risultato esatto.

Soluzione

Prepara i dati da interpolare e calcola il polinomio di interpolazione:

xi=[-2, 1, 2, 3.2];

f=inline(’exp(-x.^2)’);

fi=f(xi);

p=polyfit(xi,fi,3);

Per stimare la derivata di f , calcola la derivata di P 3, nel punto x = 3/2:

xx=3/2;

p1=3*p(1)*xx^2+2*p(2)*xx+p(3);

Infine, calcola la derivata di f in x = 3/2 e calcola l’errore:

f1=-2*xx*f(xx);

err=abs(p1-f1);

Si ottengono i risultati seguenti:

>> p1

p1 =

-0.3625

>> f1

f1 =

-0.3162

>> err

err =

0.0463

12


Considerare la funzione:f(x) = sin(2x)e−x/2.

Interpolare f con una funzione spline cubica, Sn(x), sull’intervallo [−1, 3] con una griglia equispaziatadi n + 1 punti. Studiare i due casi n = 5 e n = 10:

• disegnare il grafico di Sn e di f su una griglia grafica di 201 punti, sull’intervallo [−1, 3];

• calcolare l’errore fra Sn e f sulla griglia grafica in norma infinito.

Soluzione

Lo script che ho scritto per risolvere questo problema e:

f=inline(’sin(2*x) .*exp(-x/2)’);

xx=linspace(-1,3,201);

fxx=f(xx);

k=0;

for n=[5 10]

xi=linspace(-1,3,n+1);

fi=f(xi);

s=interp1(xi,fi,xx,’spline’);

fxx=f(xx);

k=k+1;

err(k)=norm(s-fxx,inf);

figure

plot(xx,fxx)

hold on

plot(xx,s,’:’)

plot(xi,fi,’*’)

end

In questo modo, si creano due figure diverse, che ho riportato in Fig. 2.1. Per n = 5, la funzione ela sua interpolante spline non coincidono perfettamente, mentre per n = 10 sono indistinguibili. Glierrori in norma infinito sono dati da:

>> err

err =

0.1254 0.0035

2.3 Esercizio 2, 13/2/2003

Considerare la funzione:f(x) = | cos(2x)| ex.

Costruire una griglia {xi} equispaziata con 11 punti sull’intervallo [3, 5].

• Interpolare la funzione f sulla griglia {xi} con un polinomio che interpoli i dati esattamente.

• Interpolare la funzione f sulla griglia {xi} con un polinomio lineare a tratti.

Disegnare un grafico che contenga f ed i due interpolanti. Inoltre valutare l’errore in norma infinitoper ognuno dei due problemi di interpolazione assegnati. Come si potrebbe fare per ottenere degliinterpolanti piu accurati per questa funzione?Soluzione

Calcola i dati per interpolare la funzione f :

13

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

−1.5

−1

−0.5

0

0.5

1

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

−1.5

−1

−0.5

0

0.5

1

Figura 2.1: Esercizio 2.2. Grafico di f (linea continua) e dell’interpolante spline (linea a puntini). I marker(*) indicano i punti di interpolazione. A sinistra n = 5, a destra n = 10.

f=inline(’abs(cos(2*x)).*exp(x)’);

n=10;

xi =linspace(3,5,n+1);

fi =f(xi);

Per interpolare i dati esattamente, devo usare un polinomio di grado n = 10. Valuto entrambi gliinterpolanti su una griglia fitta (in questo caso 201 punti), per poter disegnare il grafico e valutare lanorma dell’errore:

p=polyfit(xi,fi,n);

xx=linspace(3,5,201);

fxx=f(xx);

plot(xx,fxx)

hold

pxx=polyval(p,xx);

plot(xx,pxx,’:’)

plot(xi,fi,’*’)

s=interp1(xi,fi,xx,’linear’);

plot(xx,s,’--’)

err_pol=norm(fxx-pxx,inf)

err_lin=norm(fxx-s,inf)

Il grafico richiesto compare in Fig. 2.2. L’errore per il polinomio di interpolazione e:

err_pol =

32.38474255677069

mentre l’errore per l’interpolazione lineare a tratti e:

err_lin =

8.78831066653374

Notare che in questo caso l’interpolante lineare a tratti fornisce un risultato migliore. Per ottenereun buon risultato con il polinomio di interpolazione, sarebbe necessario dividere l’intervallo di inter-polazione in due tratti, in ciascuno dei quali la funzione f fosse regolare e risolvere due problemi diinterpolazione distinti su ognuno dei due intervalli.

14

3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5−20

0

20

40

60

80

100

120

140

Figura 2.2: Esercizio 2.3. Grafico di f(x) (curva continua), polinomio di interpolazione (curva a punti) einterpolante lineare a tratti (curva tratteggiata).

15

Capitolo 3

Integrazione numerica


Considerare le due funzioni:

f(x) =

∫ x

0

(cos(5s)2) ds y(x) =1.1

3x.

Disegnare un grafico delle due funzioni sull’intervallo [0, 2], in modo da determinare se esistono puntidi intersezione fra le due curve.Soluzione

Per prima cosa costruisco le funzioni richieste, e la griglia sulla quale calcolero il grafico:

f=inline(’cos((5*x).^2)’);

y=inline(’1.1*x/3’);


Per rendere piu veloce il calcolo della funzione integrale, uso le formule:

F (xk) = F (xk−1) +

∫ xk

xk−1

f(s) ds F (0) = 0.

Con questo accorgimento, uso i comandi seguenti:

ff(1) = 0;

for k=2:length(xx)

ff(k)=ff(k-1)+quad(f,xx(k-1),xx(k));

end

plot(xx,ff)

hold

plot(xx,y(xx),’--’)

Il grafico delle due funzioni si trova in Fig. 3.1. Dal grafico, si vede facilmente che ci sono due punti diintersezione: uno in x = 0, l’altro nell’intervallo (0.3, 0.4).

Diverse persone hanno interpretato questo esercizio in modo diverso, come se la funzione integrandafosse:

f(x) =

∫ x

0

(cos(5s))2 ds y(x) =1.1

3x.

In questo caso, la struttura della strategia di soluzione resta la stessa, ma i risultati sono assai diversi:

f=inline(’cos(5*x).^2’);

y=inline(’1.1*x/3’);


16

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Figura 3.1: Esercizio 3.1. Grafico delle due funzioni.

ff(1) = 0;

for k=2:length(xx)

ff(k)=ff(k-1)+quad(f,xx(k-1),xx(k));

end

plot(xx,ff)

hold

plot(xx,y(xx),’--’)

Il grafico delle due funzioni si trova sulla sinistra di Fig. 3.2. Sulla destra della figura, un ingrandimentodel grafico mostra che c’e un solo punto di intersezione, che e x = 0.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.2

0.4

0.6

0.8

1

1.2

0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65

0.05

0.1

0.15

0.2

0.25

0.3

Figura 3.2: Esercizio 3.1. A sinistra: grafico delle due funzioni. A destra: particolare.

17


Calcolare:∫ 7/3

0

|x − 1| dx, Risultato esatto: I =2

3+

2

3

√

(

4

3

)3

.

Usare la function quadl, con valore di tolleranza 10−9. Calcolare l’integrale in due modi:

• su tutto l’intervallo;

• spezzando l’intervallo in modo da avere una funzione integranda regolare.

Quale metodo da il risultato piu accurato? Quale metodo richiede meno operazioni? N.B. quadl in

output puo anche fornire il numero di operazioni effettuate.Soluzione

Calcolo l’integrale su tutto l’intervallo, con la tolleranza richiesta, contando anche il numero di opera-zioni:

f=inline(’sqrt(abs(x-1))’);

tol=1e-8;

[int,c]=quadl(f,0,7/3,tol);

Calcola lo stesso integrale, spezzando il dominio di integrazione sui due intervalli: [0, 1] e [1, 7/3]:

[int1,c1]=quadl(’sqrt(1-x)’,0,1,tol);

[int2,c2]=quadl(’sqrt(x-1)’,1,7/3,tol);

ctot=c1+c2;

int_tot=int1+int2;

L’errore ed il numero di operazioni sull’intervallo intero sono:

exa=2/3+2/3*sqrt( (4/3)^3)

err1=exa-int

c

exa =

1.69306714522600

err1 =

-3.018807426258263e-11

c =

378

L’errore ed il numero di operazioni sull’intervallo spezzato sono:

err2=exa-int_tot

ctot

err2 =

2.630707207629257e-09

ctot =

366

In questo caso, sorprendentemente, il risultato piu accurato e stato ottenuto senza spezzare l’integrale.

3.3 Esercizio 1, 13/2/2003

E’ data la funzione:

y(x) =

{ √

(x) x ≤ 1(x − 2)2 x > 1

18

Disegnare un grafico della funzione assegnata sull’intervallo [0, 2], e calcolare l’integrale della funzionefra 0 e 2, usando le functions di Matlab.Soluzione

Scrivo una function che valuti la funzione y(x). Per far questo ho bisogno di un file .m. Per esempio:

function f=funz7(x)

n=length(x);

for i=1:n

if x(i)<1

f(i) = sqrt(x(i));

else

f(i) = (x(i)-2)^2;

end

end

Notare che devo scrivere la funzione componente per componente: in questo modo, posso assegnare unvettore di ascisse in input, ed ottenere un vettore con i valori di f in output.

Per disegnare il grafico e calcolare l’integrale, uso i seguenti comandi:


fxx=funz7(xx);

plot(xx,fxx)

int=quad(@funz7,0,1) + quad(@funz7,1,2)

Ho spezzato l’integrale nella somma di due contributi, per non avere il punto angoloso che si trova inx = 1 all’interno dell’intervallo di integrazione. Il valore dell’integrale e

int =

1.0000

e corrisponde al valore esatto, che in questo caso e facile da calcolare. La figura con il grafico di y(x)si trova in Fig. 3.3.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 3.3: Esercizio 3.3. Grafico di y(x).

3.4 Esercizio 2, 1/9/03

Considerare le due funzioni:

f(x) =

∫ x

0

(10 − s)(log s)2 ds g(x) = x2.

19

Disegnare un grafico delle due funzioni sull’intervallo [1, 11], e stimare la posizione dei punti di interse-zione fra le due curve, nell’intervallo [1, 11], usando le informazioni fornite dal grafico.Soluzione

Costruisco una griglia di 101 punti sull’intervallo [1, 11], in modo da avere punti di griglia abbastanzafitti. Calcolo poi il valore di f sui punti di griglia usando l’istruzione quad. Un modo efficiente perottenere questo risultato e il seguente:

a=1; b=11;

n=100;

x=linspace(a,b,n+1);

f1=inline(’(log(x)).^2 .*(10-x)’)

% Calcola l’integrale fra x(i-1) e x(i) di f1

int(1) = 0;

for i=2:n+1

int(i) = quad(f1,x(i-1),x(i));

end

f(1) = 0;

for i=2:n+1

f(i) = f(i-1)+int(i);

end

plot(x,f)

Il grafico della seconda funzione puo essere ottenuto nel modo seguente:

% Disegna il grafico della seconda funzione

f2 = inline(’x.^2’);

f2x=f2(x);

hold

plot(x,f2x,’--’)

grid on

Il grafico delle due funzioni compare in 3.4. Si nota immediatamente che ci sono due punti di intersezionenell’intervallo [1, 11]. Facendo uno zoom sui due punti di intersezione, si conclude che c’e una soluzionenell’intervallo [3.4, 3.5], mentre la seconda soluzione si trova in [8.2, 8.3].

1 2 3 4 5 6 7 8 9 10 110

20

40

60

80

100

120

140

Figura 3.4: Esercizio 3.4: Grafico della funzione f (linea continua), e della funzione g (linea tratteggiata).

20

3.5 Esercizio 2, 8/9/03

Considerare l’integrale:

I =

∫ 2

−1

xex dx,

il cui valore esatto e I = e2 +2e−1. Per k = 1, · · · , 10, stimare il valore dell’integrale, usando la function

quad di Matlab, con tolleranza tol = 10−k. Stampare l’errore fra integrale stimato e valore esatto, perogni valore di k assegnato. Per quale valore di k si ottiene il risultato migliore? Quanto vale l’erroreminimo?Soluzione

Per calcolare gli integrali richiesti, posso usare il programma seguente,

% Calcola il valore esatto dell’integrale

iexa=exp(2)+2*exp(-1);

% Calcola ora l’integrale approssimato

f=inline(’x.*exp(x)’);

format long

for k=1:10

tol=10^(-k);

int=quad(f,-1,2,tol);

err(k) = abs(iexa-int);

end

err

Si ottengono questi risultati:

err =

Columns 1 through 4

0.01342731722420 0.00026974557270 0.00004443962649 0.00000186920833

Columns 5 through 8

0.00000026839315 0.00000001180457 0.00000000046654 0.00000000003682

Columns 9 through 10

0.00000000000257 0.00000000000013

dai quali si vede facilmente che l’errore minimo e err = 0.00000000000013, che si ottiene per k = 10.I risultati cambiano drasticamente, se considero tolleranze piu piccole. Per esempio, per k =

1, · · · , 20, si osserva che l’errore comincia a crescere per valori della tolleranza troppo piccoli. Ri-portando i valori dell’errore in funzione di k in scala logaritmica, si ottengono i risultati in Fig. 3.5.In questo grafico sono evidenti due comportamenti in competizione: dapprima, la diminuizione dellatolleranza determina una migliore precisione della formula di quadratura e l’errore diminuisce. Succes-sivamente pero diventa dominante l’errore di tipo floating point, che aumenta velocemente perche ilnumero dei conti richiesti cresce ed accumula in maniera esplosiva l’errore di origine floating point.

21

0 2 4 6 8 10 12 14 16 18 20−35

−30

−25

−20

−15

−10

−5

0

5

Figura 3.5: Esercizio 3.5: Grafico in scala logaritmica dell’errore nel calcolo dell’integrale in funzione di k.

22

Capitolo 4

Equazioni differenziali


Considerare il problema differenziale:{

y′ = f(y)y(0) = y0

,

dove:

y =

(

y1

y2

)

y0 =

(

10

)

e f(y) =

(

−2 −1a 3

)

∗

(

y1

y2

)

.

Sono assegnati tre valori per il parametro a, cioe a = 1; a = 2; a = 3.

• Disegnare su uno stesso grafico la soluzione y1, corrispondente ad ognuno dei valori di a assegnati:otterro quindi 3 funzioni sulla stessa figura.

• Trovare per quale fra i 3 valori di a y1(2) e piu vicino a −14.

Soluzione

Innanzitutto devo scrivere la function per definire il sistema differenziale. Inizialmente, assegno ilvalore a = 1:

function f=funz4(t,y)

a=1;

f=[-2 -1; a 3]*y;

Trovo la prima soluzione:

[t,y]=ode45(@funz4,[0 2],[1 0]);

Uso la function ode45, perche non ho bisogno di molta precisione. Se il sistema si rivelasse stiff,passerei ad usare la ode15s. Posso ora disegnare il grafico della prima componente della soluzione,

plot(t,y(:,1),’:’)

Ottengo la curva a puntini di Fig. 4.1. Sto usando il tratteggio, perche vorrei un grafico che si possastampare in bianco e nero. Per calcolare le altre soluzioni richieste, devo tornare al file funz4.m ecambiare il valore di a, in a = 2. Poi eseguo i comandi:

>> hold

Current plot held

>> [t,w]=ode45(@funz4,[0 2],[1 0]);

>> plot(t,w(:,1),’--’)

Infine, con il valore a = 3, eseguo i comandi:

23

>> [t,z]=ode45(@funz4,[0 2],[1 0]);

>> plot(t,z(:,1),’-.’)

Ottengo la figura di sinistra di Fig 4.1. Ingrandisco la zona vicino a T = 2, ottenendo la figura didestra. Dal dettaglio si vede facilmente che la curva che si avvicina di piu al valore -14 in t = 2, e lacurva punteggiata, che e stata ottenuta con a = 1. Oppure, posso calcolare i valori delle tre funzioni int = 2:

>> y(length(y))

ans =

-12.0750

>> z(length(z))

ans =

-19.2569

>> w(length(w))

ans =

-17.8012

Il numero fra quelli appena calcolati, che si avvicina di piu a -14 e -12.0750, che e stato ottenuto cona = 1. E’ possibile rendere piu efficiente la soluzione di questo problema, passando a come parametro

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2−20

−15

−10

−5

0

5

1.8 1.9 2

−19

−18

−17

−16

−15

−14

−13

−12

−11

−10

Figura 4.1: Esercizio 4.1. A sinistra: grafico delle tre soluzioni. A destra: particolare, intorno a t = 2.

alla function fun4.m. Poi, si chiama la function ode45, inserendo nella chiamata il parametro a.

4.2 Esercizio 3, 1/9/03

Considerare il problema di Cauchy:{

y′ = t cos(y2 + t)y(1) = 0

Risolvere il problema differenziale fino a t = 10, usando la ode-function ode45 di Matlab. Determinareil passo di integrazione minimo ed il passo di integrazione massimo utilizzati nel corso dell’integrazione.Soluzione

Innanzitutto e necessario creare una ode-function per il problema assegnato. In questo caso utilizzo ilfile fun3.m

function f=fun3(t,y)

f=t.*cos(y.^2+t);

24

A questo punto, posso integrare il problema differenziale con la function ode45:

[t,y]=ode45(’fun3’,[1,10],0);

Il vettore t contiene gli istanti in corrispondenza dei quali Matlab ha calcolato la soluzione. Se sottraggodue istanti successivi, ottengo quindi i passi di integrazione utilizzati. Quindi le informazioni richiestedall’esercizio possono essere ottenute con il seguente programma:

n=length(t);

for i=1:n-1

dt(i)=t(i+1)-t(i);

end

min(dt)

max(dt)

Il passo minimo e massimo utilizzati sono quindi:

ans =

9.2981e-05

ans =

0.1682

4.3 Esercizio 3, 8/9/03

Considerare il problema di Cauchy:{

y′ = t cos(y2 + t)y(1) = 0

Risolvere il problema differenziale fino a t = 10, usando la ode-function ode45 di Matlab. Determinarei valori massimi e minimi della soluzione nell’intervallo [1, 10].Soluzione

Calcoliamo innanzitutto la soluzione dell’equazione differenziale, disegnandone il grafico, in modo davedere dove si collocano i punti di estremo relativo della funzione y(x).

[t1,y1]=ode45(’fun3’,[1,10],0);

plot(t1,y1)

Si ottiene il grafico che compare sulla sinistra di Fig. 4.2. Dal grafico, si vede immediatamente checi sono due punti di massimo relativo, che si trovano rispettivamente negli intervalli [1, 2] e [5, 6]. Cisono inltre due punti di minimo relativi che si trovano negli intervalli [2, 3] e [8, 9]. Per determinarele coordinate dei punti di estremo, possiamo procedere in due modi. Una possibilita’ e’ di stampare ivettori t e y e localizzare i valori cercati osservando i valori numerici delle componenti di t e di y.

Un’altra possibilita e di far cercare i minimi ed i massimi a Matlab. Per fare questo, utilizziamo lefunctions min e max. Queste functions, applicate a tutto il vettore y pero forniscono soltanto il valoremassimo e minimo assoluti. Dobbiamo quindi trovare a quali componenti del vettore t corrispondanogli estremi degli intervalli nei quali sappiamo gia che si trovano i punti di estremo. Per ottenere questorisultato, disegniamo un grafico di t in funzione della sua componente. Otteniamo il grafico in Fig. 4.2a destra.

plot(t1,’*’)

grid on

Scegliendo zoom opportuni del grafico sugli intervalli che stiamo cercando, si trova facilmente chel’intervallo di tempo [1, 2] corrisponde alle componenti 1 ≤ i ≤ 25 del vettore t, mentre per gli altriintervalli si ha 2 ≤ t(i) ≤ 3 per 25 ≤ i ≤ 42; 5 ≤ t(i) ≤ 6 per 72 ≤ i ≤ 101 e infine 8 ≤ t(i) ≤ 9 per143 ≤ i ≤ 197. Notare che l’istruzione

[ymin,i]=min(y)

25

1 2 3 4 5 6 7 8 9 10−2

−1.5

−1

−0.5

0

0.5

1

1.5

2

0 50 100 150 200 2501

2

3

4

5

6

7

8

9

10

Figura 4.2: Esercizio 4.3: Grafico della soluzione dell’equazione differenziale(sinistra) e grafico di t(n) in funzione di n (a destra).

fornisce il valore minimo del vettore y e la componente in corrispondenza della quale si trova il valoredi minimo. Attenzione: l’istruzione [y1,i]=min(y(i1:i2)) crea il vettore y(i1:i2) e ne calcola ilminimo, con la coordinata relativa, quindi per avere la coordinata del minimo di y1 devo calcolare poii=i+i1-1.

Quindi, per ottenere tutte le informazioni richieste dal problema, posso usare il seguente programma:

[y1min,i1min]=min(y1(25:42))

t1min=t1(i1min+24)

[y2min,i2min]=min(y1(143:197))

t2min=t1(i2min+142)

% trovo ora i valori di massimo

[y1max,i1max]=max(y1(1:25))

t1max=t1(i1max)

[y2max,i2max]=max(y1(72:101))

t2max=t1(i2max+71)

Ottengo cosı i dati richiesti:

y1min =

-1.3416

i1min =

16

t1min =

2.9086

y2min =

-1.5767

i2min =

25

t2min =

8.5057

y1max =

0.1716

i1max =

22

t1max =

26

1.4586

y2max =

1.5024

i2max =

16

t2max =

5.5862

cioe i due valori di minimo si trovano in t ' 2.9086 e in t ' 8.5057 e valgono, rispettivamente,y ' −1.3416 e y ' −1.5767. Infine i valori di massimo relativo si trovano in t ' 1.4586 e in t ' 5.5862e valgono rispettivamente y ' 0.1716 e y ' 1.5024. Naturalmente, questi risultati sono approssimati,perche non conosco la soluzione esatta y(x), ma soltanto una sua approssimazione.

27

Introduzionecalvino.polito.it/~puppo/didattica/dispensa.pdf · Introduzione Questa dispensa ... una...

Documents

Transcript of Introduzionecalvino.polito.it/~puppo/didattica/dispensa.pdf · Introduzione Questa dispensa ... una...