Disp Matlab

85
UNIVERSIT ` A DEGLI STUDI DI FERRARA FACOLT ` A DI INGEGNERIA Corso Matlab per Analisi Matematica II c.d.l. Ingegneria Civile e Ambientale 1 Michele Miranda, Francesco Musacci, Cecilia Pasquini a.a. 2010-2011 1 versione aggiornata al 21 dicembre 2010

description

dispensa

Transcript of Disp Matlab

Page 1: Disp Matlab

UNIVERSITA DEGLI STUDI DI FERRARA

FACOLTA DI INGEGNERIA

Corso Matlab per

Analisi Matematica II

c.d.l. Ingegneria Civile e Ambientale1

Michele Miranda, Francesco Musacci, Cecilia Pasquini

a.a. 2010-2011

1versione aggiornata al 21 dicembre 2010

Page 2: Disp Matlab

2

Page 3: Disp Matlab

Indice

1 Curve 71.1 Richiami di algebra lineare . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2 Alcuni disegni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.2.1 Curve nel piano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.2 Curve in coordinate polari . . . . . . . . . . . . . . . . . . . . . . . . . 91.2.3 Curve nello spazio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.3 Vettori tangenti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.4 Lunghezza d’arco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.4.1 Grafico di funzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131.4.2 Curve parametrizzate . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.5 Curvatura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2 Funzioni di piu variabili 172.1 Funzioni e loro grafici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.1.1 Funzioni discontinue . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2 Derivate parziali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.3 Piani tangenti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.4 Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.4.1 Un esempio nel piano . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.4.2 Un esempio nello spazio . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Superfici 253.1 Superfici cartesiane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2 Superfici di rotazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.3 Superfici Parametrizzate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.3.1 Grafici di funzioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.3.2 Superfici di rotazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.4 Piano tangente ad una superficie parametrizzata . . . . . . . . . . . . . . . . 32

4 Integrali multipli 334.1 Integrali doppi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 334.2 Formule di riduzione per integrali doppi . . . . . . . . . . . . . . . . . . . . . 344.3 Integrali doppi su domini generici . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.3.1 Disegno di insiemi nel piano . . . . . . . . . . . . . . . . . . . . . . . . 364.4 Integrali tripli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.4.1 Integrali su domini generici . . . . . . . . . . . . . . . . . . . . . . . . 384.4.2 Grafico di regioni solide . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3

Page 4: Disp Matlab

4 INDICE

5 Punti stazionari e massimi e minimi 41

5.1 Classificazione di punti stazionari . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.2 Massimo e minimo di una funzione . . . . . . . . . . . . . . . . . . . . . . . . 42

5.2.1 Massimo e minimo di una funzione di due variabili su di un rettangolo 43

5.2.2 Massimi e minimi su cerchi . . . . . . . . . . . . . . . . . . . . . . . . 44

5.3 Massimi e minimi vincolati . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

6 Integrali curvilinei e di superficie 47

6.1 Integrali curvilinei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

6.1.1 Integrali curvilinei per funzioni scalari . . . . . . . . . . . . . . . . . . 47

6.1.2 Rappresentazione di campi vettoriali . . . . . . . . . . . . . . . . . . . 48

6.1.3 Integrali curvilinei per funzioni vettoriali . . . . . . . . . . . . . . . . . 49

6.2 Aree di superfici parametrizzate . . . . . . . . . . . . . . . . . . . . . . . . . . 51

6.3 Integrali di superficie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

6.3.1 Integrali di superficie per funzioni scalari . . . . . . . . . . . . . . . . 52

6.3.2 Rappresentazione di campi lungo una superficie . . . . . . . . . . . . . 53

6.3.3 Integrali di superficie per campi vettoriali: circuitazioni e flussi . . . . 53

7 Successioni e serie di funzioni 57

7.1 Successioni di funzioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

7.2 Serie di funzioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

7.3 Serie di potenze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

7.4 Serie di Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

7.4.1 Polinomi di Taylor simbolici . . . . . . . . . . . . . . . . . . . . . . . . 61

7.4.2 Moltiplicazione di serie e serie di Taylor di prodotti di funzione . . . . 61

7.5 Serie di Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

7.5.1 Fenomeno di Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

8 Equazioni differenziali 67

8.1 Linee di flusso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

8.2 Il metodo numerico ode45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

8.3 Integrazione simbolica di equazioni differenziali . . . . . . . . . . . . . . . . . 69

8.4 Sistemi di equazioni differenziali . . . . . . . . . . . . . . . . . . . . . . . . . 70

8.5 Equazioni di ordine superiore . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

8.6 Soluzioni definite implicitamente . . . . . . . . . . . . . . . . . . . . . . . . . 71

8.7 Vibrazioni meccaniche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

8.7.1 Vibrazioni non smorzate . . . . . . . . . . . . . . . . . . . . . . . . . . 73

8.7.2 Vibrazioni smorzate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

A Derivazione numerica 77

A.1 Derivata prima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

A.2 Derivata seconda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

B Il metodo di integrazione di Simpson 79

C Metodo di Eulero per le equazioni differenziali 81

Page 5: Disp Matlab

INDICE 5

D Integrazione di equazioni differenziali per serie 83D.1 Serie di Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83D.2 Serie di Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Page 6: Disp Matlab

6 INDICE

Page 7: Disp Matlab

Capitolo 1

Curve

1.1 Richiami di algebra lineare

Iniziamo col disegnare il segmento che congiunge due punti dati; ricordiamo che l’equazioneparametrica del segmento che congiunge i punti P,Q ∈ RN e data da

x = tQ+ (1− t)P ;

se i due punti sono nel piano, possiamo scrivere P = (P1, P2), Q = (Q1, Q2). Quindi un modosemplice per disegnare nel piano un segmento che congiunge P = (1, 2) con Q = (−3, 5) puoessere

t=[0 1];

P=[1 2];

Q=[-3 5];

x=(1-t)*P(1)+t*Q(1);

y=(1-t)*P(2)+t*Q(2);

plot(x,y)

A tali comandi si possono aggiungere le opzioni

axis([-4 2 0 3]); grid on;

per migliorarne l’aspetto grafico.Il precedente grafico rappresenta il segmento che congiunge i due punti, o anche il vet-

tore P −Q applicato nel punto Q. Se volgiamo invece dare una rappresentazione grafica diun vettore (quindi come di un segmento che parte dall’origine e congiunge il punto dato) edella somma tra due vettori, possiamo considerare i seguenti comandi, che danno la rapp-resentazione nel piano dei vettore a = (4, 7), b = (−2, 3), del vettore 1/2a e della sommaa+ b.

t=[0 1];

xa=4*t;

ya=7*t;

xb=-2*t;

yb=3*t;

7

Page 8: Disp Matlab

8 CAPITOLO 1. CURVE

plot(xa,ya,xb,yb,.5*xa,.5*ya,xa+xb,ya+yb);

axis([-4 5 -2 12]);

grid on;

legend(’a’,’b’,’.5*a’,’a+b’)

Se vogliamo disegnare vettori in dimensione 3 bisognera utilizzare il comando

plot3

al posto di

plot.

Ad esempio, possiamo considerare i vettori a = (2, 3,−1), b = (1,−3, 5), c = (−2, 4,−1),a− b e a+ 2b− 3c nello spazio;

t=[0 1];

xa=2*t; ya=3*t; za=-t;

xb=t; yb=-3*t; zb=5*t;

xc=-2*t; yc=4*t; zc=-t;

xs=xa-xb; ys=ya-yb; zs=za-zb;

xx=xa+2*xb-3*xc; yy=ya+2*yb-3*yc; zz=za+2*zb-3*zc;

plot3(xa,ya,za,xb,yb,zb,xc,yc,zc,xs,ys,zs,xx,yy,zz);

grid on;

legend(’a’,’b’,’c’,’a-b’,’a+2b-3c’)

Possiamo anche utilizzare il prodotto scalare per calcolare l’angolo tra due vettori dati;ricordiamo che vale la formula

x · yn∑i=1

xiyi = ‖x‖‖y‖ cosϑ

per ogni x, y ∈ RN , da cui

ϑ = arccosx · y‖x‖‖y‖

.

Osserviamo che quest’ultima formula restituisce sempre un angolo compreso tra 0 e π; quindidegli angoli individuati dai due vettori x e y, con la formula precedente si sta determinandoil minore dei due, mentre il secondo sara dato da 2π − ϑ. Un semplice esempio e dato daiseguenti comandi, in cui si determina l’angolo tra i due vettori (2, 3,−1) e (1,−3, 2).

a=[2 3 -1];

b=[1 -3 5];

La=sqrt(sum(a.*a));

Lb=sqrt(sum(b.*b));

theta=acos(sum(a.*b)/(La*Lb))

1.2 Alcuni disegni

1.2.1 Curve nel piano

Possiamo cominciare col disegnare un segmento, determinato da un punto di partenza e dauna direzione. Come punto di partenza prendiamo x = (1,−3) e come direzione fissiamo ilvettore v = (2, 4); disegniamo quindi il segmento x+ tv con t ∈ [−2, 2];

Page 9: Disp Matlab

1.2. ALCUNI DISEGNI 9

t=[-2 2];

x1=1+2*t;

x2=-3+4*t;

plot(x1,x2);

grid on;

Si provi a modificare l’insieme di definizione di t, tracciando alcuni esempi.

Esercizio 1.1 Disegnare le seguenti curve piane;

1. (t− 2, t+ sin t), t ∈ R;

2. (tet, log(1 + t)), t ∈ (−1,+∞);

3. (t|t|, t− cos t), t ∈ R;

4. (t2, t3), t ∈ [−1, 1];

5. (cos3 t, sin3 t), t ∈ [0, 2π];

6. (Rt− r sin t, R− r cos t), t ∈ [0, 4π] con diversi valori di R e r (R > r, R = r, R < r);

1.2.2 Curve in coordinate polari

Le curve in coordinate polari, cioe le curve in cui si pone il raggio in funzione dell’angolo% = %(ϑ), possono essere disegnate semplicemente considerando la parametrizzazione

r(ϑ) = %(ϑ)(cosϑ, sinϑ).

In Matlab esiste pero anche un comando speciale, polar; disegniamo come esempio la rosaa quattro petali data %(ϑ) = cos 2ϑ, ϑ ∈ [0, 2π].

theta=linspace(0,2*pi);

r=cos(2*theta);

polar(theta,r);

Possiamo anche disegnare una curva in coordinate polari definita da un parametro c perdue differenti valori di c, ad esempio c = 1 e c = 2; disegniamo quindi la curva %(ϑ) =1 + c sinϑ, ϑ ∈ [0, 2π].

theta=linspace(0,2*pi);

polar(theta,1+sin(theta),’c’);

hold on;

polar(theta,1+2*sin(theta),’r:’);

hold off;

legend(’1+sen(t)’, ’1+2sen(t)’);

Purtroppo, il comando polar non formatta bene il disegno, in quanto scala il graficobasandosi solo sul primo disegno. Possiamo ovviare questo problema nel modo seguente.

theta=linspace(0,2*pi);

r1=1+sin(theta); r2=1+2*sin(theta);

polar([theta’,theta’], [r1’,r2’]);

legend(’1+sen(t)’, ’1+2sen(t)’);

Page 10: Disp Matlab

10 CAPITOLO 1. CURVE

Esercizio 1.2 Disegnare le seguenti curve piane;

1. % = 2a(1 + cosϑ), ϑ ∈ [−π, π];

2. % = aeϑ, ϑ ∈ [0, π/6];

3. % = ϑ2, ϑ ∈ [0, 3/2];

4. % = (1− cosϑ), ϑ ∈ [0, 2π].

1.2.3 Curve nello spazio

Segmento Possiamo considerare la generalizzazione naturale dell’esempio considerato nel-la sezione 1.2.1 disegnando un segmento, determinato dal punto di partenza x = (1,−3, 2) edirezione data dal vettore v = (2, 4,−1); disegniamo quindi il segmento x+tv con t ∈ [−2, 2];

t=[-2 2];

x1=1+2*t;

x2=-3+4*t;

x3=2-t;

plot3(x1,x2,x3);

grid on;

Si provi a modificare l’insieme di definizione di t, tracciando alcuni esempi.

Elica cilindrica L’elica o spirale cilindrica e descritta dell’equazione parametrica

r(t) = (a cos t, a sin t, kt), t ∈ R

e a, k che descrivono il raggio e il passo dell’elica rispettivamente. Disegniamo l’elica cona = 2 e k = 3;

a=2; k=3;

t=linspace(0,4*pi);

plot3(a*cos(t),a*sin(t),k*t);

grid on;

Possiamo anche disegnare il sostegno dell’elica di raggio 2 e passo 1/2 e contemporanea-mente disegniamo la sua proiezione sul piano xy.

a=2; k=1/2;

t=linspace(0,4*pi);

x=a*cos(t); y=a*sin(t); z=k*t;

plot3(x,y,z,’r’,x,y,0*t,’y’);

grid on;

Il precedente grafico puo essere anche arricchito nel seguente modo:

hold on;

s=linspace(0,4*pi,30);

x=a*cos(s); y=a*sin(s); z=k*s;

stem3(x,y,z,’b’);

hold off;

Page 11: Disp Matlab

1.3. VETTORI TANGENTI 11

