La funzione (x,y) → 1/(2π)*exp(-(x^2+y^2)/2) definita su (−∞,∞)×(−∞,∞)
è una funzione densità (vedi 
Sistemi di variabili casuali).
Qui puoi vedere come tracciarne il grafico con R; prova a tracciarlo.
Con WolframAlpha
puoi verificare che si tratta di una funzione di densità battendo:
integrate 1/(2*pi)*exp(-(x^2+y^2)/2) dx dy, x=-inf..inf, y=-inf..inf
Che cosa ottieni? Come puoi verificara la cosa con R?
(A) Prova a calcolare in modo analogo l'integrale sul triangolo di vertici (0,0), (1,0) e (1,1) della funzione (x,y) → x·y.
(B) Prova a calcolare in modo analogo l'integrale sul quadrato [0,π/2]×[0,π/2] della funzione (x,y) → sin(x)+cos(y).
(C) Prova a calcolare in modo analogo l'integrale sulla figura del primo quadrante delimitata dalle curve y=x² e x=y²
della funzione (x,y) → x·y².
Con R (vedi qui) basta battere (prendendo
per le x e le y l'intervallo [-20,20] e ad esempio 500 passi):
source("http://macosa.dima.unige.it/r.R") # If I have not already loaded the library
F = function(x,y) 1/(2*pi)*exp(-(x^2+y^2)/2)
u = INTEGRAL(F,-20,20, -20,20, 500); u
# 1
(A) Procediamo prima con R:
H = function(x,y) ifelse(y<=x & x>=0, x*y, 0)
u1 = INTEGRAL(H, 0,1, 0,1, 200); u1
# 0.1258333
u2 = INTEGRAL(H, 0,1, 0,1, 400); u2; u1-u2
# 0.1254167 0.0004166621
u3 = INTEGRAL(H, 0,1, 0,1, 800); u3; u2-u3
# 0.1252083 0.0002083328
u4 = INTEGRAL(H, 0,1, 0,1, 1600); u4; u3-u4
# 0.1251042 0.0001041666
# Vedo che la differenza man mano si dimezza. Posso prendere:
D=1/2; S=D; for(i in 1:100) {D=D/2; S=S+D}; S
# 1
# anche se sapevo che 1/2 + 1/2^2 + 1/2^3 + ... = 1
u4-(u3-u4)*S
# 0.125
# Per il grafico completo la definiz. della funzione fuori dal dominio di integrazione
BF=3; HF=3 # definisco (circa in pollici) base e altezza della finestra grafica
H = function(x,y) ifelse(y<=x & x>=0 & x<=1 & y>=0, x*y, 0)
x1=-0.1; x2=1.1; y1=-0.1; y2=1.1; z1=0; z2=1
x = seq(x1,x2,len=31); y = seq(y1,y2,len=31); z = outer(x,y,H)
NW(); MARG(0.3); th=60; ph=15 # poi anche con th=140
persp(x,y,z,zlim=c(z1,z2), theta=th, phi=ph, col="cyan",
ticktype="detailed",cex.axis=0.75,cex.lab=0.8,d=2)
# Dal grafico capisco che 0.125, ossia 1/8, è attendibile come risposta.
Se introduco in WolframAlpha il seguente comando ottengo 1/8:
integrate x*y dx dy, x=0..y, y=0..1
(B) Procediamo prima con R:
K = function(x,y) sin(x)+cos(y)
u1 = INTEGRAL(K, 0,pi/2, 0,pi/2, 200); u1
# 3.141601
u2 = INTEGRAL(K, 0,pi/2, 0,pi/2, 400); u2; u1-u2
# 3.141595 6.055927e-06
u3 = INTEGRAL(K, 0,pi/2, 0,pi/2, 800); u3; u2-u3
# 3.141593 1.513979e-06
u4 = INTEGRAL(K, 0,pi/2, 0,pi/2, 1600); u4; u3-u4
# 3.141593 3.784946e-07
# Capisco che è π, comunque ...
# Vedo che la differenza man mano si divide per 4. Posso prendere:
D=1/4; S=D; for(i in 1:100) {D=D/4; S=S+D}; S
# 0.3333333 1/3
L = u4 - (u3-u4)/3; more(L)
# 3.14159265358988
more(pi)
# 3.14159265358979
# Grafico (da cui capisco che il volume è circa 1.5^3, ossia circa 3):
x1=0; x2=pi/2; y1=0; y2=pi/2; z1=0; z2=2
x = seq(x1,x2,len=21); y = seq(y1,y2,len=21); z = outer(x,y,K)
NW(); MARG(0.3); th=60; ph=15
G = persp(x,y,z,zlim=c(z1,z2), theta=th, phi=ph, col="cyan",
ticktype="detailed",cex.axis=0.75,cex.lab=0.8,d=2)
lines (trans3d(x=c(x1,x2),y=0,z=0,pmat=G), col="red",lwd=2,lty=2)
lines (trans3d(y=c(y1,y2),x=0,z=0,pmat=G), col="red",lwd=2,lty=2)
lines (trans3d(z=c(z1,z2),x=0,y=0,pmat=G), col="red",lwd=2,lty=2)
Con WolframAlpha ottengo π con uno dei seguenti due comandi:
integrate sin(x)+cos(y) dx dy, x=0..pi/2, y=0..pi/2
integrate (integrate sin(x)+cos(y), x=0..pi/2), y=0..pi/2
(C) Procediamo prima con R. Vediamo inzialmente il grafico, sia del dominio che della superficie:
# Abbiamo aggiunto anche delle curve per chiarire il contorno
W = function(x,y) ifelse(x>=0 & x<=1 & y>=0 & y<=1 & y^2<=x & y>=x^2, x*y^2, 0)
x1=0; x2=1; y1=0; y2=1; z1=0; z2=1.1
x = seq(x1,x2,len=31); y = seq(y1,y2,len=31); z = outer(x,y,W)
NW(); MARG(0.3); th=0; ph=15
G = persp(x,y,z,zlim=c(z1,z2), theta=th, phi=ph, col="cyan",
ticktype="detailed",cex.axis=0.75,cex.lab=0.8,d=2)
# Per vedere meglio il grafico traccio anche le curve che lo delimitano:
t = seq(0,1, len=201); xr=t; yr=xr^2
lines (trans3d(xr, yr, 0, pmat = G), col="red",lwd=2,lty=2)
lines (trans3d(xr, yr, W(xr,yr), pmat=G),col="brown",lwd=3)
t=seq(0,1,len=201); yr=t; xr=yr^2
lines (trans3d(xr, yr, 0, pmat = G), col="red",lwd=2,lty=2)
lines (trans3d(xr, yr, W(xr,yr), pmat = G), col="brown",lwd=3)
#
# Faccio i calcoli (dovrà venire una piccola frazione di 1). Posso trascurare
# come si comporta la funzione fuori dal "quadrato":
Q = function(x,y) ifelse(y^2<=x & y>=x^2, x*y^2, 0)
u1 = INTEGRAL(Q, 0,1, 0,1, 250); more(u1)
# 0.0535951876336641
u2 = INTEGRAL(Q, 0,1, 0,1, 500); more(u2); u1-u2
# 0.0535620065396321 3.318109e-05
u3 = INTEGRAL(Q, 0,1, 0,1, 1000); more(u3); u2-u3
# 0.0535715169843706 -9.510445e-06
# Se voglio mi fermo qui e prendo 0.05357 come approssimazione, se no:
u4 = INTEGRAL(Q, 0,1, 0,1, 2000); more(u4); u3-u4
# 0.0535723091377355 -7.921534e-07
u5 = INTEGRAL(Q, 0,1, 0,1, 4000); more(u5); u4-u5
# 0.0535716354478788 6.736899e-07
u6 = INTEGRAL(Q, 0,1, 0,1, 8000); more(u6); u5-u6
# 0.0535715329191625 1.025287e-07
# Vediamo se è una "frazione":
fraction(u6)
# 3/56
more(3/56)
# 0.0535714285714286
Posso supporre che il risultato esatto sia 3/56. Avendo tempo a disposizione potrei
avvolare l'ipotesi con N maggiore. Controlliamo con WolframAlpha: col seguente comando ottengo 3/56
integrate x*y^2 dy dx, y=0..sqrt(x), x=0..1 - integrate x*y^2 dy dx, y=0..x^2,x=0..1
(ho fatto la differenza tra l'integrale sulla superficie fra asse x e y=x^2 e quello sulla superficie fra asse x e y=√x)