Esercizio 1.3 Disegnare la curva nello spazio determinata dal luogo delle soluzioni delsistema {

y = x2

3z = 2xy.

1.3 Vettori tangenti

Possiamo disegnare la curva piana r(t) = (4 cos t, 2 sin t), t ∈ [0, 2π], e il vettore tangentealla curva in r(π/6), dato da r′(π/6). Disegneremo anche i vettori secanti la curva

r(π/6 + h)− r(π/6)

h

con h = 0.1, 0.05, 0.01.

t=linspace(0,2*pi);

x=4*cos(t); y=2*sin(t);

xp=-4*sin(pi/6); yp=2*cos(pi/6);

s=[0 1];

v1=4*cos(pi/6)+s.*xp;

v2=2*sin(pi/6)+s.*yp;

plot(x,y,v1,v2,’r’);

hold on;

for h=[.1 .05 .01];

u1=4*cos(pi/6)+s.*(4*cos(pi/6+h)-4*cos(pi/6))./h;

u2=2*sin(pi/6)+s.*(2*sin(pi/6+h)-2*sin(pi/6))./h;

plot(u1,u2);

end

hold off;

legend(’curva’,’vettore tangente’,’vettori secanti’,2);

Riprendiamo ora l’esempio dell’elica cilindrica e disegniamo il vattore tangente alla curvain alcuni suoi punti; disegniamo al contempo anche il vettore normale principale, definitodalla formula

T ′(t)

‖T ′(t)‖,

dove T (t) = r(t)/‖r(t)‖ definisce il versore tangente. Possiamo quindi considerare i seguenticomandi:

clear

a=2; k=3;

t=linspace(0,4*pi);

x=a*cos(t); y=a*sin(t); z=k*t;

plot3(x,y,z);

grid on; axis equal; hold on;

xlabel(’x’);ylabel(’y’);zlabel(’z’);

% grafico dell’elica

tt=0:pi/4:4*pi; l=length(tt);

% definisco un vettore di valori compresi tra 0 e 4*pi

Page 12: Disp Matlab

12 CAPITOLO 1. CURVE

xx=a*cos(tt); yy=a*sin(tt); zz=k*tt;

% calcolo r(t), r’(t) e r’’(t) per questi valori

dxx=-a*sin(tt); dyy=a*cos(tt); dzz=k*ones(l);

ddxx=-a*cos(tt); ddyy=-a*sin(tt); ddzz=zeros(l);

s=[0 1];

for i=1:l

tg=norm([dxx(i) dyy(i) dzz(i)]);

% per i valori in tt calcolo il vettore tangente normalizzato

tg1=xx(i)+s.*dxx(i)./tg;

% cioe’ il versore tangente

tg2=yy(i)+s.*dyy(i)./tg;

tg3=zz(i)+s.*dzz(i)./tg;

n=norm([ddxx(i) ddyy(i) ddzz(i)]);

% calcolo inoltre il vettore normale normalizzato

n1=xx(i)+s.*ddxx(i)./n;

% cioe’ il versore normale

n2=yy(i)+s.*ddyy(i)./n;

n3=zz(i)+s.*ddzz(i)./n;

plot3(tg1,tg2,tg3,’-r’,n1,n2,n3,’-g’);

pause

end

hold off;

In alernativa si puo considerare un solo giro d’elica, come segue;

clear

a=2; k=3;

t=linspace(0,2*pi);

x=a*t.*cos(t); y=a*t.*sin(t); z=k*t;

plot3(x,y,z);

grid on; axis equal; hold on;

xlabel(’x’);ylabel(’y’);zlabel(’z’);

% grafico dell’elica

tt=0:pi/8:4*pi; l=length(tt);

% definisco un vettore di valori compresi tra 0 e 4*pi

xx=a*tt.*cos(tt); yy=a*tt.*sin(tt); zz=k*tt;

% calcolo r(t), r’(t) e r’’(t) per questi valori

dxx=a*(cos(tt)-tt.*sin(tt));

dyy=a*(sin(tt)+tt.*cos(tt));

dzz=k*ones(l);

ddxx=-a*(2*sin(tt)+tt.*cos(tt));

ddyy=a*(2*cos(tt)-tt.*sin(tt));

ddzz=zeros(l);

s=[0 1];

for i=1:l

tg=norm([dxx(i) dyy(i) dzz(i)]);

Page 13: Disp Matlab

1.4. LUNGHEZZA D’ARCO 13

% per i valori in tt calcolo il vettore tangente normalizzato

tg1=xx(i)+s.*dxx(i)./tg;

% cioe’ il versore tangente

tg2=yy(i)+s.*dyy(i)./tg;

tg3=zz(i)+s.*dzz(i)./tg;

n=norm([ddxx(i) ddyy(i) ddzz(i)]);

% calcolo inoltre il vettore normale normalizzato

n1=xx(i)+s.*ddxx(i)./n;

% cioe’ il versore normale

n2=yy(i)+s.*ddyy(i)./n;

n3=zz(i)+s.*ddzz(i)./n;

plot3(tg1,tg2,tg3,’-r’,n1,n2,n3,’-g’);

pause

end

hold off;

1.4 Lunghezza d’arco

Come primo esempio di calcolo di una lunghezza, vediamo come la lunghezza di una polig-onale inscritta in una semicirconferenza puo essere usata per approssimare il valore di π.Consideriamo quindi i seguenti comandi;

t=linspace(0,pi);

for i=[1,2,3,4];

s=linspace(0,pi,3^i+1);

x=cos(s); y=sin(s);

dx=diff(x); dy=diff(y);

ds=sqrt(dx.^2+dy.^2);

L=sum(ds);

plot(cos(t),sin(t));

hold on;

plot(x,y);

T=[’Lunghezza della Poligonale= ’,num2str(L)];

text(-40,4,T);

hold off;

pause;

end;

1.4.1 Grafico di funzione

In questa sezione calcoliamo la lunghezza di una curva cartesiana, cioe la lunghezza delgrafico di una funzione y = f(x), x ∈ [a, b]. Ricordiamo che tale lunghezza e data dallaformula

L = l(r, [a, b]) =

∫ b

a

√1 + f ′(x)2dx,

Page 14: Disp Matlab

14 CAPITOLO 1. CURVE

dove r : [a, b] → R2 e la parametrizzazione r(x) = (x, f(x)). Solitamente, l’integraleprecedente non si riesce a calcolare analiticamente; se ad esempio consideriamo la funzionef(x) = x3 + 3x2 − 5x+ 2, x ∈ [−1, 2], ci ritroviamo col seguente integrale

L =

∫ 2

−1

√36x4 + 36x3 + 6x2 − 60x+ 26dx.

Determineremo il valore di L considerando le approssimazioni del grafico di f mediantepoligonali; tali poligonali avranno vertici nei punti di coordinate (x, f(x)); procederemoquindi con la suddivisione dell’intervallo [−1, 2] in segmenti della medesima ampiezza. Talecalcolo viene effettuato ad esempio dal seguente semplice script.

x=linspace(-1,2,400);

y=x.^3+3*x.^2-5*x+2;

dx=diff(x);

dy=diff(y);

ds=sqrt(dx.^2+dy.^2);

L=sum(ds)

In tale modo si ottiene il valore L = 20.8314; dire cosa si ottiene se l’ultimo comandoviene sostituito con

L=cumsum(ds)

Il vantaggio dell’approssimazione poligonale e quello di ottenere in modo abbastanzasemplice un valore approssimato del valore della lunghezza della curva; non si ha pero nes-suna informazione sull’errore commesso nell’approssimazione. Come vedremo nel prossimoparagrafo, questo errore puo essere anche consistente.

1.4.2 Curve parametrizzate

Iniziamo con il calcolo approssimato del numero π;Vediamo ora col seguente esempio quanto l’approssimazione poligonale puo non essere

accurata. Consideriamo la curva

r(t) = (4t− 2t3, t2), t ∈ [−1, 3].

Utilizziamo l’approssimazione poligonale dove viene suddiviso l’intervallo [−1, 3] in 400 partiuguali e la poligonale ha per vertici i punti r(ti); considereremo inoltre anche un metodonumerico per il calcolo del seguente integrale:

L = l(r, [−1, 3]) =

∫ 3

−1

√36t4 − 32t2 + 16dt,

che e la formula per il calcolo della lunghezza di una curva regolare.

t=linspace(-1,3,400);

x=4*t-2*t.^3;

y=t.^2;

plot(x,y);

dx=diff(x);

Page 15: Disp Matlab

1.4. LUNGHEZZA D’ARCO 15

dy=diff(y);

ds=sqrt(dx.^2+dy.^2);

L1=sum(ds);

L2=quad(fcnchk(’sqrt(36*t.^4-32*t.^2+16)’),-1,3,.0001);

T1=[’Lunghezza poligonale approssimante= ’,num2str(L1)];

T2=[’Lunghezza approssimata al decimillesimo= ’,num2str(L2)];

text(-40,4,T1);

text(-40,3.5,T2);

Si noti la differenza tra i due valori ottenuti; si puo provare a sostituire il valore 400con altri valori e vedere per quale valore di tale parametro la lunghezza della poligonale eparagonabile col valore L2.

Una variante del precedente esercizio puo essere la determinazione della lunghezza del-l’ellisse di equazione

x2

4+y2

9= 1,

che si puo parametrizzare mediante r(t) = (2 cos t, 3 sin t), t ∈ [0, 2π].

t=linspace(0,2*pi,1001);

x=2*cos(t); y=3*sin(t);

plot(x,y);

dx=diff(x); dy=diff(y); ds=sqrt(dx.^2+dy.^2);

L=sum(ds);

qT=[’Lunghezza poligonale approssimante= ’,num2str(L)];

text(-1.5,.2,T);

Ovviamente, come prima, possiamo approssimare il valore dell’integrale∫ 2π

0

√6− 5 sin2 tdt

utilizzando metodi numerici; useremo qui il comando quad8. Possiamo quindi aggiungereai precedenti comandi i seguenti;

L2=quad8(fcnchk(’sqrt(9-5*sin(t).^2)’),0,2*pi);

text(-1.5,.2,’Valore approssimato integrale= ’,num2str(L2));

Esercizio 1.4 Provare a ripetere i conti precedenti per le seguenti curve; si provi anche avedere se la lunghezza della curva puo essere determinata utilizzando il symbolic toolbox.

1. y = senx per x ∈ [0, π/2];

2. x = cot t, y = sen 2t per t ∈ [π/4, π/2];

3. y = x3 per x ∈ [0, 2];

4. x = t cos t, y = t sen t per t ∈ [0, 2π].

Page 16: Disp Matlab

16 CAPITOLO 1. CURVE

1.5 Curvatura

Consideriamo qui la curva parametrizzata da

r(t) = (2 sin(t), cos t, cos t);

creeremo una figura divisa in due parti, dove nella prima parte tracceremo il sostegno di r,mentre nella seconda parte disegneremo il grafico della funzione curvatura

k(t) =‖r′(t)× r′′(t)‖‖r′(t)‖3

=1

(1 + cos2 t)3/2.

t=linspace(0,2*pi,200);

subplot(2,1,1);

plot3(2*sin(t),cos(t),cos(t));

grid on;

text(0,1,1-.2,’\uparrow t=0’)

text(2,0,0-.2,’\uparrow t=\pi/2’)

text(-2,0,0-.2,’\uparrow t=3\pi/2’)

subplot(2,1,2);

plot(t,1./(1+cos(t).^2).^(3/2));

set(gca,’XTick’,[0:pi/2:2*pi])

Tramite il calcolo della curvatura, possiamo anche ricavare il cerchio osculatore allacurva. Prenderemo in considerazione il caso di una curva piana descritta come grafico dellafunzione f(x) = x2, x ∈ [−2, 2]. Disegneremo il sostegno della curva e il cerchio osculatorela curva nel punto (1/2, 1/4). La curvatura per una curva cartesiana e data da

k(x) =|f ′′(x)|

(1 + |f ′(x)|2)3/2=

2

(1 + 4x2)3/2.

Per x = 1/2 otteniamo k(1/2) = 1/√

2; dato che il versore tangente in tale punto e dato da1√2(1, 1) mentre il versore normale e dato da 1√

2(−1, 1), il cerchio osculatore avra raggio

√2

e sara centrato in (− 12 ,

54 ).

x=linspace(-2,2,400);

t=linspace(0,2*pi,600);

plot(x,x.^2,-.5+sqrt(2)*cos(t),1.25+sqrt(2)*sin(t));

axis equal;

Esercizio 1.5 1. Far variare il cerchio osculatore nell’ultimo esempio considerato.

2. Si approssimi la lunghezza di r(t) = (t, t2, t3), t ∈ [0, 2] sia con l’approssimazione polig-onale dividendo l’intervallo in 1000 sottointervalli uguali, quindi usando il comandoquad.

3. Determinare k(t) per r(t) = (sin t, cos t, sin2 t), t ∈ [0, 2π] e, mediante il comandosubplot, disegnare la curva e il grafico della sua curvatura.

4. Disegnare la curva cartesiana y = log x, x ∈ [2, 4] e il cerchio osculatore alla curva nelpunto relativo ad x = 2.

5. Disegnare il sostegno della curva r(t) = (t, t2, t3), t ∈ [0, 2] e il cerchio osculatorerelativo al valore t = 1.

Page 17: Disp Matlab

Capitolo 2

Funzioni di piu variabili

2.1 Funzioni e loro grafici

Iniziamo col disegnare il grafico di alcune funzioni di due variabili z = f(x, y); possiamo

ad esempio considerare la funzione z =√

4− 2x2 − y2, definita nella regione 2x2 + y2 ≤ 4(ellisse di semiassi

√2 e 2). Per definire quindi il dominio nelle variabili (x, y) tramite il

comando meshgrid abbiamo dei problemi; conviene quindi passare alle coordinate polarix = r

√2 cosϑ, y = 2r sinϑ, con (r, ϑ ∈ [0, 1] × [0, 2π], regione rettangolare). Se si prova

pero a far variare r fino al valore 1 si riscontrano problemi, in quanto la funzione si annullae Matlab genera problemi; per questo faremo variare r fino a 0.95.

[r,theta]=meshgrid(0:.05:.95,0:2*pi/30:2*pi);

x=sqrt(2)*r.*cos(theta);

y=2*r.*sin(theta);

z=sqrt(4-2.*x.^2-y.^2);

surf(x,y,z); axis equal;

Le righe che si vedono in figura altro non sono che le restrizioni della funzione sulle curve ar costante e ϑ costante. In particolare le linee che dividono le regioni con colori differenti (i“paralleli“) sono ottenute tramite ϑ 7→ f(r, ϑ), con f(r, ϑ) = f(

√2r cosϑ, 2r sinϑ), mentre

le linee che escono a raggiera dal punto di massima elevazione (i ”meridiani“) sono ottenutitramite r 7→ f(r, ϑ).

Se non si volessero utilizzare le coordinate polari, si puo restare nelle coordinate (x, y),moltiplicando la funzione per 0 quando 4− 2x2 − y2 < 0. Si puo procedere come segue.

[x,y]=meshgrid(-sqrt(2):.1:sqrt(2),-2:.1:2);

z=sqrt((4-2.*x.^2-y.^2 >= 0).*(4-2.*x.^2-y.^2));

surfl(x,y,z);

Notiamo qui che l’operazione (4− 2. ∗ x.2 − y.2 >= 0) e un operatore logico, che restituisceil valore 1 se la condizione e verificata, altrimenti restituisce 0. Otteniamo quindi unamatrice di valori che sono o 1 o 0 a seconda che il punto (x, y) appartenga o meno all’ellisse.L’operazione consiste quindi semplicemente nel disegnare il grafico della funzione

f(x, y) =√χE(x, y)(4− 2x2 − y2), E = {(x, y) ∈ R2 : 2x2 + y2 ≤ 4}

17

Page 18: Disp Matlab

18 CAPITOLO 2. FUNZIONI DI PIU VARIABILI

e χE e la funzione caratteristica dell’insieme E, cioe

χE(x, y) =

1 se (x, y) ∈ E

0 se (x, y) 6∈ E.

In questo modo la funzione χE(x, y)(4 − 2x2 − y2) e definita per ogni (x, y) ∈ R2 e nonabbiamo problemi a disegnarne il grafico. Concludiamo questa sezione con un commento; lelinee che sono visualizzate nel grafico ottenuto altro non sono che le restrizioni della funzionex 7→ f(x, y) e y 7→ f(x, y).

2.1.1 Funzioni discontinue

Iniziamo con una funzione che abbiamo gia visto durante il corso di teoria, e cioe la funzione

f(x, y) =xy

x2 + y2, (x, y) 6= (0, 0).

Utilizziamo le coordinate polari per disegnare tale grafico;

[r,theta]=meshgrid(0:.02:1,0:2*pi/50:2*pi);

x=r.*cos(theta);

y=r.*sin(theta);

z=x.*y./(x.^2+y.^2);

surf(x,y,z); axis equal;

Dal grafico si intuisce la discontinuita della funzione. Possiamo visualizzare meglio questofenomeno tramite le restrizioni y = mx, che corrisponde in coordinate polari, alle sezionir 7→ f(r, ϑ), dove f(r, ϑ) = f(r cosϑ, r sinϑ).

Le sezioni ad angolo costante si ottiene mediante;

r=linspace(0,1);

plot3(r,r.*0,r.*0);

grid on;

axis([-1 1 -1 1 -.5 .5]);

hold on;

for i=[0:2*pi/50:2*pi]

plot3(r.*cos(i),r.*sin(i),r.^2.*cos(i)*sin(i)./r.^2);

pause(.1)

end;

hold off;

mentre le sezioni a raggio fissato si ottengono tramite;

theta=linspace(0,2*pi);

plot3(cos(theta),sin(theta),cos(theta).*sin(theta),’r’); grid on;

hold on

for r=[0:.1:1]

plot3((1-r)*cos(theta),(1-r)*sin(theta),cos(theta).*sin(theta),’r’);

pause(.5);

end;

hold off;

Page 19: Disp Matlab

2.1. FUNZIONI E LORO GRAFICI 19

Si provi a ripetere i comandi precedenti, cioe disegnando prima le sezioni ad angolo costantee poi quelle a raggio costante, togliendo il comando hold off; alla fine delle sezioni ad angolocostante.

Possiamo anche ricostruire il grafico della funzione utilizzando il metodo degli insiemi dilivello. Ad esempio possiamo scrivere;

[x,y]=meshgrid(-2:.1:2);

contour(x,y,x.*y./(x.^2+y.^2+eps),10);

Il numero 10 nel comando precedente specifica quante linee di livello si desiderano disegnare;tale parametro puo essere omesso. La variabile eps serve per evitare la divisione per 0; euna variabile Matlab. Possiamo mettere anche le etichette sui livelli; ad esempio possiamoconsiderare:

[x,y]=meshgrid(-2:.1:2);

[C h]=contour(x,y,x.*y./(x.^2+y.^2+eps));

clabel(C,h);

Si puo infine selezionare il livello o mettere l’etichetta solo al livello che si desidera. Vedremoquesto nel prossimo paragrafo.

Possiamo visualizzare il comando contour come sezione del grafico della funzione con ipiano orizzontali z = c; possiamo ad esempio considerare i seguenti comandi.

[x,y]=meshgrid(-2:.1:2);

z=x.*y./(x.^2+y.^2+eps);

[u,v]=meshgrid([-2 2]);

subplot(2,4,1);

surf(x,y,z); hold on; surf(u,v,-0.3*ones(2,2)); title(’z=-0.3’);

subplot(2,4,2);

[C h]=contour(x,y,z,[-0.3]); clabel(C,h);

subplot(2,4,3);

surf(x,y,z); hold on; surf(u,v,-0.1*ones(2,2)); title(’z=-0.1’);

subplot(2,4,4);

[C h]=contour(x,y,z,[-0.1]); clabel(C,h);

subplot(2,4,5);

surf(x,y,z); hold on; surf(u,v,0.1*ones(2,2)); title(’z=0.1’);

subplot(2,4,6);

[C h]=contour(x,y,z,[0.1]); clabel(C,h);

subplot(2,4,7);

surf(x,y,z); hold on; surf(u,v,0.3*ones(2,2)); title(’z=0.3’);

subplot(2,4,8);

[C h]=contour(x,y,z,[0.3]); clabel(C,h);

Per mettere assieme il grafico della funzione e il disegno dei livelli, possiamo considerarei seguenti comandi:

[x,y]=meshgrid(-2:.05:2);

surfc(x,y,x.*y./(x.^2+y.^2+eps));

Page 20: Disp Matlab

20 CAPITOLO 2. FUNZIONI DI PIU VARIABILI

Si provi per esercizio a ripetere la procedura precedente considerando le sezioni a x = ce y = c. Si creino quindi dei subplot in cui da una parte viene disegnato il grafico dellafunzione f e il piano verticale x = c, e da una parte il grafico della funzione f(c, y); si ripetaancora l’esercizio con y = c e f(x, c).

Un secondo esempio Proseguiamo ora disegnando il grafico della funzione, discontinuain (0, 0), definita da

f(x, y) =xy2

x2 + y4.

Possiamo provare direttamente a disegnare il grafico di f :

[x,y]=meshgrid(-2:.1:2);

surf(x,y,(x.*y.^2)./(x.^2+y.^4));

Il disegno viene comunque prodotto, anche se con un buco attorno all’origine (0, 0). Talebuco puo essere eliminato mediante la variabile eps, che e il piu piccolo numero macchinapositivo Matlab.

[x,y]=meshgrid(-2:.1:2);

surf(x,y,(x.*y.^2)./(x.^2+y.^4+eps));

Con i seguenti comandi disegneremo alcuni livelli della funzione data e metteremo l’etichettasolo ad una selezione speciale di livelli:

[x,y]=meshgrid(-2:.05:2);

L=[-.45:.05:.45];

[C,h]=contour(x,y,(x.*y.^2)./(x.^2+y.^4+eps),L);

v=[-.4:.2:.4]; clabel(C,h,v);

Si ripetano infine come esercizio le cose viste nella sezione precedente.

Esercizio 2.1 1. Disegnare i grafici delle funzioni

xy

x2 + y2,

x2y

x2 − y2,

y2

x2 + y2,

x2y

x4 + y2.

2. Tracciare le curve di livello delle precedenti funzioni.

3. Scrivere le superfici di livello delle funzioni di tre variabili

x2 + y2

z2, |x|+ |y|+ |z|.

4. Disegnare il grafico di f(x, y) = xy3 − x3y sul dominio (x, y) ∈ [−3, 3]2; creare poiseparatamente una mappa dei livelli di tale funzione su tale dominio. Usare infine ilcomando surfc.

5. Usare il comando contour per disegnare la curva definita implicitamente dall’equazione:

2x2 + xy + y4 =1

2.

Page 21: Disp Matlab

2.2. DERIVATE PARZIALI 21

2.2 Derivate parziali

Possiamo utilizzare il symbolic toolbox di Matlab per calcolare le derivate parziali diuna funzione. Useremo ad esempio il comando ezplot (su veda anche l’appendice A permaggiori dettagli sulla derivazione numerica); i comandi seguenti calcolano le derivate dif(x, y) = sin(x2 + y2) e tracciamo il grafico di f(1, y), ∂xf(x, 1), ∂2

xxf(1, y) e ∂2xyf(x, 1).

syms x y

f=sin(x^2+y^2); fx=diff(f,’x’);

fxx=diff(f,’x’,2); fxy=diff(fx,’y’);

subplot(2,2,1); ezplot(subs(f,’x’,1));

subplot(2,2,2); ezplot(subs(fx,’y’,1));

subplot(2,2,3); ezplot(subs(fxx,’x’,1));

subplot(2,2,4); ezplot(subs(fxy,’y’,1));

Supponiamo ora di voler disegnare il grafico della funzione ∂2xyf utilizzando il comando

surf. Dobbiamo convertire la funzione ∂2xyf definita precedentemente come una espressione

simbolica in una funzione buona per poter utilizzare surf. Possiamo procedere, a partire daicomandi precedenti, come segue;

Fxy=fcnchk(char(fxy),’x’,’y’,’vectorized’);

[x,y]=meshgrid(-pi:2*pi/40:pi);

z=Fxy(x,y);

figure(2); clf; surf(x,y,z);

Esercizio 2.2 Si ripeta quanto fatto ora per la funzione f(x, y) = sin2(2x+3y); si traccino igrafici di f(x, 2), ∂xf(x, 2), ∂2

xyf(x, 2) e ∂3xyyf(x, 2). Si disegni infine il grafico di ∂3

xyyf(x, y),(x, y) ∈ [−π, π]2.

2.3 Piani tangenti

Vediamo ora come disegnare il grafico di una funzione e del piano tangente al grafico in unpunto fissato; consideriamo ad esempio la funzione f(x, y) = x2 + y2 e il punto (1, 2, 5); laformula per il piano tangente e data da

z = f(x0, y0) +∇f(x0, y0) · (x− x0, y − y0).

Possiamo quindi considerare i seguenti comandi:

[x,y]=meshgrid(-3:.1:3);

surf(x,y,x.^2+y.^2);

hold on;

[u,v]=meshgrid(-2:2:2);

L=5+2*u+4*v;

surf(u+1,v+2,L); hold off;

view(-121,20);

Si noti che la griglia per le variabili (u, v) e meno dettagliata; infatti, essendo L un piano,servono meno punti per avere un buon grafico.

Possiamo anche considerare la funzione x2 − y2 nel punto (1, 2,−3). Con i seguenti co-mandi disegneremo il grafico della funzione per (x, y) ∈ [−3, 3]2 assieme al suo piano tangentein (1, 2,−3) ed attiveremo il comando rotate3d; si provino gli effetti di tale comando.

Page 22: Disp Matlab

22 CAPITOLO 2. FUNZIONI DI PIU VARIABILI

[x,y]=meshgrid(-3:.1:3);

surf(x,y,x.^2-y.^2); hold on;

[u,v]=meshgrid(-2:2:2);

L=-3+2*u-4*v;

surf(1+u,2+v,L); hold off;

rotate3d on;

2.4 Gradiente

Consideriamo in questa sezione il campo gradiente, che come si e visto e un campo ortogonaleagli insiemi di livello della funzione. Inizieremo con una funzione di due variabili per passarepoi ad una funzione di tre variabili.

2.4.1 Un esempio nel piano

Data la funzione f(x, y) = x2 + 5y2 con (x, y) ∈ [−2, 2]2, tracceremo le sue curve di livelloe nello stesso disegno tracceremo il campo gradiente utilizzando la funzione quiver;

[x,y]=meshgrid(-2:.05:2);

L=5*[0:.2:2].^2; % questo vettore definisce la scelta dei livelli

[c,h]=contour(x,y,x.^2+5*y.^2,L); hold on;

[x,y]=meshgrid(-2:.2:2);

dx=2*x; dy=10*y;

quiver(x,y,dx,dy); hold off;

Un secondo esempio e dato dalla funzione f(x, y) = 2x2 + xy + y2 nel dominio [−2, 2]×[−3, 3]; tracciamo quindi nello stesso grafico il campo gradiente di f e un insieme di livello;

[x,y]=meshgrid(-2:.2:2,-3:.2:3);

fx=4*x+y; fy=x+2*y;

quiver(x,y,fx,fy); hold on;

contour(x,y,2*x.*y+y.^2); hold off;

2.4.2 Un esempio nello spazio

Consideriamo qui la funzione f(x, y, z) = x2 + y2 − z2 e vogliamo disegnare il livello E1 ={f = 1}; nello stesso disegno considereremo il campo gradiente di f dato da ∇f(x, y, z) =(2x, 2y,−2z). Il livello E1 e descritto dall’equazione

x2 + y2 = z2 + 1,

che significa che per ogni z fissato i punti (x, y) appartengono ad una circonferenza di raggio√z2 + 1,

[theta,z]=meshgrid(0:pi/20:2*pi,-2:.1:2);

x=sqrt(z.^2+1).*cos(theta);

y=sqrt(z.^2+1).*sin(theta);

surf(x,y,z); hold on;

quiver3(x,y,z,2*x,2*y,-2*z); hold off;

rotate3d on;

Page 23: Disp Matlab

2.4. GRADIENTE 23

Un problema che si nota e che la superficie copre i vettori normali; potremmo anche disegnaresolo i vettori gradiente, la superficie sottesa si dedurra dai punti di applicazione di tali vettori;

[theta,z]=meshgrid(0:pi/20:2*pi,-2:.1:2);

x=sqrt(z.^2+1).*cos(theta);

y=sqrt(z.^2+1).*sin(theta);

quiver3(x,y,z,2*x,2*y,-2*z);

rotate3d on;

Esercizio 2.3 1. Si disegni, nello stesso grafico, i livelli di f(x, y) = x2 − y2 a il campogradiente di f .

2. Si disegni la curva x2 − y2 = 1 parametrizzandola (si usino seno e coseno iperbolico)e si disegni nello stesso grafico il campo gradiente di f .

3. Disegnare il livello f(x, y, z) = x2 + 3y2 + 4z2 = 9 utilizzando le relazioni;{x =√

9− 4z2 cosϑ

y =√

9− 4z2 sinϑ

e nello stesso grafico si disegni il campo gradiente di f .

Page 24: Disp Matlab

24 CAPITOLO 2. FUNZIONI DI PIU VARIABILI

Page 25: Disp Matlab

Capitolo 3

Superfici

Iniziamo con un semplice esempio e disegniamo un piano nello spazio; considereremo il pianodi equazione x+2y+z = 3, che puo essere considerato come il grafico z = 3−x−2y = f(x, y)di una funzione di due variabili. Disegneremo il grafico di questa funzione sul dominiorettangolare x ∈ [−4, 4], y ∈ [−2, 2]; per fare questo utilizzeremo il comando meshgrid.

[x,y]=meshgrid(-4:8:4,-2:4:2);

z=3-x-2*y;

surf(x,y,z); view(120,20);

Il comando view serve per modificare il punto da cui viene visualizzata la superficie (i valoridi default sono −37.5o di azimuth, 30o di elevazione). In alternativa al comando surf per lavisualizzazione si puo utilizzare il mesh, come mostra il seguente esempio;

[x,y]=meshgrid(-4:4:4,-2:2:2);

z=3-x-2*y;

mesh(x,y,z); view(120,20);

Per visualizzare un piano possiamo anche utilizzare l’equazione parametrica del pi-ano; disegniamo quindi il piano che passa per il punto (1, 3,−2) e determinato dai vettori(1/5, 3/10,−1/2) e (2/5,−1/3,

√2), con i parametri s, t ∈ [−2, 2]. Possiamo quindi scrivere:

[s,t]=meshgrid(-2:4:2);

x=1+.2*s+.4*t;

y=3+.3*s-1/3*t;

z=-2-.5*s+sqrt(2)*t;

figure(1); mesh(x,y,z); view(120,60);

figure(2); surf(x,y,z); view(120,60);

Per aiutare la visualizzazione, si provi ad utilizzare il comando rotate3d invece di utilizzareil comando view.

Visualizziamo inoltre il grafico di due piani; con i seguenti comandi determiniamo anchel’angolo con cui i due piani si incontrano.

[x,y]=meshgrid(-2:.1:2,-1:.1:1);

z1=1-x-y; z2=(1-x+2*y)/3;

surf(x,y,z1); hold on;

25

Page 26: Disp Matlab

26 CAPITOLO 3. SUPERFICI

surf(x,y,z2); hold off;

view(110,20);

n1=[1 1 1]; n2=[1 -2 3];

L1=sqrt(sum(n1.*n1)); L2=sqrt(sum(n2.*n2));

angle=acos(sum(n1.*n2)/L1/L2);

text(3,-1,-3,[’Angolo tra i due piani= ’,num2str(angle)]);

Esercizio 3.1 1. Disegnare il piano passante per (1, 3, 2) determinato dai due vettori

(1, 2, 1/2), (3,−1, 2).

2. Disegnare il piano determinato dall’equazione x+ y − z = 2.

3. Disegnare nello stesso grafico i piani x + y + z = 2 e 3x − 4y + 5z = −8, calcolandol’angolo con cui tali piani si incontrano.

3.1 Superfici cartesiane

Tra le superfici possiamo ovviamente considerare le superfici cartesiane, cioe le superfici chesono grafici di funzioni di due variabili. Richiamiamo quindi alcuni concetti gia visti nelcapitolo 2.

Disegniamo ad esempio il grafico della funzione f(x, y) = 2x2 − y2, (x, y) ∈ [−2, 2] ×[−4, 4]; di tale funzione ne considereremo anche le sezioni lungo i piani x = c, y = c e z = c,con c una costante fissata. Tracciamo anzitutto il grafico;

[x,y]=meshgrid(-2:.1:2,-4:.1:4);

surf(x,y,2*x.^2-y.^2);

shading interp;

Le sezioni con x = c o y = c si possono visualizzare semplicemente passando il parametro

shading faceted;

Le linee che si vedono altro non sono che le intersezioni del grafico di f con i piani verticalix = c e y = c. Per ottenere le intersezioni con i piani z = c bastera considerare il comandocountour3.

Il seguente esempio e dato da una superficie che e cartesiana ma come funzione y =f(x, z) = x2; tale superficie e un cilindro la cui base e descritta dalla parabola y = x2;

[x,z]=meshgrid(-2:.1:2,[-2,2]);

surf(x,x.^2,z);

Disegniamo infine il piano verticale x + 3y = 3, che e una superficie cartesiana se vistanella forma y = 1− x/3 oppure x = 3− 3y.

[y,z]=meshgrid(-1:2,-1:1);

x=3-2*y;

surf(x,y,z);

Esercizio 3.2 1. Disegnare il grafico di f(x, y) = (x2 − y2)2, tracciando anche i livellidi tale funzione.

2. Disegnare il grafico di f(x, y) = xe−x2−y2 ; affiancare a tale grafico il grafico costituito

da 25 curve di livello.

Page 27: Disp Matlab

3.2. SUPERFICI DI ROTAZIONE 27

3.2 Superfici di rotazione

L’esempio piu semplice di superficie di rotazione e dato dalla sfera; utilizziamo qui il comandosphere per disegnare la sfera x2 + y2 + z2 = 9.

[x,y,z]=sphere(30);

surf(3*x,3*y,3*z);

axis equal;

shading interp;

Il comando sphere crea la sfera di raggio 1 centrata nell’origine; possiamo quindi utilizzarela stessa funzione per disegnare un ellissoide;

[x,y,z]=sphere(30);

surf(4*x,3*y,2*z);

axis equal;

Possiamo utilizzare il comando cylinder per disegnare la superficie ottenuta ruotandouna funzione x = f(z) attorno all’asse z, con z ∈ [a, b]. La sintassi sara la seguente:

z=linspace(a,b,40);

x=f(z);

[X,Y,Z]=cylinder(x,30);

surf(X,Y,Z);

Per migliorarne l’effetto visivo, si puo usare il comando shading interp. Ma vediamo unesempio concreto; ruotiamo il grafico della funzione z2 + 1 attorno all’asse z. Il grafico chevogliamo ruotare sara il seguente;

z=linspace(1,3,50);

x=z.^2+1;

plot(x,z);

La rotazione si ottiene quindi nel seguente modo;

z=linspace(1,3,50);

x=z.^2+1;

[X,Y,Z]=cylinder(x,50);

figure(1);

surf(X,Y,Z);

zlabel(’z’); xlabel(’x’); ylabel(’y’);

figure(2);

surf(X,Y,Z);

shading interp;

Chiaramente si puo ottenere anche la rotazione della stessa funzione ruotandola attornoall’asse x invece che attorno all’asse z. Per far questo si ricava z in funzione di x e si ripetonoi comandi precedenti.

z=linspace(2,10,50);

x=sqrt(z-1);

[X,Y,Z]=cylinder(z,50);

surf(X,Y,Z);

zlabel(’z’); xlabel(’x’); ylabel(’y’);

Page 28: Disp Matlab

28 CAPITOLO 3. SUPERFICI

Si puo provare a sostituire il comando

surf(X,Y,Z);

col comando

surf(Z,Y,X);

e vedere cosa succede (si provi a commentare il risultato).Possiamo anche disegnare superfici di rotazione utilizzando le coordinate cilindriche e

sferiche, a seconda delle simmetrie della superficie che si vuole disegnare. Ad esempio, ilcono stesso z2 = a2(x2 +y2), a > 0 arbitrario, puo essere disegnato ponendo z = ±ar, r ≥ 0e ϑ ∈ [0, 2π] (di seguito consideriamo il caso a = 1/4).

a=.25;

[r,theta]=meshgrid(0:.2:2,0:2*pi/20:2*pi);

x=a*r.*cos(theta);

y=a*r.*sin(theta);

z=a*r;

surf(x,y,z); axis equal;

Analogamente, possiamo disegnare la superficie laterale di un cilindro come segue;

a=2;

[theta,z]=meshgrid(0:2*pi/30:2*pi,-3:6:3);

x=2*cos(theta); y=2*sin(theta);

surf(x,y,z); axis equal;

Per disegnare una sfera di raggio R possiamo invece utilizzare le coordinate sferiche, conangoli ϕ ∈ [0, π], ϑ ∈ [0, 2π];

R=2;

[theta,phi]=meshgrid(0:2*pi/30:2*pi,0:2*pi/30:pi);

x=R*cos(theta).*sin(phi);

y=R*sin(theta).*sin(phi);

z=R*cos(phi);

surf(x,y,z); axis equal;

ed in modo del tutto analogo si puo anche disegnare un’ellissoide; i seguenti comandidisegnano l’ellissoide di semiassi 5, 10 e 5/3 centrato in (1,−1, 0).

[theta,phi]=meshgrid(0:2*pi/30:2*pi,0:2*pi/30:pi);

x=1+5*cos(theta).*sin(phi);

y=-1+10*sin(theta).*sin(phi);

z=5/3*cos(phi);

surf(x,y,z); axis equal;

Esercizio 3.3 1. Disegnare il grafico di x2 + y2 = z2

2. Disegnare il grafico di x2 + y2

4 −z2

9 = 1.

3. Disegnare il grafico del cilindro z = sin(3x) con x ∈ [0, 2π], z ∈ [−1, 1].

Page 29: Disp Matlab

3.3. SUPERFICI PARAMETRIZZATE 29

4. Disegnare la superficie ottenuta ruotando la funzione

z = 1 +1

2sin2 x, x ∈ [0, π],

attorno all’asse x.

5. Disegnare la superficie ottenuta ruotando la funzione

z =1

x, x ∈ [0.1, 2],

attorno all’asse z.

6. Disegnare la parte di sfera di raggio 3 la cui latitudine e compresa tra [−π/3, π/3].

7. Disegnare il cono z2 = 19 (x2 + y2) e il cilindro y2 + z2 = 1 nello stesso grafico.

8. Disegnare la sfera di raggio 3 centrata in (1,−1, 0); usare hold on e quindi disegnarel’ellissoide (x− 1)2 + (y + 1)2/4 + 9z2 = 25 nello stesso grafico.

3.3 Superfici Parametrizzate

Disegneremo qui superfici definite parametricamente, cioe tramite funzioni

f(u, v) = (f1(u, v), f2(u, v), f3(u, v)),

dove f : [a, b]× [c, d]→ R3. Utilizzeremo i seguenti comandi

[u,v]=meshgrid(a:du:b,c:dv:d);

surf(f1(u,v),f2(u,v),f3(u,v));

dove tipicamente du e dv sono presi pari a .1 o .2; in alternativa, potremo usare anche ilcomando surfl a seconda dei casi (provare a capire le differenze tra surf e surfl utilizzandola guida di Matlab). Un esempio di superficie parametrizzata lo abbiamo gia incontratocome grafico di una funzione; vedremo anche esempi di superfici di rotazione.

3.3.1 Grafici di funzioni

Siamo qui nel caso gia visto di z = g(x, y), cioe con parametrizzazione data da f(u, v) =(u, v, g(u, v)), g : [a, b]× [c, d]→ R; avremo quindi

[x,y]=meshgrid(a:.1:b,c:.1:d);

surf(x,y,g(x,y));

Ad esempio, se consideriamo la funzione g(x, y) = x3 + y3 sul dominio (x, y) ∈ [−2, 2]2,avremo

[x,y]=meshgrid(-2:.1:2);

surfl(x,y,x.^3+y.^3);

con la seguente variante se si vuole modificare le coordinate della sorgente di luce:

Page 30: Disp Matlab

30 CAPITOLO 3. SUPERFICI

[x,y]=meshgrid(-2:.1:2);

S=[120,30];

surfl(x,y,x.^3+y.^3,S);

Si potrebbe aggiungere ai precedenti comandi alcuni comandi che disegnano il campo divettori normali alla superficie; ricordiamo che il campo normale e dato da

n(u, v) =fu(u, v)× fv(u, v)

‖fu(u, v)× fv(u, v)‖,

dove fu e fv rappresentano le derivate parziali della parametrizzazione f . Nel caso di ungrafico avremo quindi che fu(u, v) = (1, 0, gu(u, v)) e fv(u, v) = (0, 1, gv(u, v)), da cui

n(u, v) =1√

1 + |∇g(u, v)|2(−gu(u, v),−gv(u, v), 1).

Come esercizio, si provi ad utilizzare, come fatto nel capitolo 2, il comando quiver perdisegnare il campo normale alla superficie, cioe i vettori n(u, v).

3.3.2 Superfici di rotazione

Ricordiamo che una superficie di rotazione puo essere ottenuta ruotando attorno all’asse zil grafico di una funzione z = g(x), x ∈ [a, b], oppure ruotando una curva parametrizzata(x, z) = r(t), t ∈ [a, b]. La parametrizzazione della superficie si ottiene quindi mediantef : [a, b]× [0, 2π] data da, nel primo caso

f(u, v) = (u cos v, u sin v, g(u)),

mentre nel secondo caso

f(u, v) = (r1(u) cos v, r1(u) sin v, r2(v)).

Utilizzeremo quindi in generale i seguenti comandi:

[t,theta]=meshgrid(a:.1:b,0:2*pi/30:2*pi);

x=r1(t).*cos(theta);

y=r1(t).*sin(theta);

z=r2(t);

surf(x,y,z);

Possiamo ad esempio considerare il toro, ottenuto ruotando attorno all’asse z la cir-conferenza (x(t), z(t)) = (3 + cos t, sin t), t ∈ [0, 2π]. Fissiamo il punto di vista tramite ilcomando view e imponiamo che gli assi di riferimento siano uguali mediante axis equal.

[t,theta]=meshgrid(0:2*pi/30:2*pi);

x=(3+cos(t)).*cos(theta);

y=(3+cos(t)).*sin(theta);

z=sin(t);

surfl(x,y,z);

view(-10,20);

axis equal;

Page 31: Disp Matlab

3.3. SUPERFICI PARAMETRIZZATE 31

Come esercizio, si provi a dire quali sono, nel grafico appena ottenuto, le linee u 7→ f(u, v)con v costante e quali le linee v 7→ f(u, v) con u costante. Si usino inoltre i comandi shadinginterp e shading faceted.

Come esercizio, si utilizzi anche in questo caso il comando quiver per disegnare il camponormale alla superficie, cioe i vettori

n(u, v) =fu(u, v)× fv(u, v)

‖fu(u, v)× fv(u, v)‖,

che nel nostro caso, dato che f(u, v) = ((3 + cosu) cos v, (3 + cosu) sin v, sinu), sara dato da

n(u, v) =(sinu cos v,− sinu sin v, cosu)× (−(3 + cosu) sin v, (3 + cosu) cos v, 0)

‖(sinu cos v,− sinu sin v, cosu)× (−(3 + cosu) sin v, (3 + cosu) cos v, 0)‖=(− cosu cos v,− cosu sin v.− sinu).

Esercizio 3.4 1. Disegnare la superficie parametrizzata da

f(u, v) = (u sinu cos v, u cosu cos v, u sin v), (u, v) ∈ [0, 2π]2.

Si individuino le curve u = costante e v = costante.

2. Disegnare il grafico del cilindro circolare x2 +y2 = 1 di altezza 5 e si aggiustino gli assiin modo da farlo apparire circolare e si modifichi il punto di vista in modo da vederel’interno del cilindro.

3. Disegnare il nastro di Mobius parametrizzato da

f(u, v) = (2 cos v + u cosv

2, 2 sin v + u cos

v

2, u sin

v

2), u ∈ [−1/2, 1/2], v ∈ [0, 2π].

Utilizzare il comando rotate3d per aiutare la visualizzazione della superficie. Individ-uare, anche tracciando un secondo grafico, le curve con r = −1/2 e r = 1/2.

4. Disegnare il grafico della porzione di iperboloide −x2 − y2 + z2 = 1 individuata dax ∈ [−2, 2] e y ∈ [−3, 3].

5. Possiamo utilizzare le funzioni iperboliche per parametrizzare l’iperboloide del puntoprecedente; si disegni la superficie parametrizzata da

f(u, v) = (cosu sinh v, sinu sinh v,± cosh v), u ∈ [0, 2π], v ∈ [−2, 2].

6. Modificare il punto precedente per disegnare l’iperboloide x2 + y2 − z2 = 1.

Esercizio 3.5 Si rappresenti la proiezione stereografica della sfera, cioe si consideri laparametrizzazione della sfera f : R2 → R3 data da

f(u, v) =

(2u

u2 + v2 + 1,

2v

u2 + v2 + 1,u2 + v2 − 1

u2 + v2 + 1

)e nella stessa figura si disegni la sfera di raggio 1, il piano z = 0 e alcuni punti di tale piano(u, v, 0) insieme al segmento che li congiunge col punto di coordinate (0, 0, 1).

Page 32: Disp Matlab

32 CAPITOLO 3. SUPERFICI

3.4 Piano tangente ad una superficie parametrizzata

Ricordiamo che data una superficie parametrizzata da

f(u, v) = (x(u, v), y(u, v), z(u, v)), (u, v) ∈ [a, b]× [c, d],

l’equazione del piano tangente la superficie in f(u0, v0) e data da

(x(s, t), y(s, t), z(s, t) = f(u0, v0) + (s− u0)fu(u0, v0) + (t− v0)fv(u0, v0).

Prendiamo ad esempio l’iperboloide parametrizzato da

f(u, v) = (cosu sinh v, sinu sinh v, cosh v)

e il punto f(0, 1). Il piano tangente sara quindi dato da

(x(s, t), y(s, t), z(s, t) = (sinh 1 + t cosh 1, s sinh 1, cosh 1 + t sinh 1).

I seguenti comandi disegnano l’iperboloide per u ∈ [0, 2π], v ∈ [−2, 2] e il suo piano tangente.

[u,v]=meshgrid(0:pi/40:2*pi,-2:.1:2);

x=cos(u).*sinh(v); y=sin(u).*sinh(v); z=cosh(v);

surf(x,y,z); hold on;

[s,t]=meshgrid(-2:2:2);

X=sinh(1)+cosh(1)*t;

Y=sinh(1)*s;

Z=cosh(1)+sinh(1)*t;

surf(X,Y,Z); hold off;

Esercizio 3.6 1. Disegnare il grafico di z = sin(x2 + y2) e il suo piano tangente in(1, 1, sin 2).

2. Disegnare il grafico di

f(u, v) = (2 cosu sin v, sinu sin v, 3 cos v)

e il suo piano tangente nel punto f(0, π/3).

3. Disegnare il toro

f(u, v) = ((3 + cosu) cos v, (3 + cosu) sin v, sinu), (u, v) ∈ [0, 2π]2

e il suo piano tangente nel punto f(π/4, 0) (utilizzare anche il comando axis equal).

Page 33: Disp Matlab

Capitolo 4

Integrali multipli

In questo capitolo trattiamo gli integrali multipli, e piu precisamente gli integrali doppi nellasezione 4.1 e gli integrali tripli nella sezione 4.4.

4.1 Integrali doppi

Ricordiamo che l’integrale doppio di una funzione definita su di un rettangolo f : [a, b] ×[c, d]→ R e dato passando al limite nelle somme di Riemann∑

h,k

f(x∗h, y∗k)(xh+1 − xh, yk+1 − yk)

dove a ≤ x0 ≤ · · · ≤ xN ≤ b, c ≤ y0 ≤ · · · ≤ yM ≤ d induce una partizione del rettangoloin rettangolini [xh, xh+1] × [yk, yk+1] e (x∗h, y

∗k) ∈ [xh, xh+1] × [yk, yk+1] (sara il punto di

minimo di f su tale rettangolino se stiamo costruendo le somme integrali inferiori, mentresara il massimo di f sul rettangolino se stiamo considerando le somme integrali superiori).Il punto (x∗h, y

∗k) puo essere preso arbitrario nel caso in cui f sia una funzione continua (e

questo sara il caso per tutte le funzioni che qui considereremo).

Tale costruzione si adatta molto bene a Matlab, in quanto possiamo considerare par-tizioni del rettangolo fatte da rettangoli congruenti, cioe ad esempio

xh =h(b− a)

N, yk =

k(c− d)

M, h = 0, . . . , N, k = 0, . . . ,M.

Potremmo quindi considerare i seguenti comandi

dx=(b-a)/N; dy=(c-d)/M;

[x,y]=meshgrid(a+dx2:dx:b-dx/2,c+dy/2:dy:d-dy/2);

f=f(x,y);

Riemann=sum(sum(f))*dx*dy;

Ad esempio, se f(x, y) = sin(x + y) e il rettangolo e dato da [0, π/2]2, allora possiamoconsiderare i seguenti comandi, dove si disegna il grafico di sin(x+ y) e si calcola la sommadi Riemann;

33

Page 34: Disp Matlab

34 CAPITOLO 4. INTEGRALI MULTIPLI

[x,y]=meshgrid(0:pi/60:pi/2);

surf(x,y,sin(x+y));

N=50; M=50; dx=pi/100; dy=dx;

[x,y]=meshgrid(0+dx/2:dx:pi/2-dx/2);

f=sin(x+y); R=sum(sum(f))*dx*dy;

R=num2str(R); text(.1,.1,1.4,[’Somma di Riemann= ’,R]);

Si utilizzi il comando rotate3d per visualizzare meglio la figura. Si provi poi a ripeterequanto fatto ora con la funzione f(x, y) = cos2(x+ 2y) sul rettangolo [0, π]× [0, π/2].

4.2 Formule di riduzione per integrali doppi

Vediamo qui una applicazione della formula di riduzione di un integrale doppio con f :[a, b]× [c, d]→ R ∫

[a,b]×[c,d]

f(x, y)dxdy =

∫ b

a

dx

∫ d

c

f(x, y)dy.

In Matlab esiste la funzione dblquad che da una quadratura numerica dell’integrale doppio.Vediamo come funziona con un esempio; la funzione f deve essere data in pasto al comandodblquad o tra il comando fcnchk oppure deve essere definita utilizzando un m–file. Inpratica, f deve essere riconosciuta da Matlab come una function, che puo essere definitasia direttamente col comando fcnchk, oppure esternamente tramite un m–file. Calcoliamoad esempio l’integrale di

f(x, y) =1√

x2 + 2y2 + 1

con (x, y) ∈ [0, 1]× [0, 2]; useremo il comando fcnchk.

Stima=dblquad(fcnchk(’1./sqrt(x.^2+2*y.^2+1)’),0,1,0,2);

Stm=num2str(Stima);

[x,y]=meshgrid(0:.5:1,0:.1:2);

z=1./sqrt(x.^2+2*y.^2+1); surf(x,y,z);

rotate3d;

axis([0 1 0 2 0 1]);

text(.2,1.5,.2,[’Integrale= ’,Stm]);

Si puo anche usare il comando Matlab del symbolic toolbox int; tale comando serve percalcolare gli integrali di funzioni di una variabile, quindi per il calcolo di un integrale doppioandra utilizzato due volte.

syms x y real

Int1=int(1/sqrt(x^2+2*y^2+1),x,0,1);

R=int(Int1,y,0,2);

Si noti che con tali comandi il risultato esplicito non si riesce ad ottenere (il compilatoreo restituisce un’espressione complicata o scrive un messaggio di attenzione comunicandoche l’integrale non si riesce a calcolare esplicitamente). Possiamo pero ottenere il risultatonumerico dell’integrazione mediante il comando;

R=double(R);

Page 35: Disp Matlab

4.3. INTEGRALI DOPPI SU DOMINI GENERICI 35

Esercizio 4.1 1. Si utilizzi sia l’integrale iterato che il comando dblquad per calcolare∫[0,π/6]×[0,π/3]

x sin(x+ y)dxdy.

Si tracci anche il grafico della funzione integranda e mediante il comando text si riportiil valore delle integrazioni effettuate.

2. Trovare il volume del solido sottostante il grafico dell’iperboloide

z = x2 − y2, (x, y) ∈ [−1, 1]× [1, 3].

Si disegni il grafico della funzione f(x, y) = x2 − y2 e si riportino i valori delleintegrazioni effettuate sia con il dblquad che con l’integrazione iterata.

4.3 Integrali doppi su domini generici

Per poter integrare funzioni su domini limitati generici, procederemo in analogia a quanto sie fatto nel corso della teoria; piu precisamente, si considerera un rettangolo che contiene ilnostro dominio ed estenderemo la funzione ponendola zero nella regione del rettangolo chenon appartiene al dominio.

Prendiamo ad esempio la funzione f(x, y) = 3−x2−2y2 sul dominio E = {x2 +y2 ≤ 1};l’integrale di f su E rappresenta il volume del sottografico di f che e un cilindro compresotra il piano z = 0 e il paraboloide z = f(x, y). Calcoleremo l’integrale utilizzando le sommedi Riemann e anche il comando dblquad.

z=fcnchk(’(x.^2+y.^2<=1).*(3-x.^2-2*y.^2)’);

dx=2/50; dy=2/50;

[x,y]=meshgrid(-1+dx/2:dx:1-dx/2);

Riemann=sum(sum(z(x,y)))*dx*dy;

R=num2str(Riemann);

Stima=dblquad(z,-1,1,-1,1);

S=num2str(Stima);

[x,y]=meshgrid(-1:.2:1);

surf(x,y,z(x,y));

text(-1,.8,3.2,[’Somma di Riemann= ’,R]);

text(-1,.8,3.6,[’Integrazione numerica= ’,S]);

text(-1,.8,4,[’Valore preciso=9*pi/4=7.0686’]);

Una variazione ai comandi precedenti si puo effettuare quando si richiama la funzionedblquad: piu precisamente, la funzione puo essere definita tramite un m–file. Cosı adesempio potremmo considerare il comando

Stima=dblquad(z,-1,1,-1,1,Puntomedio);

dove Puntomedio.m e un m–file descritto ad esempio dai seguenti comandi;

function [A,n]=Puntomedio(fun,a,b,E);

fun=fcnchk(fun,’vectorized’);

x=a:.001:b;

y=fun(x);

Page 36: Disp Matlab

36 CAPITOLO 4. INTEGRALI MULTIPLI

y2=diff(y,2)/.001^2;

K=max(abs(y2));

n=ceil(sqrt(K*(b-a)^3/E/24));

h=(b-a)/n;

xstart=a+h/2:h:b-h/2;

f=fun(xstart);

A=sum(f)*h;

Atri modi per passare la funzione da integrare al comando dblquad possono essere i seguenti;

dblquad(@(x,y) funzione,a,b,c,d);

dove al posto della stringa funzione va inserita la definizione della nostra funzione e a, b, ce d sono gli estremi di integrazione, cioe nel nostro caso

funzione= 3-x.^2-2*y.^2

a=-1; b=1; c=-1; d=1;

In alternativa si puo considerare la sintassi:

dblquad(@funzione.m,a,b,c,d);

dove funzione.m e un m–file che definisce la funzione da integrare, e deve quindi esserestrutturato come segue;

function z=funzione(x,y);

...

4.3.1 Disegno di insiemi nel piano

A volte puo essere utile, prima di calcolare un integrale, avere un’idea della regione su cuisi vuole calcolare l’integrale. Ad esempio, se si chiede di integrare f(x, y) = xy sul dominiolimitato dai grafici y = x4 e y = 3x − x3, possiamo preliminarmente utilizzare i seguenticomandi;

x=linspace(-2,2);

plot(x,x.^4,x,3*x-x.^3);

Si nota quindi che i due grafici si incontrano per x = 0 e in un secondo punto x = b compresotra 1 e 1.5; avremo che l’integrale sara quindi dato da∫ b

0

∫ 3x−x3

x4

xydy =1

2

(b8

8+

9b4

4− b10

10− b6

).

Il punto b andrebbe determinato come soluzione dell’equazione

b3 + b2 − 3 = 0,

ma per questo tipo di equazioni non esistono formule risolutive. Con i seguenti comandi,individueremo tale soluzione numericamente con il comando fzero;

Page 37: Disp Matlab

4.4. INTEGRALI TRIPLI 37

x=linspace(-2,2); plot(x,x.^4,x,3*x-x.^3);

axis([-.25 1.5 -1 4]); grid on;

b=fzero(’x.^4-3*x+x.^3’,1);

syms x y real

syms x y real;

A=int(int(x*y,y,x^4,3*x-x^3),x,0,b);

A=num2str(double(A));

T=[’\int_D xy dx dy= ’,A];

text(.01,3.25,T); text(.7,1,’D’);

set(gca,’XTick’,[0:.5:1 b]);

text(-.1,3*(-.1)-(-.1)^3,’\leftarrow 3x-x^3’);

text(1,1,’\leftarrow x^4’);

Qui i comandi nell’ambiente text che vengono preceduti da un backslash sono riconosciutidal formattatore LaTeX, che e un formattatore per testi scientifici, molto usato in ambitomatematico. Si provi ad esempio ad inserire dentro un comando text i seguenti comandi

\frac{x^2}{\sqrt{1-x}}

\int_0^2 e^\sqrt{5-6x^2} dx

Esercizio 4.2 Si consideri l’integrale∫D

√x3 + 1dxdy

con D dominio limitato dai grafici y = x2, x = 1 e l’asse delle x.

1. Calcolare le somme di Riemann con N = M = 50 ed usare il comando fcnchk perestendere f(x, y) =

√x3 + 1 nulla all’esterno di D.

2. Disegnare il grafico di f , estesa a zero fuori da D, sul quadrato (x, y) ∈ [0, 1]2 ed usareil comando text per riportare il valore della somma di Riemann.

3. Calcolare a mano l’integrale.

4. Usare il comando dblquad per stimare l’integrale. Se il comando non funziona a causadelle troppe discontinuita della funzione, utilizzare come metodo alternativo l’m–filePuntomedio.m o creare un m–file alternativo per il calcolo dell’integrale (ad esempioscrivere un m–file che implementi il metodo di Simpson per il calcolo dell’integrale;per informazioni sul metodo di Simpson, rimandiamo all’appendice B).

5. Tracciare il grafico delle curve che definiscono il dominio D ed usare il comando textper mettere le etichette sulle curve, sul dominio D e per riportare il valore delle sommedi Riemann ottenute, del valore esatto calcolato a mano e il valore numerico ottenutocon dblquad (o con un qualsiasi altro metodo adottato per il calcolo).

4.4 Integrali tripli

La somma di Riemann per un integrale triplo sara definita mediante una griglia tridimen-sionale, cosı possiamo considerare comandi come segue;

Page 38: Disp Matlab

38 CAPITOLO 4. INTEGRALI MULTIPLI

dx=(b-a)/N; dy=(c-d)/M; z=(e-f)/P;

[x,y,z]=meshgrid(a:dx:b,c:dy:d,e:dz:f);

dove si e preparata la griglia per il calcolo di un integrale triplo su di un parallelepipedo[a, b]× [c, d]× [e, f ].

Se ad esempio vogliamo calcolare l’integrale della funzione f(x, y, z) = x2 + yz su E =[0, 2] × [−3, 0] × [−1, 1] con N = M = P = 30, possiamo considerare i seguenti comandi,dove a fianco della somma di Riemann utilizziamo anche l’integrazione iterata, cioe il fattoche E e un insieme semplice e∫

E

f(x, y, z)dxdydz =

∫ 2

0

dx

∫ 0

−3

dy

∫ 1

−1

(x2 + yz)dz.

dx=2/30; dy=3/30; dz=2/30;

px=dx/2:dx:2-dx/2;

py=-3+dy/2:dy:-dy/2;

pz=-1+dz/2:dz:1-dz/2;

[x,y,z]=meshgrid(px,py,pz);

Riemann=sum(sum(sum(x.^2+y.*z)))*dx*dy*dz

syms x y z real

A=int(int(int(x^2+y*z,z,-1,1),y,-3,0),x,0,2)

4.4.1 Integrali su domini generici

Consideriamo ora il problema dell’integrazione della funzione f(x, y, z) = ex sull’insieme

E = {(x, y, z) : 0 ≤ y ≤ 1, 0 ≤ x ≤ y, 0 ≤ z ≤ x+ y};

possiamo considerare i comandi seguenti;

w=fcnchk(’(0<=y).*(y<=1).*(0<=x).*(x<=y).*(0<=z).*(z<=x+y).*exp(x)’);

dx=1/50; dy=dx; dz=2/100;

[x,y,z]=meshgrid(dx/2:dx:1-dx/2,dy/2:dy:1-dy/2,dz/2:dz:2-dz/2);

Riemann=sum(sum(sum(w(x,y,z))))*dx*dy*dz;

syms x y z real;

A=int(int(int(exp(x),z,0,x+y),0,y),0,1)

Si noti come sono stati utilizzati gli operatori logici di confronto. Infine, si noti che il risultatoviene riportato scritto in forma esponenziale decimale; per convertirlo in un numero decimalebasta inserire la seguente riga di comando;

double(A);

4.4.2 Grafico di regioni solide

Puo essere utile prima di calcolare un integrale triplo a mano avere un’idea di come e fatto ildominio di integrazione, in modo da avere un aiuto per la scelta del metodo di integrazione.Possiamo quindi utilizzare Matlab per questo scopo, disegnando i bordi dell’insieme diintegrazione. Se ad esempio volessimo disegnare la regione solida della sezione precedentedescritta analiticamente da

E = {(x, y, z) : 0 ≤ y ≤ 1, 0 ≤ x ≤ y, 0 ≤ z ≤ x+ y},

potremmo procedere come segue;

Page 39: Disp Matlab

4.4. INTEGRALI TRIPLI 39

y=linspace(0,1); x=0*y; z=0*y;

plot3(x,y,z);

hold on; grid on;

x=0*y; z=y; plot3(x,y,z);

x=y; z=0*y; plot3(x,y,z);

x=y; z=x+y; plot3(x,y,z);

x=linspace(0,1); y=1+0*x; z=0*x; plot3(x,y,z);

z=x+1; plot3(x,y,z);

z=linspace(0,2); x=1+0*z; y=1+0*z;

plot3(x,y,z);

z=linspace(1,0); x=0*z; y=1+0*z;

plot3(x,y,z); hold off;

Esercizio 4.3 1. Si scriva l’integrale triplo∫E

e−%3

dxdydz

con E = {(x, y, z) : x2 + y2 + z2 ≤ 1} in coordinate sferiche; qui si intende % =√x2 + y2 + z2. Si scrivano le somme di Riemann per tale integrale triplo con N =

M = P = 40 e si riporti il risultato ottenuto in un grafico dove viene rappresentata lasfera di raggio 1. Si riporti inoltre il risultato del valore esatto dell’integrale ottenutoa mano o con il calcolo simbolico di Matlab.

2. Si calcolino le somme di Riemann per l’integrale∫E

ydxdydz

con E = {0 ≤ z ≤ x+ 2y, 0 ≤ x ≤ 1, 0 ≤ y ≤ x2}. Si disegni la regione solida E.

Page 40: Disp Matlab

40 CAPITOLO 4. INTEGRALI MULTIPLI

Page 41: Disp Matlab

Capitolo 5

Punti stazionari e massimi eminimi

5.1 Classificazione di punti stazionari

Il problema della classificazione dei punti stazionari consiste nella classificazione di unaforma quadratica; il primo metodo per la classificazione di una forma quadratica consistenello studio del segno degli autovalori di una matrice. Possiamo quindi utilizzare il comandoeig che, data una matrice, restituisce come risultato due matrici, una che contiene, nellecolonne, gli autovettori, e una seconda in forma diagonale sulla cui diagonale sono contenutigli autovalori. Cosı ad esempio possiamo considerare i seguenti comandi;

A=[-6 -6; -6 2];

[V,D]=eig(A)

Si deduce quindi che A ha un autovalore negativo e uno positivo, quindi A e una matriceindefinita. Analogamente;

A=[6 0 -2; 0 4 0; -2 0 2];

[V,D]=eig(A)

da cui si deduce che A e definita positiva.La verifica si puo anche fare applicando il Teorema di Sylvestre, e cioe calcolando i

determinanti dei minori principali; possiamo quindi considerare i seguenti comandi

A=[-6 6; -6 2];

A1(1,1)

det(A)

da cui si vede che i numeri ottenuti sono entrambi negativi per la prima matrice, e quindila matrice e indefinita, mentre

A=[6 0 -2; 0 4 0; -2 0 2];

A(1,1)

det([6 0; 0 4])

det(A)

41

Page 42: Disp Matlab

42 CAPITOLO 5. PUNTI STAZIONARI E MASSIMI E MINIMI

da cui si vede che tutti i determinanti dei minori principali sono positivi, e quindi A edefinita positiva.

Tipicamente il problema consiste nel cercare prima i punti stazionari e poi nel classificarlie per questo pu essere utilizzato il symbolic toolbox. Consideriamo per semplicit la funzionedi due variabili

f(x, y) = 4x2 + 3y2 + 12xy + 7x+ 8y + 5

che, come si vedr, ha un solo punto stazionario.

syms x y real

f = 4*x^2+3*y^2+12*x*y+7*x+8*y+5;

e calcoliamone i punti stazionari attraverso il comando solve applicato alle derivate parziali

fx=diff(f,x); fy=diff(f,y);

S=solve(fx,fy,x,y)

Si ottiene una struttura S che contiene un valore di x e uno di y, cio f(x, y) ha un solopunto stazionario (xm, ym).

xm=S.x

ym=S.y

Calcoleremo ora le derivate seconde e le valuteremo nel punto (xm, ym) con il comando subs,ottenendo la matrice Hessiana H(xm, ym).

fxx=diff(fx,x); fxy=diff(fx,y);

fyx=diff(fy,x); fyy=diff(fy,y);

H=[subs(fxx,{x,y},{xm,ym}) subs(fxy,{x,y},{xm,ym});

subs(fyx,{x,y},{xm,ym}) subs(fyy,{x,y},{xm,ym})]

Per determinare se (xm, ym) un punto di minimo, di massimo o di sella, troviamo gliautovalori di H;

eig(H)

Gli autovalori sono di segno opposto, dunque la matrice Hessiana indefinita e (xm, ym)un punto di sella. Andiamo a controllare nel grafico:

x=linspace(-1,.5);

[X,Y]=meshgrid(x);

f = 4.*X.^2+3.*Y.^2+12.*X.*Y+7.*X+8.*Y+5;

surf(X,Y,f)

5.2 Massimo e minimo di una funzione

I massimi e i minimi, assieme ai loro punti di massimo e di minimo di una funzione possonoessere individuati utilizzando i comandi Matlab min e max; l’accuratezza del calcolo deimassimi e dei minimi dipende dalla griglia scelta per definire la funzione.

Iniziamo con un esempio di una funzione di una variabile; supponiamo di voler deter-minare il massimo e il minimo della funzione

f(x) = 2 cos(2x)− 3 cosx, x ∈ [0, 2π].

Possiamo considerare molto semplicemente i seguenti comandi;

Page 43: Disp Matlab

5.2. MASSIMO E MINIMO DI UNA FUNZIONE 43

x=linspace(0,2*pi,700);

y=2*cos(2*x)-3*cos(x);

M=max(y)

m=min(y)

Si provi a calcolare analiticamente il massimo e il minimo della funzione data e si commentinoi risultati ottenuti.

Esercizio 5.1 Si scrivano comandi Matlab che calcolino il massimo e il minimo dellafunzione

f(x) = 2 sin(2x)− 3 cos(x

2),

ed utilizzino subplot per disegnare il grafico della funzione e della sua derivata (calcolatanumericamente utilizzando incrementi pari a .001) nell’intervallo [0, 2π]. Come variante diquesto esercizio, si ripeta il procedimento utilizzando il comando subplot per disegnare lafunzione assieme alle sue derivate prima e seconda calcolata analiticamente.

5.2.1 Massimo e minimo di una funzione di due variabili su di unrettangolo

Supponiamo ora di voler determinare il massimo e il minimo di

f(x, y) = sinx+ sin y + sin(x+ y), (x, y) ∈ [−π/4, π/4]2.

Possiamo quindi considerare i seguenti comandi, in cui la determinazione del massimo edel minimo viene fatta sia guardando il grafico della funzione, sia studiando le sue linee dilivello.

[x,y]=meshgrid(-pi/4:pi/80:pi/4);

z=sin(x)+sin(y)+sin(x+y);

figure(1); surf(x,y,z);

figure(2); [C,h]=contour(x,y,z,20);

clabel(C,h);

Studiando i due grafici ottenuti, si deduce subito che il massimo viene assunto nel punto(π/4, π/4), mentre il punto di minimo viene assunto in (−π/4,−π/4); possiamo modificarei precedenti comandi in modo da individuare i valori ed i punti di massimo e minimo dellafunzione;

[x,y]=meshgrid(-pi/4:pi/80:pi/4);

z=sin(x)+sin(y)+sin(x+y);

figure(1); surf(x,y,z);

[M,ind]=max(z(1:prod(size(z))));

Mx=num2str(x(ind)); My=num2str(y(ind));

M=num2str(M);

TM=[’Massimo= ’,M,’, assunto in (’,Mx,’,’,My,’)’];

text(-1,1,5,TM);

Si provi a modificare i precedenti comandi in modo da individuare anche il minimo e il puntodi minimo. Si provi quindi ad usare il comando help per capire come funzionano i comandimax e min, sia nel caso in cui l’argomento e un vettore, che nel caso in cui l’argomento e una

Page 44: Disp Matlab

44 CAPITOLO 5. PUNTI STAZIONARI E MASSIMI E MINIMI

matrice. Se ne dovrebbe dedurre che nel caso in cui z e una matrice, max(z) e un vettoreriga, i cui elementi sono determinati dai massimi di z in ogni colonna. In questo modo,il valore massimo e individuato semplicemente con il comando max(max(z)). Nei comandiconsiderati sopra, quello che e stato fatto e stato praticamente ridurre la matrice z ad unvettore di lunghezza il prodotto delle dimensioni della matrice z; in questo modo fare maxuna volta sola coincide con il comando max(max(z)).

5.2.2 Massimi e minimi su cerchi

Supponiamo ora di voler calcolare il massimo e il minimo della funzione

f(x, y) = x2 + 2y2 − 4y − 2, x2 + y2 ≤ 4.

Ci sono sostanzialmente due modi di procedere; uno e l’utilizzo delle coordinate polari, ilsecondo e estendere la funzione su di un dominio rettangolare ponendola zero all’esterno delcerchio mediante gli operatori logici; in questo secondo metodo c’e pero da tener presenteche se la funzione e tutta strettamente positiva o strettamente negativa nel cerchio, allorala risposta ottenuta sarebbe incorretta, e quindi il metodo andrebbe modificato.

Il primo metodo viene trattato con i seguenti comandi;

[r,theta]=meshgrid(0:.1:2,0:pi/20:2*pi);

x=r.*cos(theta); y=r.*sin(theta);

f=x.^2+2*y.^2-4*y-2; surfl(x,y,f);

[m,imin]=min(f(1:prod(size(f))));

m=num2str(m);

mx=num2str(x(imin)); my=num2str(y(imin));

Tm=[’Minimo= ’,m,’, assunto in (’,mx,’,’,my,’)’];

text(0,2,10,Tm);

Nel secondo caso possiamo considerare i seguenti comandi:

[x,y]=meshgrid(-2:.05:2);

f=(x.^2+y.^2<=4).*(x.^2+2*y.^2-4*y-2);

surfl(x,y,f);

[M,imax]=max(f(1:prod(size(f))));

M=num2str(M);

Mx=num2str(x(imax)); My=num2str(y(imax));

TM=[’Massimo= ’,M,’, assunto in (’,Mx,’,’,My,’)’];

text(0,2,10,TM);

Provare ad impostare il punto di vista con view(-100,10) per avere una buona visualizzazionedei punti di massimo e di minimo.

Esercizio 5.2 1. Disegnare il grafico e le linee di livello della funzione

f(x, y) = 2x3 + xy2 + 5x2 + y2

nel rettangolo −2.5 ≤ x ≤ 0.6, −3 ≤ y ≤ 3; inserire a mano nel grafico i puntistazionari ed etichettarli come punti di sella o punti di massimo o minino locali.

Page 45: Disp Matlab

5.3. MASSIMI E MINIMI VINCOLATI 45

2. Utilizzare i metodi grafici (grafico e linee di livello) per individuare i punti estremalidi

f(x, y) = e−xy, 4x2 + y2 ≤ 2.

Utilizzare il comando text per riportare i valori massimo e minimo e i punti dove ven-gono assunti. Si ripeta l’esercizio utilizzando le coordinate polari (ellittiche). Si utilizzianche il metodo di estensione di f ad un dominio rettangolare mediante operatori logicie si commenti il risultato ottenuto.

5.3 Massimi e minimi vincolati

Quando il vincolo e parametrizzabile, possiamo individuare il massimo e il minimo pervia grafica; il grafico del vincolo e delle linee di livello possono inoltre essere utilizzati perl’individuazione dei punti estremali (saranno individuati dai punti di tangenza tra il vincoloe il livello della funzione), dando cosı una visualizzazione del metodo dei moltiplicatori diLagrange.

Supponiamo ad esempio di voler trovare gli estremi della funzione

f(x, y) = x3 + y3 + 3xy

soggetta al vincolo (x − 3)2 + (y − 3)2 = 9. Tracciamo quindi il grafico di f , la curva chedefinisce il vincolo e la curva nello spazio che si ottiene restringendo la funzione al vincolo:

[x,y]=meshgrid(-1:.1:7);

z=x.^3+y.^3+3*x.*y;

surf(x,y,z); hold on;

t=linspace(0,2*pi,300);

x=3+3*cos(t); y=3+3*sin(t);

f=x.^3+y.^3+3*x.*y;

plot3(x,y,f,’r’);

plot3(x,y,0*f);

hold off;

Si utilizzi il comando rotate3d per avere una visualizzazione migliore del risultato ottenuto.Utilizzeremo ora una procedura grafica che ci dara una buona stima dei valori estremali

di f soggetta al vincolo dato. Nel primo caso parametrizzeremo la curva e quindi tracceremoil grafico della funzione ristretta a tale curva; gli estremali saranno individuati guardandotale grafico. Il vincolo e una circonferenza con parametrizzazione x = 3 + 3 cos t,

y = 3 + 3 sin t, t ∈ [0, 2π].

Possiamo quindi utilizzare tale parametrizzazione per ottenere una funzione di una variabile,e possiamo quindi individuare il massimo e il minimo utilizzando i comandi max e min.

t=linspace(0,2*pi,300);

x=3+3*cos(t); y=3+3*sin(t);

f=x.^3+y.^3+3*x.*y; plot(t,f);

Page 46: Disp Matlab

46 CAPITOLO 5. PUNTI STAZIONARI E MASSIMI E MINIMI

Quello ottenuto e il grafico della funzione ristretta al vincolo; si potrebbe darne una visual-izzazione tridimensionale, considerando ad esempio il comando

plot3(x,y,f);

Come esercizio, si tracci il grafico cosı ottenuto, assieme alla curva che definisce il vincolovista nello spazio tridimensionale. Come variante e/o abbellimento della figura, invece delcomando plot3(x,y,f) si potrebbe utilizzare il comando stem3(x,y,f).

Osservando il grafico ottenuto, si intuisce che il massimo e circa 350 ed e ottenuto pert = 0.8, mentre il minimo e pari a 0 ottenuto circa per t = 3.95. Per avere una stimamigliore di tali risultati, si puo aggiungere la seguente lista di comandi a quanto fattoprecedentemente;

[m,mt]=min(f); tmin=num2str(t(mt));

m=num2str(m);

T=[’Minimo= ’,m,’ ottenuto per t= ’,tmin];

text(3.5,50,T);

Con il secondo metodo, disegneremo la curva che definisce il vincolo e i livelli di f . Ovvi-amente, le limitazioni sulla griglia scelta vanno scelte in modo da contenere tutta la curvache definisce il vincolo; anche la scelta dei livelli dovra essere fatta in maniera opportuna edin genere si procede per tentativi successivi.

[x,y]=meshgrid(-1:.1:7);

f=x.^3+y.^3+3*x.*y;

L=[0:10 20:10:100 150:50:350];

[C,h]=contour(x,y,f,L);

v=[0:2:10 150:50:350];

clabel(C,h,v); hold on;

t=linspace(0,2*pi,300);

x=3+3*cos(t); y=3+3*sin(t);

plot(x,y); hold off;

Come variante, per una visualizzazione del metodo dei moltiplicatori di Lagrange, si potreb-be disegnare il campo gradiente della funzione che definisce il vincolo

g(x, y) = (x− 3)2 + (y − 3)2 − 9

assieme al campo gradiente della funzione f ; tali campi gradiente vanno disegnati sulla curvag(x, y) = 0. I punti estremali si avranno in corrispondenza dei punti nei quali tali campigradiente sono paralleli tra loro.

Esercizio 5.3 Si utilizzi il metodo grafico per stimare il massimo e il minimo della funzione

f(x, y) = xy

soggetta al vincolo 4x2 + y2 = 9. Inoltre, si tracci nello stesso disegno il grafico di f , lacurva che definisce il vincolo e la curva dello spazio ottenuta restringendo f al vincolo. Siparametrizzi inoltre il vincolo e si tracci il grafico della funzione di una variabile ottenutarestringendo la funzione al vincolo e da tale grafico si dia una stima del massimo e delminimo della funzione. Si utilizzino quindi i comandi max e min per la determinazione deipunti estremali; usare infine il comando text per scrivere i valori ottenuti sul grafico.

In conclusione, si usi il metodo delle linee di livello per individuare i punti di massimo edi minimo, etichettando i livelli.

Page 47: Disp Matlab

Capitolo 6

Integrali curvilinei e disuperficie

In questa sezione considereremo gli integrali curvilinei di prima e seconda specie (rispetti-vamente, per funzioni scalari e vettoriali) e gli integrali di superficie per funzioni scalari evettoriali. Come applicazioni considereremo le formule di Gauss–Green e Stokes.

6.1 Integrali curvilinei

6.1.1 Integrali curvilinei per funzioni scalari

Una prima applicazione degli integrali curvilinei di prima specie, come visto nelle lezioni diteoria, e ad esempio il calcolo della massa di un filo la cui densita di massa sia descritta dauna data funzione. Una volta calcolata la massa totale del filo, e possibile calcolare quindile coordinate del baricentro del filo. Supponiamo ad esempio di avere un filo della formadi una semicirconferenza di raggio 2 nel piano la cui densita di massa sia descritta dallafunzione

f(s) =√

1 + s2;

qui stiamo supponendo che il parametro s sia il parametro d’arco calcolato a partire dalpunto (2, 0). Stiamo quindi considerando la parametrizzazione r : [0, 2π]→ R2,

r(s) = (2 cos(s/2), 2 sin(s/2)).

L’intervallo e la parametrizzazione sono stati calcolati tenendo conto che la lunghezza delfilo e L = 2π. La massa totale del filo sara quindi data da

M =

∫ 2π

0

√1 + s2ds.

Si noti che nella formula precedente non compare ‖r′(s)‖, in quanto, avendo scelto s comeparametro d’arco, tale quantita e pari ad 1. Calcoliamo il precedente integrale usando ilsymbolic toolbox;

syms s

47

Page 48: Disp Matlab

48 CAPITOLO 6. INTEGRALI CURVILINEI E DI SUPERFICIE

densita=sqrt(1+s^2);

int(densita,0,2*pi)

massa=double(ans)

Calcoliamo ora le coordinate del baricentro del filo; per effettuare questo calcolo, invecedi utilizzare il parametro d’arco, considereremo la paremtrizzazione equivalente data daγ : [0.π]→ R2,

γ(t) = (2 cos t, 2 sin t).

Dato che ‖γ′(t)‖ = 2, si avra che s(t) = 2t, da cui r(s) = r(2t) = γ(t) e la densita diventa√1 + 4t2. Quindi le coordinate del baricentro saranno date da

x =1

M

∫ π

0

2 cos t√

1 + 4t2dt, y =1

M

∫ π

0

2 sin t√

1 + 4t2dt.

Per il calcolo di tali coordinate utilizziamo il comando quad8, in quanto il calcolo simboliconon si riesce a fare.

f=inline(’4*sqrt(1+4*t.^2).*cos(2*t)’);

g=inline(’4*sqrt(1+4*t.^2).*sin(2*t)’);

xbar=(1/massa)*quad8(f,0,pi)

ybar=(1/massa)*quad8(g,0,pi)

6.1.2 Rappresentazione di campi vettoriali

Prima di affrontare il problema del calcole di integrali curvilinei per funzioni vettoriali,proviamo a dare una rappresentazione grafica dei campi di vettori che vogliamo integrare.

Campi nel piano Un campo di vettori nel piano e determinato da una coppia di funzioniscalari,

F (x, y) = (F1(x, y), F2(x, y))

e il metodo migliore per visualizzare un campo di vettori, come gia visto nei capitoliprecedenti, e l’utilizzo del comando quiver.

Ad esempio, supponiamo di voler disegnare il campo

F (x, y) = (1, x+ y2), (x, y) ∈ [−2, 3]× [−1, 2].

La cosa importante e non disegnare troppi vettori del campo in considerazione, altrimenti ildisegno rischia di diventare illeggibile. Una griglia divisa tra le 10 e le 15 parti e solitamentepiu che sufficiente.

f1=inline(’0*x+1’,’x’,’y’);

f2=inline(’x+y.^2’,’x’,’y’);

x=linspace(-2,3,11);

y=linspace(-1,2,11);

[X,Y]=meshgrid(x,y);

F1=f1(X,Y); F2=f2(X,Y);

quiver(X,Y,F1,F2)

axis image

Page 49: Disp Matlab

6.1. INTEGRALI CURVILINEI 49

Il comando axis image serve solamente i numeri sulle coordinate x e y (confrontare conil comando axis equal). Se vogliamo disegnare il grafico delle direzioni del campo F ,possiamo modificare i comandi precedenti, semplicemente normalizzando il vettore F (x, y).

f1=F1./sqrt(F1.^2+F2.^2);

f2=F2./sqrt(F1.^2+F2.^2);

quiver(X,Y,f1,f2)

axis image

Potremmo anche disegnare il campo gradiente di una funzione assegnata; se prendiamoad esempio la funzione

f(x, y) = e−x2

2 −y2−xy, (x, y) ∈ [−2, 2];

[x,y]=meshgrid(-2:.2:2);

f=exp(-x.^2/2-y.^2+x.*y);

[fx,fy]=gradient(f,.2,.2);

quiver(x,y,fx,fy);

Campi nello spazio In modo analogo possiamo disegnare campi vettoriali nello spazio;ad esempio, consideriamo il campo

F (x, y, z) = (1, x+ y2, z), (x, y, z) ∈ [−2, 3]× [−1, 2]× [−1, 1].

u=inline(’1+0*x’,’x’,’y’,’z’);

v=inline(’x+y.^2’,’x’,’y’,’z’);

w=inline(’z’,’x’,’y’,’z’);

x=linspace(-2,3,6);

y=linspace(-1,2,6);

[X,Y]=meshgrid(x,y);

for z = -1:.4:1

Z=z+0*X;

V=v(X,Y,Z);

W=w(X,Y,Z);

quiver3(X,Y,Z,U,V,W)

hold on

end

hold off

6.1.3 Integrali curvilinei per funzioni vettoriali

Ricordiamo che data una funzione F : D ⊂ R3 → R3 e una curva r : [a, b]→ D, si definiscel’integrale curvilineo di seconda specie come:∫

r

Fd~r :=

∫ b

a

F (r(t)) · r′(t)dt.

Cosı ad esempio se consideriamo la funzione F (x, y, z) = (x, z, ex+y) e la curva

r(t) = (t cos t, t sin t, t2), t ∈ [0, 2π],

possiamo calcolare l’integrale curvilineo come segue;

Page 50: Disp Matlab

50 CAPITOLO 6. INTEGRALI CURVILINEI E DI SUPERFICIE

u=inline(’x’,’x’,’y’,’z’);

v=inline(’z’,’x’,’y’,’z’);

w=inline(’exp(x+y)’,’x’,’y’,’z’);

n=200;

t=linspace(0,2*pi,n+1); dt=2*pi/n;

s=simpvec(n);

x=t.*cos(t); y=t.*sin(t); z=t.^2;

xdot=cos(t)-t.*sin(t);

ydot=sin(t)+t.*cos(t);

zdot=2*t;

I1=u(x,y,z).*xdot;

I2=v(x,y,z).*ydot;

I3=w(x,y,z).*zdot;

integrale=dot(s,(I1+I2+I3))*dt/3

Qui il comando simpvec e una functions che possiamo definire mediante un m–file comesegue;

function s=simpvec(n)

s=2*ones(1,n+1);

s(2:2:n)=4*ones(1,n/2);

s(1)=1; s(n+1)=1;

Puo succedere che la curva su cui calcolare il nostro integrale curvilineo non sia definitain modo parametrico ma sia semplicemente descritta da una tabella di valori. Ad sempiopotremmo avere la seguente tabella di valori;

x 0 .1 .25 .4 .54 .76 .82 .93 1y 0 .005 .0312 .0800 .1458 .2888 .3362 .4352 .5

Se la funzione da integrare e data da F (x, y) = (x cos y, x+y), per il calcolo dell’integralecurvilineo utilizziamo l’approssimazione che considera come curva la poligonale che congingei punti (Pj)j=1,...,n definiti dalla precedente tabella (quindi con n = 9):∫

r

F · d~r =

n∑j=1

∫ 1

0

F (rj(t)) · r′j(t)dt

con rj(t) = (1− t)Pj + tPj+1

x=[0 .1 .25 .4 .54 .76 .82 .93 1];

y=[0 .005 .0312 .0800 .1458 .2888 .3362 .4352 .5];

u=inline(’x.*cos(y)’,’x’,’y’);

v=inline(’x+y’,’x’,’y’);

n=50;

s=simpvec(n);

dt=1/n;

lavoro=0;

for j=1:8

xx=linspace(x(j),x(j+1),n+1);

yy=linspace(y(j),y(j+1),n+1);

Page 51: Disp Matlab

6.2. AREE DI SUPERFICI PARAMETRIZZATE 51

I=(x(j+1)-x(j))*u(xx,yy)+(y(j+1)-y(j))*v(xx,yy);

integrale=dot(s,I)*dt/3;

lavoro=lavoro+integrale;

end

lavoro

6.2 Aree di superfici parametrizzate

Ricordiamo che la formula per il calcolo dell’area di una superficie Σ = r(D) parametrizzatada r : D ⊂ R2 → R3 e data dalla formula

Area(Σ) =

∫D

‖ru(u, v)× rv(u, v)‖dudv;

nel caso di una superficie cartesiana, r(u, v) = (u, v, f(u, v)), tale formula diventa

Area(Σ) =

∫D

√1 + ‖∇f(u, v)‖2dudv.

Se volessimo ad esempio calcolare l’area di un elicoide

r(u, v) = (u cos(2v), u sin(2v), v), (u, v) ∈ [0, 1]× [0, π],

l’area da calcolare diventa

Area(Σ) =

∫ π

0

dv

∫ 1

0

√1 + 4u2du.

Per il calcolo di tale area utilizzeremo il symbolic toolbox come segue;

[u,v]=meshgrid(0:.05:1,0:pi/20:pi);

surf(u.*cos(2*v),u.*sin(2*v),v);

du=1/50; dv=pi/50;

[u,v]=meshgrid(du/2:du:1-du/2,dv/2:dv:pi-dv/2);

Riemann=sum(sum(sqrt(1+4*u.^2)))*du*dv;

T2=[’Somma di Riemann= ’,num2str(Riemann)];

text(-.9,1,3.4,T2);

Stima=pi*quad(fcnchk(’sqrt(1+4*u)’),0,1);

T3=[’Area stimata= ’,num2str(Stima)];

text(-.9,1,3,T3);

syms u real

Area=pi*int(sqrt(1+4*u^2),u,0,1);

T1=[’Area= ’,num2str(double(Area))];

text(-.9,1,3.8,T1)

Se invece volessimo calcolare l’area della superficie cartesiana data dal grafico dellafunzione

f(x, y) = x2 + y2, (x, y) ∈ [0, 1]2,

l’integrale da calcolare diventa

Area(Σ) =

∫[0,1]2

√1 + 4x2 + 4y2dxdy.

Possiamo considerare i seguenti comandi;

Page 52: Disp Matlab

52 CAPITOLO 6. INTEGRALI CURVILINEI E DI SUPERFICIE

[x,y]=meshgrid(0:.05:1);

surfl(x,y,x.^2+y.^2);

dx=1/50; dy=1/50;

[x,y]=meshgrid(dx/2:dx:1-dx/2);

Riemann=sum(sum(sqrt(1+4*x.^2+4*y.^2)))*dx*dy;

T1=[’Somma di Riemann= ’,num2str(Riemann)];

text(0,1,2.3,T1);

Stima=dblquad(fcnchk(’sqrt(1+4*x+4*y)’),0,1,0,1);

T2=[’Area mediante dblquad= ’,num2str(Stima)];

text(0,1,2.1,T2);

Esercizio 6.1 1. Calcolare l’area della superficie laterale del toro, ottenuto ruotando lacirconferenza di raggio 2 centrata in (3, 0). Si disegni il grafico della superficie e siriportino i seguenti valori;

(a) l’area calcolata mediante l’integrale e il valore ottenuto calcolandola esplicita-mente a mano;

(b) la somma di Riemann con m = n = 50 usando il metodo del punto medio;

(c) una stima dell’area ottenuta mediante dblquad.

2. Si trovi l’area della porzione di sfera di raggio 1 con angolo ϕ ∈ [0, π/4] e si riportinosul grafico della superficie i seguenti valori:

(a) l’area ottenuta calcolando a mano l’integrale;

(b) una stima dell’area ottenuta mediante integrazione numerica;

(c) la somma di Riemann ottenuta con n = m = 50 mediante il metodo del puntomedio.

3. Si ripeta l’esercizio precedente con la porzione di sfera che giace sopra il quadrato[−1/

√2, 1/√

2]2.

4. Si disegni il grafico della superficie definita parametricamente in coordinate sfericheda

%(ϑ, ϕ) = 1 +1

5sin(6ϑ) cos(5ϕ)

utilizzando i vari metodi visti sopra, in particolare dblquad e le somme di Riemann.

6.3 Integrali di superficie

Cosı come per gli integrali curvilinei, anche per gli integrali di superficie si distinguono duecasi, il caso scalare e il caso vettoriale.

6.3.1 Integrali di superficie per funzioni scalari

Se Σ e una superficie regolare parametrizzata da una funzione r : D ⊂ R2 → R3, si definiscel’integrale di una funzione continua f : Σ→ R ponendo∫

Σ

fdΣ :=

∫D

f(r(u, v))‖ru(u, v)× rv(u, v)‖dudv.

Page 53: Disp Matlab

6.3. INTEGRALI DI SUPERFICIE 53

Se poi la superficie e cartesiana, cioe r(u, v) = (u, v, g(u, v)), la formula precedente diventa∫Σ

fdΣ :=

∫D

f(u, v, g(u, v))√

1 + ‖∇g(u, v)‖2dudv.

Quindi, se ad esempio volessimo esempio calcolare l’integrale di f(x, y, z) = x+ y + z sullasuperficie grafico della funzione g(x, y) = sin(x + y2), (x, y) ∈ [0, 2]2, possiamo considerarei seguenti comandi;

f=fcnchk(’x+y+sin(+y.^2).*sqrt(1+cos(x+y.^2).^2.*(1+4*y.^2))’,’x’,’y’);

dblquad(f,0,2,0,2)

Come esercizio, seguendo l’interpretazione di f come una densita di massa, si calcoli ilbaricentro della superficie appena considerata.

6.3.2 Rappresentazione di campi lungo una superficie

Tracciamo ora il grafico del campo vettoriale F (x, y, z) = (y, z, x) lungo la superficie x +y + z = 0. Useremo un meshgrid per definire [x,y] e ricaveremo la z di conseguenza;

[x,y]=meshgrid(-2:2:2);

surfl(x,y,-x-y); hold on

[x,y]=meshgrid(-2:.5:2);

z=-x-y;

quiver3(x,y,z,y,z,x);

hold off

Potremmo anche solo disegnare il campo senza disegnare la superficie,

[x,y]=meshgrid(-2:.5:2);

z=-x-y;

quiver3(x,y,z,y,z,x);

6.3.3 Integrali di superficie per campi vettoriali: circuitazioni eflussi

Rotore e Teorema di Stokes Una prima applicazione puo essere il calcolo della cir-cuitazione di un campo; come esempio considereremo il campo

F (x, y, z) = (−z, x, y)

e come curva considereremo l’intersezione del paraboloide z = x2 + y2 e il piano z = y. Talecurva puo essere parametrizzata da

r(t) =

(1

2cos t,

1

2+

1

2sin t,

1

2+

1

2sin t

), t ∈ [0, 2π].

Vogliamo quindi verificare che vale la seguente identita;∮r

F · d~r =

∫Σ

rotF · nΣdΣ,

Page 54: Disp Matlab

54 CAPITOLO 6. INTEGRALI CURVILINEI E DI SUPERFICIE

dove Σ e la porzione del piano z = y contenuta all’interno della curva r. Nella formula prece-dente il versore nΣ deve essere preso in maniera coerente, cioe in modo tale che l’orientazioneindotta da tale campo su Σ induca la giusta orientazione sulla curva r. Avendo fissato laparametrizzazione r, il versore nΣ, ortogonale al piano z = y, deve avere componento zpositiva, quindi

nΣ =1√2

(0,−1, 1),

Con i seguenti comandi, disegneremo la curva r, il piano z = y assieme al versore normalenΣ(0, 1/2, 1/2) e calcoliamo l’integrale curvilineo∫ 2π

0

(−1

2− 1

2sin t,

1

2cos t,

1

2+

1

2sin t

)·(−1

2sin t,

1

2cos t,

1

2cos t

)dt

e l’integrale di superficie

(6.1)

∫Σ

(1,−1, 1) · 1√2

(0,−1, 1)dΣ

con Σ = {(x, y, z) ∈ R3 : z = y, x2 + y2 ≤ y}. Con i seguenti comandi calcoliamo l’integralecurvilineo;

syms t real

Ft=[-1/2-1/2*sin(t) 1/2*cos(t) 1/2+1/2*sin(t)];

T=[-1/2*sin(t) 1/2*cos(t) 1/2*cos(t)];

fz_int=dot(Ft,T);

I=double(int(fz_int,0,2*pi))

Si scrivano come esercizio i comandi che servono a calcolare l’integrale di superficie (6.1).

Divergenza e flussi Chiudiamo questo capitolo con il calcolo di un flusso di un campovettoriale. Consideriamo quindi il campo

F (x, y, z) = (y, x, z√x2 + y2)

e come superficie Σ consideriamo il bordo del cilindro

E = {(x, y, z) ∈ R3 : x2 − 2x+ y2 ≤ 0, 0 ≤ z ≤ 1}.

Calcoleremo tale flusso in due modi differenti, verificando quindi la validita della formula diGauss–Green ∫

E

divF (x, y, z)dxdydz =

∫Σ

F · nΣdΣ,

con nΣ normale uscente dal cilindro. Se vogliamo calcolare l’integrale della divergenza,abbiamo che

divF (x, y, z) =√x2 + y2,

mentre il dominio di integrazione in coordinate cilindriche

x = 1 + % cosϑ, y = % sinϑ, z = t

con % ∈ [0, 1], ϑ ∈ [0, 2π) e t ∈ [0, 1]. Si tratta quindi di calcolare il seguente integrale;

(6.2)

∫E

divF (x, y, z)dxdydz =

∫E

√x2 + y2dxdydz =

∫ 2π

0

∫ 1

0

%√

1 + %2 + 2% cosϑdϑ,

che cacleremo utilizzando Matlab tramite i seguenti comandi.

Page 55: Disp Matlab

6.3. INTEGRALI DI SUPERFICIE 55

f=fcnchk(’sqrt(1+rho.^2+2*rho.*cos(theta)).*rho’,’rho’,’theta’);

dblquad(f,0,1,0,2*pi)

L’integrale di superficie va invece diviso in tre parti, tenendo conto la decomposizione

Σ = Σ1 ∪ Σ2 ∪ Σ3

doveΣ1 = {x2 − 2x+ y2 = 0, 0 ≤ z ≤ 1}, Σ2 = {z = 0, x2 − 2x+ y2 < 0},

Σ3 = {z = 1, x2 − 2x+ y2 < 0}.

La prima superficie la parametrizziamo in coordinate cilindriche r : [0, 2π)× [0, 1],

r(ϑ, t) = (1 + cosϑ, senϑ, t),

da cui si ricava chenΣ1

(ϑ, t) = (cosϑ, senϑ, 0).

Otteniamo quindi che∫Σ1

F · nΣ1dΣ =

∫ 2π

0

∫ 1

0

( senϑ, 1 + cosϑ,√

1 + 2 cosϑ) · (cosϑ, senϑ, 0)dt

=

∫ 2π

0

(2 senϑ cosϑ+ senϑ)dϑ = 0.

La superficie Σ2 e il grafico z = 0 con normale nΣ2(x, y, 0) = (0, 0,−1); dato che F (x, y, 0) =

(y, x, 0), si vede subito che ∫Σ2

f · nΣ2dΣ = 0.

Per l’ultimo integrale, vista Σ3 come grafico z = 1 con normale nΣ3(x, y, 1) = (0, 0, 1),troviamo l’integrale∫

Σ3

F · nΣ3dΣ =

∫{x2−2x+y2≤0}

F (x, y, 1) · (0, 0, 1)dxdy =

∫{x2−2x+y2≤0}

√x2 + y2dxdy,

che e lo stesso integrale in (6.2).

Page 56: Disp Matlab

56 CAPITOLO 6. INTEGRALI CURVILINEI E DI SUPERFICIE

Page 57: Disp Matlab

Capitolo 7

Successioni e serie di funzioni

In questo capitolo considereremo le successioni e le serie di funzioni, tra cui principalmentetratteremo le serie di potenze, le serie di Taylor e le serie di Fourier.

7.1 Successioni di funzioni

Cominciamo con una successione di funzione data da

fn(x) =nx

n2x2 + 1,

definite su tutto R. Possiamo considerare i seguenti comandi, che servono per disegnare igrafici delle prime dieci funzioni della successione;

x=linspace(-.5,.5,100);

[X,N]=meshgrid(x,n);

fn=(N.*X)./(N.^2.*X.^2);

plot(x,fn)

h=legend(’n=0’,’n=1’,’n=2’,’n=3’,’n=4’,’n=5’,’n=6’,...

’n=7’,’n=8’,’n=9’,’n=10’);

Possiamo anche inserire un ciclo for per vedere l’evoluzione dei grafici della successionedata;

for cnt=1:length(n)

plot(x,fn(cnt,:))

hold on

pause(.5)

end

hold off

Possiamo ripetere le operazioni precedenti nel caso delle due funzioni

fn(x) =sin(nt)

n, gn(x) =

√x2 +

1

n.

Prendiamo in considerazione la seconda funzione, lasciando come esercizio il caso dellaprima.

57

Page 58: Disp Matlab

58 CAPITOLO 7. SUCCESSIONI E SERIE DI FUNZIONI

x=linspace(-1,1,100);

n=[0:20];

[X,N]=meshgrid(x,n);

gn=sqrt(X.^2+1./N);

plot(x,abs(x),’r’);

legend(’|x|’)

hold on

for cnt=1:length(n)

plot(x,gn(cnt,:))

pause(.5)

end

hold off

7.2 Serie di funzioni

Consideriamo qui l’esempio della serie di funzioni associata alla successione di funzioni

un(x) =1

1 + x2n;

le funzioni date sono pari, quindi considereremo solo il caso x ≥ 0. Possiamo considerare iseguenti comandi;

x=linspace(0,4,100);

n=[0:5];

[X,N]=meshgrid(x,n);

un=1./(1+x.^(X*N));

serie=cumsum(un);

subplot(2,1,1);

plot(x,un,’LineWidth’,2);

set(gca,’FontSize’,14);

xlabel(’Asse x’,’FontSize’,14);

ylabel(’Asse y’,’FontSize’,14);

h=legend(’n=0’,’n=1’,’n=2’,’n=3’,’n=4’,’n=5’)

grid on

subplot(2,1,2)

plot(X,serie,’LineWidth’,2)

set(gca,’FontSize’,14)

xlabel(’Asse x’,’FontSize’,14);

ylabel(’Asse y’,’FontSize’,14);

h=legend(’n=0’,’n=1’,’n=2’,’n=3’,’n=4’,’n=5’)

Esercizio 7.1 1. Si ripeta l’esercizio con la serie associata alla successione

un(x) = nxxn.

2. Si consideri la serie associata alla successione

un(x) =xn

(1 + x2)n,

discutendone, riportandola ad una serie geometrica, la convergenza.

Page 59: Disp Matlab

7.3. SERIE DI POTENZE 59

7.3 Serie di potenze

Formalmente una serie di potenze e un polinomio di grado eventualmente infinito e viceversa,i polinomi sono serie di potenze in cui i coefficienti sono definitivamente nulli. Vediamo cometrattare il caso della serie di potenze

∞∑n=0

xn

n!;

per ogni k ∈ N, la somma parziale di grado k di tale serie, data da

sk(x) =

k∑n=0

xn

n!;

e un polinomio di grado k. Cosı ad esempio se k = 3 otteniamo

s3(x) = 1 + x+x2

2+x3

6.

In Matlab, il modo di determinare un polinomio e assegnare i coefficienti del polinomio, conla convenzione che il vettore dei coefficienti rappresenta le potenze di x in ordine decrescentedi grado. Ad esempio, il polinomio s3(x) e determinato dal vettore [1/6 1/2 1 1]. Coni seguenti comandi determiniamo il polinomio s10(x) e ne disegna il grafico sull’intervallo[−2, 3].

n=[1 1:10];

a=1./cumprod(n); % coefficienti dei primi 11 termini della serie

b=a(length(a):-1:1); % b e’ il vettore a letto da destra a sinistra

% e rappresenta il polinomio somma parziale

% di grado 10 della serie

x=linspace(-2,3);

plot(x,polyval(b,x))

Il comando polyval(b,x) serve in Matlab per costruire il polinomio, nella variabile x cheha il vettore b come coefficienti. Si provi come esercizio a sovrapporre al grafico precedenteil grafico, tracciato in rosso, della funzione ex; si modifichi anche l’intervallo [−2, 3] con altriintervalli.

Serie con solo termini dispari. Come ulteriore esempio possiamo considerare la seriedi potenze

∞∑n=0

(−1)nx2n+1

(2n+ 1)!;

determineremo il polinomio di grado 101 ne tracceremo il grafico sull’intervallo [0, 7]. Lasomma parziale di grado 101 si ottiene considerando la somma parziale con n = 50. Sipresenta qualche difficolta; la prima consiste nel fatto che il fattoriale (2n + 1)! non esemplicemente cumprod(2*n+1), in quanto questo prodotto mancherebbe i termini pari, chevanno quindi aggiunti (inoltre tale vettore va modificato in modo che i primi due terminirappresentino 0! = 1 e 1! = 1). La seconda difficolta consiste nel fatto che nelle sommeparziali relative alla serie data, le potenze pari sono tutte nulle e tali zeri vanno inseriti nelvettore dei coefficienti. Possiamo procedere come segue, semplicemente costruendo il vettoreb che sia lungo il doppio del vettore a e che inizialmente contenga tutti 0, e i coeffcienti dia vengono inseriti nelle posizioni pari di b,

Page 60: Disp Matlab

60 CAPITOLO 7. SUCCESSIONI E SERIE DI FUNZIONI

n=[0:50];

a=(-1).^n./(cumprod(2*n+1).*cumprod([1 2*n(2:51)]));

b=zeros(1,2*length(a));

b(2:2:length(b))=a;

c=b(length(b):-1:1); % c e’ il vettore che definisce il polinomio

% somma parziale di grado 101

x=0:.01:7;

plot(x,polyval(c,x))

Esercizio 7.2 1. Tracciare il grafico della somma parziale di grado 10 della serie dipotenze

∞∑n=0

xn

sull’intervallo [−.9, .9]. Tracciare nella stessa figura il grafico della funzione y = 11−x ,

inserendo anche una legenda.

2. Disegnare la somma parziale di grado 50 della serie di potenze

∞∑n=0

(−1)nx2n+1

2n+ 1

sull’intervallo [−.5, 1].

3. Disegnare la somma parziale di grado 6 della serie di potenze

∞∑n=0

(−1)nx2n

(2n)!;

si ripeta quanto fatto per la somma parziale di grado 100, tracciando anche il graficodella funzione cosx, considerando vari intervalli di definizione.

7.4 Serie di Taylor

Disegnamo ora nella figura il grafico di sinx e i suoi polinomi di Taylor di grado 1, 3, 11 e51 nell’intervallo [−2π, 2π]; possiamo considerare i seguenti comandi;

n=[0:25]; % length(n)=26

a=(-1).^n./(cumprod(2*n+1).*cumprod([1 2*n(2:length(n))])); %length(a)=26

b=zeros(1,2*length(a)); % length(b)=52

b(2:2:length(b))=a;

c51=b(length(b):-1:1); % length(c51)=52

c11=c51(41:52); % length(c11)=12

c3=c11(9:12); % length(c3)=4

c1=c3(3:4); % length(c1)=2

x=-2*pi:.01:2*pi;

y51=polyval(c51,x);

y11=polyval(c11,x);

y3=polyval(c3,x);

Page 61: Disp Matlab

7.4. SERIE DI TAYLOR 61

y1=polyval(c1,x);

plot(x,sin(x),x,y1,x,y3,x,y11,x,y51);

legend(’sin x’,’Grado 1’,’Grado 3’,’Grado 11’,’Grado 51’);

axis([-2*pi 2*pi -1.5 1.5])

Si noti che il grafico di sinx, che dovrebbe comparire in blu nel grafico, in realta e coper-to completamente dal grafico della funzione y51, mostrando che il polinomio di grado 51rappresenta quasi completamente la funzione sinx nell’intervallo [−2π, 2π].

7.4.1 Polinomi di Taylor simbolici

Possiamo anche utilizzare il symbolic toolbox per trovare l’espansione di Taylor di una datafunzione attorno ad un arbitrario punto x0. A tal proposito e comodo il comando Matlabtaylor; possiamo cosı trovare il polinomio di Taylor di grado 8 di sinx attorno a x0 = 0 equello di grado 10 attorno a x0 = π/4;

syms x

T7=taylor(sin(x),8)

T10=taylor(sin(x),pi/4,10)

Possiamo utilizzare il comando sym2poly per convertire una espressione simbolica di unpolinomio nel vettore dei suoi coefficienti;

t7=sym2poly(T7)

t10=sym2poly(T10)

e viceversa possiamo anche utilizzare il comando poly2sym per il viceversa.

7.4.2 Moltiplicazione di serie e serie di Taylor di prodotti di fun-zione

Se abbiamo due serie di potenze

f(x) =

∞∑n=0

anxn, g(x) =

∞∑n=0

bnxn,

possiamo considerare il loro prodotto h(x) = f(x)g(x) che e ancora una serie di potenze

h(x) =

∞∑n=0

cnxn, cn =

∑h = 0nahbn−h.

Questo prodotto, che e il caso discreto della convoluzione tra due funzioni, viene implemen-tato in Matlab dal comando conv;

f=[1 6 1];

g=[2 -3 5];

h=conv(f,g)

Questo esempio ha come risultato i coefficienti della seguente espressione:

h(x) = f(x)g(x) = (x2 + 6x+ 1)(2x2 − 3x+ 5) = 2x4 + 9x3 − 11x2 + 27x+ 5.

Page 62: Disp Matlab

62 CAPITOLO 7. SUCCESSIONI E SERIE DI FUNZIONI

Possiamo applicare quanto appena detto per il calcolo della serie di Taylor per unafunzione prodotto; ad esempio, possiamo considerare la funzione

h(x) = (x2 − 2x+ 1) sinx

e ne calcoliamo il polinomio di Taylor attorno ad x0 = 0 di grado 51, tracciando i grafici in[−2, 2].

n=[0:25];

a=(-1).^n./(cumprod(2*n+1).*cumprod([1 2*n(2:26)]));

b=zeros(1,2*length(a));

b(2:2:length(b))=a;

c51=b(length(b):-1:1); % Polinomio di Taylor di sin(x)

p=[1 -2 1]; % il polinomio x^2-2x+1

c=conv(p,c51); % il polinomio del prodotto, pero’ troppo lungo

q=c(length(c)-51:length(c)); % fino al grado 51

plot(-2:.01:2,polyval(q,-2:.01:2))

Si provi infine il comando taylortool che racchiude in una interfaccia grafica quantoappena visto.

Esercizio 7.3 1. Nella stessa figura, tracciare il grafico di ex e i suoi polinomi di Taylordi grado 1, 3, 8 e 12 attorno ad x0 = 0 nell’intervallo [−2, 2] (si metta una legenda, siriscalino gli assi se necessario).

2. Nella stessa figura, si tracci il grafico di cosx e i suoi polinomi di Taylor centrati inx0 = 0 di grado 2, 6, 12 e 50 nell’intervallo [−2π, 2π] (si metta una legenda, si riscalinogli assi se necessario, si inserisca una griglia).

3. Nella stessa figura, si tracci il grafico di log(x) e i suoi polinomi di Taylor centrati inx0 = 2 di grado 1, 3, 10 e 20 nell’intervallo [.1, 4] (si metta una legenda, si riscalinogli assi se necessario).

4. Utilizzare il symbolic toolbox per determinare il polinomio di grado 9 attorno a x0 = 0della funzione arctanx; si utilizzi il comando sym2poly per convertire quanto trovatonella rappresentazione vettoriale del polinomio. Tracciare infine nella stessa figura ilgrafico di arctanx e del polinomio trovato nell’intervallo [−1.3, 1.3].

5. Trovare il polinomio di grado 15 attorno a x0 = 0 di (x2 − 1) sinx cosx e si tracci ilgrafico di quanto ottenuto nell’intervallo [−2, 2].

7.5 Serie di Fourier

Ricordiamo che data una funzione f : [0, T ]→ R (estesa poi T–periodica a tutto R, possiamocostruire la serie di Fourier associata definendo i coefficienti

ak =2

T

∫ T

0

f(x) cos(kωx)dx, k ≥ 0,

bk =2

T

∫ T

0

f(x) sin(kωx)dx, k ≥ 1,

Page 63: Disp Matlab

7.5. SERIE DI FOURIER 63

con ω = 2πT frequenza o pulsazione della funzione. La serie sara data da

Sf (x) =a0

2+

∞∑k=1

(ak cos(kωx) + bk sin(kωx)) .

Come esercizio, calcoliamo il polinomio di Fourier di grado 7 associato alla funzione

f(x) = xχ[0,π)(x) + x2χ[−π,0)(x),

definita in [−π, π) e prolungata 2π–periodica. Possiamo ad esempio considerare i seguenticomandi;

x=linspace(-pi,pi,100);

fx=x.^2.*(x<0)+x.*(x>=0);

k=[1:7];

[X,K]=meshgrid(x,k);

[FX,K]=meshgrid(fx,k);

a0=1/pi*trapz(x,fx);

ak=FX.*cos(K.*X);

ak=trapz(x,ak’);

ak=ak./pi;

bk=FX.*sin(K.*X);

bk=trapz(x,bk’);

bk=bk./pi;

[FX,AK]=meshgrid(fx,ak);

[FX,BK]=meshgrid(fx,bk);

SF=a0/2*ones(size(x));

SF(2:(length(k)+1),:)=AK.*cos(K.*X)+BK.*sin(K.*X);

SF=sum(SF);

subplot(2,1,1)

p=plot(x,fx,’b’,x,SF,’r’,...

x-2*pi,fx,’b’,x-2*pi,SF,’r’,...

x+2*pi,fx,’b’,x+2*pi,SF,’r’);

set(p,’LineWidth’,2);

grid on;

set(gca,’FontSize’,14)

title(’f(x) e Serie di Fourier’,’FontSize’,14);

xlabel(’Asse x’,’FontSize’,14)

ylabel(’Asse y’,’FontSize’,14)

legend(’f(x)’,’SF’)

subplot(2,1,2)

bar([ak;bk]’)

grid on

set(gca,’FontSize’,14)

title(’Piano delle frequenze’,’FontSize’,14)

xlabel(’k’,’FontSize’,14)

ylabel(’ak,bk’,’FontSize’,14);

legend(’ak’,’bk’)

Puo essere utile creare un function m–file, che chiameremo CoeffFourier, per il calcolodei coefficienti di Fourier di una data funzione.

Page 64: Disp Matlab

64 CAPITOLO 7. SUCCESSIONI E SERIE DI FUNZIONI

function [a0,ak,bk]=CoeffFourier(x,fx)

%

global T omega k

%

%

[X,K] = meshgrid(x,k);

[FX,K] = meshgrid(fx,k);

%

a0 = 2/T * trapz(x,fx);

%

ak = FX .* cos(omega*X.*K);

ak = trapz(x,ak’);

ak = 2/T * ak;

%

bk = FX .* sin(omega * X.*K);

bk = trapz(x,bk’);

bk = 2/T * bk;

Si noti la definizione delle variabili globali global k T omega, che quindi dovranno esseredefinite prima di richiamare la funzione CoeffsFourier.

Esercizio 7.4 1. Modificare i comandi precedenti per comprendere, nel grafico dellefrequenze, anche il termine a0/2; considerare anche la variante in cui si riportano nelpiano delle frequenze i valori di |ak| e |bk|.

2. Ripetere quanto fatto per la funzione f : (−π, π]→ R, f(x) = x

7.5.1 Fenomeno di Gibbs

Il fenomeno di Gibbs e quel fenomeno che si verifica in prossimita delle discontinuita di unafunzione; quello che si nota, anche visivamente osservandone i grafici, e che in prossimita diun salto della funzione (in prossimita di una discontinuita), si nota che la serie di Fourierproduce un errore che si sposta sempre piu vicino al punto di discontinuita, ma che inintensita non diminuisce se si aumenta il grado di approssimazione. Tale errore resta semprepari ad un decimo il salto della funzione. Verifichiamo questo fenomeno con un esempio.

Consideriamo la funzione f : [−π, π]→ R definita da

f(x) = χ[−π/2,π/2](x)

ed introduciamo l’errore relativo

En(x) =1

|f(x)|

∣∣∣∣∣f(x)− a0

2−

n∑k=1

(ak cos kx+ bk sin kx)

∣∣∣∣∣ .x=linspace(-pi,pi,100);

fx=1.*(x>=-pi/2 & x<=pi/2);

k=[1:30];

[X,K]=meshgrid(x,k);

[FX,K]=meshgrid(fx,k);

Page 65: Disp Matlab

7.5. SERIE DI FOURIER 65

a0=1/pi*trapz(x,fx);

ak=FX.*cos(K.*X);

ak=trapz(x,ak’);

ak=ak./pi;

bk=FX.*sin(K.*X);

bk=trapz(x,bk’);

bk=bk./pi;

[FX,AK]=meshgrid(fx,ak);

[FX,BK]=meshgrid(fx,bk);

SF=a0/2*ones(size(x));

SF(2:(length(k)+1),:)=AK.*cos(K.*X)+BK.*sin(K.*X);

SF=sum(SF);

En=abs((SF-fx))./abs(fx);

subplot(2,1,1);

p=plot(x,fx,’b’,x,SF,’r’,...

x-2*pi,fx,’b’,x-2*pi,SF,’r’,...

x+2*pi,fx,’b’,x+2*pi,SF,’r’);

set(p,’LineWidth’,2);

grid on;

set(gca,’FontSize’,14)

title(’f(x) e Serie di Fourier’,’FontSize’,14);

xlabel(’Asse x’,’FontSize’,14)

ylabel(’Asse y’,’FontSize’,14)

legend(’f(x)’,’SF’)

axis([0,1.6,.5,1.2]) % ci concentriamo sul salto

subplot(2,1,2);

plot(x,En,’g’,’Linewidth’,2)

axis([0,1.6,0,.1])

grid on

set(gca,’FontSize’,14)

title(’Errore relativo’,’FontSize’,14)

xlabel(’Asse x’,’FontSize’,14)

ylabel(’Asse y’,’FontSize’,14)

Page 66: Disp Matlab

66 CAPITOLO 7. SUCCESSIONI E SERIE DI FUNZIONI

Page 67: Disp Matlab

Capitolo 8

Equazioni differenziali

In questo capitolo considereremo alcune equazioni differenziali ordinare e vedremo alcunimetodi per lo studio delle soluzioni, anche in casi in cui il calcolo esplicito della soluzionenon si puo fare.

8.1 Linee di flusso

Possiamo utilizzare il comando quiver per disegnare le linee di flusso di una equazionedifferenziale. L’idea sostanziale risiede nella considerazione che cercare una funzione u :I → R soluzione dell’equazione

u′(t) = f(t, u(t))

equivale a carcare una funzione per la quale nel punto (t, u(t)) la direzione tangente algrafico sia determinata da f , cioe nel punto (t, u(t)) sia assengata la direzione tangente(1, u′(t)) = (1, f(t, u(t))).

Possiamo ad esempio considerare l’equazione

u′(t) = cos(2t+ u(t))

con t ∈ [−5, 5] e u ∈ [.5, 4.5];

[t,u]=meshgrid(-5:.5:5,.5:.5:4.5);

dt=ones(size(t));

du=cos(2*t+u);

quiver(t,u,dt,du);

title(’u’’=cos(2t+u)’);

Per migliorare la resa grafica di quanto ottenuto ad esempio riducendo della meta la lunghez-za dei vettori ottenuti, considerando ad esempio il comando

quiver(t,u,dt,du,.5);

Possiamo anche considerare il Problema di Cauchy associato all’equazione differenziale data;{u′(t) = cos(2t+ u(t))u(t0) = u0

67

Page 68: Disp Matlab

68 CAPITOLO 8. EQUAZIONI DIFFERENZIALI

Un comando Matlab utile in questi casi e il comando streamline, che possiamo utilizzarecome segue. Supponiamo di partire all’istante t0 = −2 e di cercare le quattro soluzioniassociate ai valori iniziali 1, 2, 3, 4; possiamo quindi aggiungere ai comandi precedenti quantosegue;

hold on

u0=[1,2,3,4];

t0=-2*ones(size(u0));

S=streamline(t,u,dt,du,t0,u0);

set(S,’Color’,’r’)

Le linee in rosso saranno quindi i grafici delle soluzioni dei Problemi di Cauchy associati.

Esercizio 8.1 1. Utilizzare il comando quiver per tracciare le linee di flusso del campo

u′ = 2t2 − u3;

si utilizza la finestra [−1, 1]2 con incrementi .2 in ogni direzione.

2. Utlizzare il comandott quiver per tracciare le linee di flusso di

u′ = sin(t+ u)

su di un rettangolo contenente l’origine (si facciano varie prove con differenti valoriper le dimensioni della finestra e degli incrementi).

8.2 Il metodo numerico ode45

In Matlab esistono vari pacchetti per risolvere le equazioni differenziali, tra questi segna-liamo qui il pacchetto ode45. Come esempio consideriamo il Problema di Cauchy{

y′(t) = t+2 cos ty2(t)

y(0) = 2

la cui soluzione analitica e data da u(t) = 3√

32 t

2 + 6t+ 8. Per la soluzione numerica, che

sara una soluzione approssimata, possiamo procedere come segue;

[t,u]=ode45(fcnchk(’(t+2*cos(t))./u.^2’),[0 5],2);

U=(3*t.^2/2+6*sin(t)+8).^(1/3);

plot(t,u,’b’,t,U,’r’)

[u(length(t)),U(length(t))]

I precedenti comandi danno le soluzioni numerica ed analitica in [0, 5] e ne disegna il grafico;i due grafici sono sovrapposti, mostrando quanto la soluzione numerica sia accurata. L’ul-timo comando confronta infatti le due soluzioni in t = 5. Nell’eventualita che la soluzionenumerica non sia cosı accurata, si puo anche modificare la griglia di calcolo, che nei coman-di precedenti viene scelta di default dal comando ode45; possiamo ad esempio sostituire laprima riga di comando con la seguente;

[t,u]=ode45(fcnchk(’(t+2*cos(t))./u.^2’),[0:.2:5],2);

Page 69: Disp Matlab

8.3. INTEGRAZIONE SIMBOLICA DI EQUAZIONI DIFFERENZIALI 69

Nel comando ode45 non e necessario utilizzare tt fcnchk, ma possiamo anche utilizzare unfunction m–file, con la sintassi

[t,u]=ode45(’f’,[a b],u0);

dove la funzione f viene definita in un m–file che avra la forma

function uprimo=f(x,y)

uprimo=(t+2*cos(t))./u.^2;

8.3 Integrazione simbolica di equazioni differenziali

In Matlab esiste anche il symbolic toolbox dsolve per trovare la soluzione generale di unaequazione differenziale, che viene utilizzato con la seguente sintassi

syms t y

g=dsolve(’Dy=f(t,y)’,’t’)

per l’integrale generale, mentre useremo la sintassi

syms t y

g=dsolve(’Dy=f(t,y),y(a)=b’,’t’)

per la soluzione del Problema di Cauchy.

Ad esempio, se vogliamo risolvere l’equazione y′ = t2 + y, possiamo considerare;

syms t y

g=dsolve(’Dy=t^2+y’,’t’)

Il risultato sara una funzione con un parametro C2; se vogliamo determinare la soluzione inun dato punto, ad esempio t = 0, possiamo considerare il comando

a=subs(g,0)

da cui si puo capire il valore della costante C2 per cui si ha che la soluzioni soddisfi unacerca condizione iniziale per t = 0.

Se invece vogliamo considerare direttamente un Problema di Cauchy, ad esempio{u′ = u− tu(−1.5) = −1

possiamo considerare i seguenti comandi;

syms t y

g=dsolve(’Dy=y-t,y(-1.5)=-1’,’t’)

ezplot(g,[-2 2])

Page 70: Disp Matlab

70 CAPITOLO 8. EQUAZIONI DIFFERENZIALI

8.4 Sistemi di equazioni differenziali

Si possono anche considerare sistemi di equazioni differenziali, come ad esempio{x′(t) = 2x(t)− 3y(t), x(0) = 1y′(t) = 2x(t) + y(t), y(0) = −1;

utilizzeremo ode45 per determinare una soluzione approssimata per t ∈ [0, 1]. L’equazionee un’equazione del prim’ordine u′ = F (t, u), con u vettore di componenti (x, y) ed F (t, u) =(2u1 − 3u2, 2u1 + u2); quindi la procedura non sara molto differente da quella vista per leequazioni del prim’ordine. Iniziamo col creare un m–file

function xprimo=F(t,x)

xprimo=[2*x(1)-3*x(2);2*x(1)+x(2)];

e poi richiamare tale m–file coi seguenti comandi;

t0=0; t1=1;

[t,x]=ode45(’F’,[t0 t1], [1,-1]);

figure(1); clf;

plot(x(:,1),x(:,2));

n=length(t);

T=[’(’,num2str(x(n,1)),’, ’,num2str(x(n,2)),’)’];

text(x(n,1),x(n,2),T);

figure(2); clf;

plot(t,x(:,1),t,x(:,2));

legend(’x(t)’,’y(t)’);

8.5 Equazioni di ordine superiore

Possiamo utilizzare quanto visto per i sistemi per risolvere equazioni di ordine superiore alprimo; ad esempio, supponiamo di voler risolvere l’equazione

x′′(t) + 6x(t)x′(t)− x2(t) = sin(x(t)), x(0) = 1, x′(0) = 2

con t ∈ [0, 5]. Possiamo quindi introdurre la variabile y(t) = x′(t) ed ottenere quindi ilsistema {

x′(t) = y(t), x(0) = 1y′(t) = −6x(t)y(t) + x2(t) + sin(x(t)), y(0) = 2.

Definiamo quindi l’m–file

function xprimo=higher(t,x)

xprimo=[x(2); -6*x(1).*x(2)+sin(x(1))];

e quindi costruiamo la soluzione;

[t,x]=ode45(’higher’,[0 5], [1 2]);

plot(t,x(:,1),t,x(:,2));

legend(’x(t)’,’y(t)’)

Page 71: Disp Matlab

8.6. SOLUZIONI DEFINITE IMPLICITAMENTE 71

Esercizio 8.2 1. Si determini la soluzione dell’equazione differenziale;

y′(t) =ty(t) + 2t

t2 + 1

e tra queste si trovi quella che soddisfa y(1) = 2. Utilizzare ode45 per trovare lasoluzione approssimata con t ∈ [1, 6], disegnando in blu la soluzione esatta ed inrosso quella approssimata (si usi anche il comando text per dare i valori esatti edapprossimati della soluzione in t = 6).

2. Si utilizzi ode45 per determinare la soluzione del problema

y[(t) = t− y(t) + 1, y(0) = 3

con t ∈ [0, 3]. Si disegni la soluzione, indicando il valore della soluzione per t = 3.

3. Si utilizzi ode45 per trovare una soluzione approssimata del sistema{x′(t) = 3x(t) + 4y(t), x(0) = 1y′(t) = 3x(t) + 2y(t), y(0) = 2

con t ∈ [0, 1]. Disegnare nello stesso grafico le soluzioni x = x(t) e y = y(t), ed in unsecondo grafico la curva y = x.

4. Si consideri l’equazione differenziale del secondo ordine

x′′(t) + tx′(t) + t2x(t) = 0;

utilizzare ode45 per trovare la soluzione in t ∈ [0, 3] con i seguenti dati iniziali,specificando per la variabile indipendente i seguenti valori t=[0:.01:3];

a) x(0) = 1, x′(0) = 0 e si chiami tale soluzione f ;

b) x(0) = 0, x′(0) = 1 e si chiami tale soluzione g;

c) x(0) = 2, x′(0) = −3 e si chiami tale soluzione h.

Si noti che, dalla teoria delle equazioni differenziali lineari (a coefficienti non necessari-amente costanti), la soluzione generica del problema e data da una combinazione di fe g; verificare che effettivamente la soluzione h e data da h = 2f − 3g, dove 2 = h(0)e −3 = h′(0). Si disegni quindi il grafico di 2f − 3g e di h nella stessa figura.

8.6 Soluzioni definite implicitamente

Puo succedere a volte che la soluzione di un’equazione differenziale la si trova ma non definitain modo esplicito, ma implicitamente nella forma

F (t, y(t)) = c,

per una qualche funzione F e c ∈ R. In tal caso puo essere utile utilizzare il comandocontour per lo studio di tali soluzioni.

Supponiamo ad esempio di voler risolvere l’equazione

y′(t) =t√t2 + 1

y(t)ey(t);

Page 72: Disp Matlab

72 CAPITOLO 8. EQUAZIONI DIFFERENZIALI

la soluzione di tale equazione si trova in modo diretto (l’equazione e a variabili separabili)ed e determinata da

(8.1) (y(t)− 1)ey(t) − (t2 + 1)32

3= c,

con c che dipende dalle condizioni iniziali. Disegneremo nel quadrato [−3, 3]2 le curve chesono i grafici di tale equazione con c che assume i valori da −4 a 2 saltando di 1/2;

[t,y]=meshgrid(-3:.1:3);

z=(y-1).*exp(y)-1/3.*(t.^2+1).^(3/2);

c=[-4:.5:2]; % i livelli

contour(t,y,z,c)

Possiamo aggiungere anche le etichette al grafico precedente, sostituendo l’ultima riga con

[C,h]=contour(t,y,z,c);

clabel(C,h);

Se abbiamo un Problema di Cauchy, cioe se imponiamo una condizione iniziale y(t0) = y0,allora il livello c sara univocamente determinato; ad esempio, se nell’esempio precedenteimponiamo y(0) = 1/10, allora da (8.1) ricaviamo che c = − 9

10e1/10 − 1

3 ;

[t,y]=meshgrid(-3:.1:3);

z=(y-1).*exp(y)-1/3.*(t.^2+1).^(3/2);

c=-.9*exp(.1)-1/3;

contour(t,y,z,[c c]);

Si noti che vengono disegnate due curve; si provi a capire come mai e soprattutto, qualedelle due curve e quella che definisce il grafico della funzione.

Esercizio 8.3 Trovare la soluzione generale dell’equazione

t2y′(t) =t2 + 1

3y2(t) + 1

ed utilizzare il comando contour per tracciare il grafico di dieci curve che definiscono lesoluzioni; in un altra figura, si tracci il grafico della curva che definisce la soluzione delproblema y(1) = 2.

8.7 Vibrazioni meccaniche

In questa sezione considereremo il problema delle vibrazioni meccaniche; con questo in-tendiamo il problema della determinazione e dello studio delle soluzioni dell’equazionedifferenziale del secondo ordine a coefficienti costanti

(8.2) y′′(t) + 2δy′(t) + ω20y(t) = B cos(ωt),

dove δ ∈ R rappresenta un fattore di smorzamento, ω0 la frequenza propria del sistema, ωe B rispettivamente la frequenza e l’intensita del termine forzante. Nel caso B = 0 si parladi vibrazioni libere, mentre per B 6= 0 si parla di vibrazioni forzate; il caso δ = 0 siamo nelcaso di vibrazioni non smorzate, con δ 6= 0 di vibrazioni smorzate. Tratteremo prima il casodi vibrazioni non smorzate, quindi il caso delle vibrazioni smorzate.

Page 73: Disp Matlab

8.7. VIBRAZIONI MECCANICHE 73

8.7.1 Vibrazioni non smorzate

Vibrazioni libere Possiamo considerare ad esempio il problema di Cauchy associato a(8.2) con B = 0 e dati iniziali y(0) = y′(0) = 1, la cui soluzione e data da

y(t) = cos(ω0t)− sin(ω0t);

per disegnare il grafico della soluzione possiamo considerare i seguenti comandi:

t=linspace(0,40,100);

y0=1;

v0=1;

omega0=1;

A=sqrt(y0^2+v0^2/omega0^2);

c=y0/A;

s=-v0/(A*omega0);

yt=A*(c*cos(omega0*t)-s*sin(omega0*t));

plot(t,yt)

grid on

Vibrazione forzata; caso ω 6= ω0 e battimenti Come variazione del problema prece-dente consideriamo il caso in (8.2) con B = 1 e ω = 0; i dati iniziali che consideriamo sonosempre y(0) = y′(0) = 1, e quindi la soluzione sara data, se ω0 = 1, da

y(t) = 1 + sin(t).

Tracciamo il grafico di tale soluzione nel caso di ω0 = 1;

t=linspace(0,40,100);

y0=1;

v0=1;

omega0=1;

B=1;

A=sqrt((y0-B/omega0^2)^2+v0^2/omega0^2);

c=(y0-B/omega0^2)/A;

s=-v0/(A*omega0);

yt=A*(c*cos(omega0*t)-s*sin(omega0*t))+B/omega0^2;

ht=B*ones(length(t));

plot(t,yt,t,ht,’r’);

grid on

Qui si e disegnato in rosso il grafico del termine forzante, mentre in blu si ha la soluzione; sinota quindi semplicemente un’oscillazione libera, traslata di uno rispetto alla posizione diequilibrio.

Possiamo modificare il precedente caso considerando ω = 1/√

2; in tal caso la soluzionesara data da

y(t) = 2 cos

(t√2

)− cos t+ sin t;

siamo in presenza del fenomeno noto come battimenti, e possiamo visualizzare tale fenomenoconsiderando i seguenti comandi;

Page 74: Disp Matlab

74 CAPITOLO 8. EQUAZIONI DIFFERENZIALI

t=linspace(0,200,1000);

y0=1;

v0=1;

omega0=1;

B=1;

A=sqrt((y0-(B/omega0)^2)^2+(v0^2/omega0)^2);

c=(y0-(B/omega0)^2)/A;

s=-v0/(A*omega0);

omega=1/sqrt(2);

yt=A*(c*cos(omega0*t)-s*sin(omega0*t))+B*cos(omega*t)/(omega0^2-omega^2);

ht=B*cos(omega*t);

plot(t,yt,t,ht,’r’);

grid on

Ancora, in rosso abbiamo disegnato il grafico del termine forzante, mentre in blu abbiamo ilgrafico della soluzione; si nota quindi nei battimenti un’oscillazione con ampiezza variabileme in qualche modo con ampiezza periodica. Si provi come esercizio a considerare differentivalori di ω 6= ω0, considerando in particolare valori di ω sempre piu vicini, ma non uguali,a quelli di ω0.

Vibrazione forzata; caso ω = ω0 e risonanza Se il termine forzante e ha una pulsazioneω = ω0, allora si verificare il fenomeno della risonanza; in questo caso, la soluzione particolareva moltiplicata per t, e quindi la soluzione di (8.2) con y(0) = y′(0) = 1 e data da

y(t) = sin(t) +t

2sin(t).

t=linspace(0,250,1000);

y0=1;

v0=1;

omega0=1;

B=1;

A=sqrt((y0-B/omega0^2)^2+v0^2/omega0^2);

c=(y0-B/omega0^2)/A;

s=-v0/(A*omega0);

yt=A*(c*cos(omega0*t)-s*sin(omega0*t))+B*t.*sin(omega0.*t)/(2*omega0);

hx=B*cos(omega0*t);

plot(t,yt,t,ht,’r’)

grid on

8.7.2 Vibrazioni smorzate

Trattiamo qui il caso in cui δ 6= 0; il polinomio caratteristico associato sara dato da p(λ) =λ2 + 2δλ + ω2

0 e il comportamento delle soluzioni dipendera dal valore del discriminante∆/4 = δ2−ω2

0 . Si parla di smorzamento supercritico se |δ| > ω0 (in questo caso il polinomiocaratteristico avra due radici reali distinte), smorzamento critico se δ = ω0 (una sola radicereale del polinomio caratteristico) ed infine smorzamento subcritico se |δ| < ω0 (avremo quidue radici complesse coniugate).

Page 75: Disp Matlab

8.7. VIBRAZIONI MECCANICHE 75

Vibrazione smorzata: caso supercritico Consideriamo l’equazione (8.2) con δ = 1.2,ω0 = 1, B = 0 e il problema di Cauchy con dati iniziali y(0) = 1, y′(0) = −3; in tal caso lasoluzione e data da

y(t) =9−√

11

2√

11e−δt

(e

2√

1110 t − e−

2√

1110 t

)e il suo grafico lo possiamo disegnare come segue;

t=linspace(0,10,100);

y0=1;

v0=-3;

omega0=1;

delta=1.2;

r1=-delta-sqrt(delta^2-omega0^2);

r2=-delta+sqrt(delta^2-omega0^2);

c1=(v0-y0*r2)/(r1-r2);

c2=y0-c1;

yt=c1*exp(r1*t)+c2*exp(r2*t);

plot(t,yt);

grid on

Vibrazione smorzata: caso critico Consideriamo come esempio il caso ω0 = δ = 1,B = 0 e il problema di Cauchy con dati iniziali y(0) = 1, y′(0) = −3; la soluzione e data da

y(t) = e−t(1− 2t),

il cui grafico sara dato da

t=linspace(0,10,100);

y0=1;

v0=-3;

omega0=1;

yt=y0*exp(-omega0*t)+(v0+omega0*y0)*t.*exp(-omega0*t);

plot(t,yt);

grid on

Vibrazione smorzata: caso subcritico Consideriamo il caso ω0 = 1, δ = 0.1, B = 0con dati iniziali y(0) = y′(0) = 1; la soluzione e

y(t) = e−δt

(cos($t)−

√11

3sin($t)

),

dove $ =√ω2

0 − δ2. Il grafico di tale soluzione sara dato da

t=linspace(0,80,1000);

y0=1;

v0=1;

omega0=1;

delta=0.1;

omega=sqrt(omega0^2-delta^2);

Page 76: Disp Matlab

76 CAPITOLO 8. EQUAZIONI DIFFERENZIALI

A=sqrt(y0^2+((v0+delta*y0)/omega)^2);

c=y0/A;

s=-(omega0+delta*y0)/(A*omega);

yt=A*exp(-delta*t).*(c*cos(omega*t)+s*sin(omega*t));

plot(t,yt);

grid on

Vibrazione smorzata forzata Chiudiamo con un esempio di vibrazione smorzata contermine forzante; consideriamo solo il caso subcritico con ω0 = 1, δ = 0.1, termine forzantecon B = 1/2, ω = 1/3 e risolviamo il problema di Cauchy con dati iniziali y(0) = y′(0) = 1.La soluzione e data da

y(t) = Ae−δt (c cos($t) + s sin($t)) +2δωB

a2sin(ωt) +

B(ω20 − ω2)

a2cos(ωt)

dove

a =√

(ω20 − ω2)2 + 4δ2ω2, $ =

√ω2

0 − δ2, A =

√1 +

(1 + δ

$

)2

,

c =1

A, s = −

1A + cδ

$;

tracciamo il grafico di y(t) assieme al grafico del termine forzante (in rosso);

t=linspace(0,80,1000);

y0=1;

v0=1;

omega0=1;

delta=0.1;

omega=1/3;

B=1/2;

a=sqrt((omega0^2-omega^2)^2+4*delta^2*omega^2);

nu=sqrt(omega0^2-delta^2);

A=sqrt(y0^2+((v0+delta*y0)/nu)^2);

c=y0/A;

s=-(v0/A+c*delta)/nu;

yt=A*exp(-delta*t).*(c*cos(nu*t)+...

s*sin(nu*t))+(2*delta*omega*B/a^2)*sin(omega*t)...

+(B*(omega0^2-omega^2)/a^2)*cos(omega*t);

plot(t,yt,t,B*cos(omega*t),’r’);

grid on

Si noti come la soluzione ha un certo andamento per un certo intervallo di tempo, mentrepoi a regime si comporta come il termine forzante.

Page 77: Disp Matlab

Appendice A

Derivazione numerica

Si puo utilizzare Matlab per il calcolo approssimato delle derivate di una funzione data;sostanzialmente, una funzione altro non sara che un vettore di valori. Possiamo quindi ri-cavare il valore approssimato delle derivate semplicemente utilizzando i rapporti incrementalidella funzione.

A.1 Derivata prima

Per il calcolo della derivata prima della funzione f possiamo utilizzare la formula;

∆f

∆x(x) =

f(x+ ∆x)− f(x)

∆x.

Il valore ottenuto sara necessariamente un valore approssimato di f ′(x) e tale errore saratanto piu piccolo quanto piu piccolo sara il valore di ∆x. Si puo quindi calcolare la derivatadi sinx in [0, 2π] mediante questo metodo:

x=linspace(0,2*pi,1000);

y=sin(x);

deltax=diff(x); deltay=diff(y);

derivata=deltay./deltax;

Ovviamente, dato che ∆x e costante lungo la griglia, avremmo anche potuto considerare ilcomando;

derivata=deltay/(2*pi*1000):

Possiamo comparare il risultato ottenuto tracciando il grafico del vettore derivata assiemeal grafico della funzione cosx;

plot(x(1:999),derivata,’-’,x,cos(x),’:’);

legend(’Derivazione Numerica’,’Derivata’);

title(’Derivazione numerica ed analitica di y=sin(x)’);

xlabel(’x’); ylabel(’y’);

Il comando diff non fa altro che ottenere da un vettore dato un vettore che contienele differenze tra due entrate consecutive del vettore, e viene quindi prodotto un vettore la

77

Page 78: Disp Matlab

78 APPENDICE A. DERIVAZIONE NUMERICA

cui lunghezza e di uno inferiore rispetto alla lunghezza del vettore di partenza. Per questomotivo abbiamo dovuto considerare, nel grafico della derivata, il vettore x(1:999).

Come varianti del comando diff si puo considerare il comando diff(x,2), che, dato adesempio il vettore x=[a,b,c,d], produce il vettore [c-2b+a,d-2c+b]; tale comando coincidecol comando diff(diff(x)), e servira in seguito per il calcolo delle derivate seconde. Il comandodiff(x,1) coincide con diff(x).

Allo stesso modo, si possono calcolare numericamente le derivate parziali di una funzionedi piu variabili; cosı ad esempio, per una funzione di due variabili f(x, y), possiamo ottenerele seguenti approssimazioni

∂f

∂x(x, y) ∼ f(x+ ∆x, y)

∆x,

∂f

∂y(x, y) ∼ f(x, y + ∆y)

∆y.

A.2 Derivata seconda

Per il calcolo della derivata seconda si puo procedere essenzialmente come nella sezioneprecedente, e cioe utilizzando il rapporto incrementale della derivata prima approssimata

f ′′(x) ∼f′(x+ ∆x)− f ′(x)

∆x∼

f(x+2∆x)−f(x+∆x)∆x − f(x+∆x)−f(x+)

∆x

∆x

=f(x+ 2∆x)− 2f(x+ ∆x) + f(x)

(∆x)2.

Quindi, se volessimo ad esempio calcolare la derivata seconda numerica di f(x) = ex in [0, 3]potremmo considerare i seguenti comandi;

x=linspace(0,3,1000);

deltax=3/1000;

derivata2=diff(exp(x),2)/(deltax)^2;

Ovviamente, allo stesso modo potremmo ottenere le derivate seconde numeriche per unafunzione di due variabili;

∂2f

∂x2(x, y) ∼ f(x+ 2∆x, y)− 2f(x+ ∆x, y) + f(x, y)

(∆x)2,

e∂2f

∂y2(x, y) ∼ f(x, y + 2∆y)− 2f(x, y + ∆y) + f(x, y)

(∆y)2.

mentre la derivata seconda mista sara data da

∂2f

∂x∂y(x, y) ∼ f(x+ ∆x, y + ∆y)− f(x+ ∆x, y)− f(x, y + ∆y) + f(x, y)

∆x∆y.

Come esercizio, si scrivano comandi Matlab che calcolino le derivate seconde di una fun-zione di due variabili, dove si utilizzi il comando diff(x,2) e diff(y,2) per il calcolo dellederivate seconde non miste.

Page 79: Disp Matlab

Appendice B

Il metodo di integrazione diSimpson

Con tale nome di denota il metodo di approssimazione dell’integrale di una funzione f :[a, b]→ R ottenuto suddividendo l’intervallo in n = 2m parti uguali, prendendo cioe h = b−a

ncome lunghezza dei sottointervalli e sostituendo le somme di Riemann con le somme

Sn =h

3

(f(a) + 4

m∑k=1

f(x2k−1) + 2

m−1∑k=1

f(x2k) + f(b)

)

dove xk = kh, k = 0, . . . , n sono i punti estremi degli intervalli che definiscono la partizione.Le somme di Simpson possono quindi essere calcolate utilizzando in Matlab i seguenticomandi;

h=(b-a)/n;

xdispari=a+h:2*h:b-h;

xpari=a+2*h:2*h:b-2*h;

Sn=h/3*(f(a)+4*sum(f(xdispari))+2*sum(f(xpari))+f(b));

Si potrebbe verificare una tale formula, confrontandola anche con le somme di Riemann econ il valore esatto dell’integrale, sull’integrale∫ 2

0

1

1 + x2dx

con n = 20. Una possibile successione di comandi diventa quindi;

a=0; b=2; n=20;

h=(b-a)/n;

xdispari=a+h:2*h:b-h;

xpari=a+2*h:2*h:b-2*h;

fdispari=1./(1+xdispari.^2);

fpari=1./(1+xpari.^2);

S20=h/3*(1/(1+a^2)+4*sum(fdispari)+2*sum(fpari)+1/(1+b^2));

79

Page 80: Disp Matlab

80 APPENDICE B. IL METODO DI INTEGRAZIONE DI SIMPSON

Page 81: Disp Matlab

Appendice C

Metodo di Eulero per leequazioni differenziali

Consideriamo qui il metodo noto come metodo di Eulero per il calcolo approssimato dellesoluzioni delle equazioni differenziali del prim’ordine; prendiamo in considerazione l’inter-vallo [−5, 5], n = 40 e il problema ai valori iniziali{

y′(t) = cos(2t+ y(t))y(−5) = 4.

Si usa quindi il seguente metodo:

t0=-5; y0=5; tf=5;

n=40; % numero di passi

h=(tf-t0)/n; % ampiezza del passo

t=t0:h:tf;

y=y0*ones(size(t));

for i=1:n

y(i+1)=y(i)+h*cos(2*t(i)+y(i));

end

plot(t,y)

Si noti che il metodo si basa semplicemente sulla scrittura della derivata discreta

y(ti+1)− y(ti)

h∼ y′(ti) = cos(2ti + y(ti)).

Possiamo anche creare un function m–file per l’implementazione del metodo di Eulero; talefunction avra in ingresso la funzionef che compare in y′ = f(t, y), il punto iniziale t0 e ilpunto finale t1, il dato iniziale y0 e il numero di passi n. Possiamo quindi scrivere

function [t,y]=Euler(f,t0,t1,y0,n)

% Metodo di Eulero per il problema

%

% y’=f(t,y)

% y(t0)=y0

81

Page 82: Disp Matlab

82 APPENDICE C. METODO DI EULERO PER LE EQUAZIONI DIFFERENZIALI

%

% nell’intervallo [t0,t1] con n passi

%

% Utilizzo; si scrivere

%

% y=Euler(’f’,t0,t1,y0,n)

%

% dove f deve essere una funzione \textsc{Matlab} built-in

% o il nome, racchiuso tra apici, di un function m--file.

% In uscita si ottiene un vettore t dei valori della variabile

% indipendente e il vettore y dei valori della variabile dipendente.

% Il comando plot(t,y) produce il grafico della soluzione.

if nargin<5

error(’Mancano dei parametri. Usare: y=Euler(f,t0,t1,y0,n)’)

end

f=fcnchk(f,’vectorized’);

h=(t1-t0)/n;

t=t0:h:t1;

y=y0*ones(size(t));

for i=1:n

y(i+1)=y(i)+h*f(t(i),y(i));

end

A questo punto, se vogliamo utilizzare il file appena creato per risolvere il problema diCauchy precedente, possiamo semplicemente considerare

[t,y]=Euler(’cos(2*t+y)’,-5,5,4,40);

plot(t,y)

Come variante, possiamo anche considerare disegnare la soluzione trovata in rosso sulcampo direzionale definito dall’equazione y′ = cos(2t+ y);

[t,y]=meshgrid(-5:.5:5,.5:.5:4.5);

dt=ones(size(t));

dy=cos(2*t+y);

quiver(t,y,dt,dy);

hold on;

[t,y]=Euler(’cos(2*t+y)’,-5,5,4,40);

plot(t,y,’r’); hold off;

title(’y’’=cos(2t+y)’);

Page 83: Disp Matlab

Appendice D

Integrazione di equazionidifferenziali per serie

Si possono utilizzare anche le serie di funzioni per determinare le soluzioni di determinateequazioni differenziale; presentiamo qui due esempi, uno in cui utilizziamo le serie di Taylor,ed uno in cui utilizzeremo le serie di Fourier.

D.1 Serie di Taylor

A volte puo essere utile la teoria delle serie di Taylor per la determinazione di una soluzionedi un’equazione differenziale; infatti, se ad esempio consideriamo l’equazione

(D.1) y′′(t)− 2ty′(t) + y(t) = 0, y(0) = 1, y′(0) = 0,

con t ∈ [0, 1], allora possiamo cercare la soluzione nella forma

(D.2) y(t) =

∞∑n=0

antn.

Le condizioni iniziali sono condizioni sui primi due coefficienti, a0 = 1, a1 = 0; per ladeterminazione degli altri coefficienti, basta inserire (D.2) in (D.1) per trovare l’equazionealgebrica

an+2 =2n− 1

(n+ 1)(n+ 2)an, n ≥ 0.

Si noti che quindi a partire da a0 si costruiranno tutti i cefficienti di indice pari, mentrepartendo da a1 si costruiranno quelli di indice dispari. Possiamo quindi considerare i seguenticomandi, che ci daranno il polinomio di Taylor di grado 10; l’unica accortezza e che nellacostruzione dei coefficienti col comando a=b(11:-1:1), avremo la relazione a(n) = an−1

b=[1 0]; % Condizioni iniziali

% coefficienti di indice pari

for n=0:2:8

b(n+3)=(2*n-1)/((n+1)*(n+2))*b(n+1);

end

83

Page 84: Disp Matlab

84 APPENDICE D. INTEGRAZIONE DI EQUAZIONI DIFFERENZIALI PER SERIE

% coefficienti di indice dispari

for n=1:2:7

b(n+3)=(2*n-1)/((n+1)*(n+2))*b(n+1);

end

a=b(11:-1:1); % invertiamo i coefficienti

plot(-2:.01:2,polyval(a,-2:.01:2))

Si provi con gradi maggiori, come ad esempio 30 o 100 e si notera che la velocita di calcoloe comunque elevata.

Esercizio D.1 1. Disegnare il polinomio di Taylor di grado 10 della soluzione del prob-lema di Cauchy {

y′′(t)− 2ty′(t) + y(t) = 0,y(0) = 0, y′(0) = 1,

con t ∈ [−2, 2].

2. Disegnare il polinomio di Taylor di grado 10 della soluzione del problema di Cauchy{y′′(t)− 2ty′(t) + y(t) = 0,y(0) = 1, y′(0) = 1,

con t ∈ [−2, 2].

D.2 Serie di Fourier

Supponiamo di voler risolvere la seguente equazione differenziale lineare del secondo ordinea coefficienti costanti

y′′(t) + ω20y(t) = b(t)

e supponiamo che b sia una funzione T–periodica sviluppabile in serie di Fourier; allora possi-amo cercare la soluzione particolare dell’equazione utilizzando le serie di Fourier. Cerchiamoquindi una funzione y sviluppata in serie di Fourier mediante

y(t) =α0

2+

∞∑k=1

[αk cos(kωt) + βk sin(kωt)] ,

dove ω = 2πT , in modo tale che

y′′(t) + ω20 y(t) = b(t) =

a0

2+

∞∑k=1

[ak cos(kωt) + bk sin(kωt)] ,

con

ak =2

T

∫ T/2

−T/2b(t) cos(kωt)dt, k ≥ 0,

bk =2

T

∫ T/2

−T/2b(t) sin(kωt)dt, k ≥ 1.

Page 85: Disp Matlab

D.2. SERIE DI FOURIER 85

Siamo ricondotti quindi ad un sistema (infinito) che coinvolge i coefficienti in quanto larichiesta che y sia soluzione particolare si traduce nell’equazione

α0ω20

2+

∞∑k=1

[αk(ω2

0 − k2ω2) cos(kωt) + βk(ω20 − k2ω2) sin(kωt)

]=a0

2+

∞∑k=1

[ak cos(kωt) + bk sin(kωt)] .

Troviamo quindi che

α0 =a0

ω20

,

mentre per k ≥ 1, se ω20 6= k2ω2, allora

αk =ak

ω20 − k2ω2

, βk =bk

ω20 − k2ω2

.

Se infine esiste m ≥ 1 per cui ω20 = m2ω2, allora siamo nella situazione in cui va trovata la

soluzione particolare dell’equazion2

y′′(t) + ω20y(t) = am cos(ω0t) + bm sin(ω0t);

tale soluzione e data da

− am2ω0

t cos(ω0t) +bm2ω0

t sin(ω0t).

In definitiva avremo che la soluzione particolare sara data da

y(t) =a0

2ω20

− am2ω0

t cos(ω0t) +bm2ω0

t sin(ω0t)+

+

∞∑k=1,k 6=m

[ak

ω20 − k2ω2

cos(kωt) +bk

ω20 − k2ω2

sin(kωt)

].

Come esempio possiamo considerare la funzione b : [−π, π) → R studiata nella Sezione7.5 associata all’equazione differenziale

y′′(t) + y(t) = b(t),

dove su [−π, π) la funzione b e data da b(t) = tχ[0,π)(t) + t2χ[−π,0)(t), e poi viene estesa2π–periodica. Abbiamo quindi T = 2π e ω = 1. Abbiamo anche ω0 = 1, e quindi lo sviluppoin serie di Fourier della soluzione sara

y(t) =a0

2− a1

2t cos t+

b12t sin t+

+

∞∑k=2

[ak

1− k2cos(kt) +

bk1− k2

sin(kt)

].

Come esercizio, si scriva in Matlab dei comandi che calcolino i coefficienti di y e si disegniil grafico della soluzione particoalre cosı ottenuta.