source("http://macosa.dima.unige.it/r.R")    # If I have not already loaded the library
---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
Integrals of two variables functions
 
# For the integral of a function of two variables (f in [a,b]×[c,d])
# we can use the command  INTEGRAL(f, a,b, c,d, N)  that calculates
# the "volume" of a N×N parallelepiped whose upper part approximates
# the surface of f.
     
  
# When the graph of the function is a plane INTEGRAL approximates very quickly. # An example: f = function(x,y) 5 INTEGRAL(f, 0,5, 0,5, 1) # 125 INTEGRAL(f, 0,5, 0,5, 10) # 125 # Another (the volume of the part under z=0 equalize the one above): g = function(x,y) (y-x)/2 INTEGRAL(g, -5,5, -5,5, 1) # 0 INTEGRAL(f, -5,5, -5,5, 10) # 0 Other examples: # [1] A bivariate gaussian (just a small N value is enough): F = function(x,y) 1/(2*pi)*exp(-(x^2+y^2)/2) u = INTEGRAL(F,-20,20, -20,20, 10); u # 0.1865616 u = INTEGRAL(F,-20,20, -20,20, 20); u # 0.9714394 u =INTEGRAL(F,-20,20, -20,20, 50); u # 1 # # [2] The half-sphere of radius 1 (I put 0 out of the domain): g = function(x,y) ifelse(x^2+y^2>1, 0, sqrt(1-x^2-y^2) ) u = INTEGRAL(g,-1,1, -1,1, 500); u # 2.094412 u = INTEGRAL(g,-1,1, -1,1, 1000); u # 2.094397 u/pi # 0.6666674 fraction(0.6666666666666) # 2/3 The integral is 2/3·π # # [3] The function in the 3rd graph (I put 0 out of the domain): F = function(x,y) ifelse(y > 1-x^2, 0, 1-x+y) u = INTEGRAL(F,-1,1, 0,1, 250); more(u) # 1.86677023999999 u = INTEGR(F,-1,1, 0,1, 500); more(u) # 1.86705564800002 u = INTEGR(F,-1,1, 0,1, 1000); more(u) # 1.86669718 fraction(1.866666666666666) # 28/15 # Using WolframAlpha: # integrate(integrate 1-x+y, y =0..1-x^2), x = -1..1 # 28/15 # # I N is very large, approximation errors can become preponderant. # # Another example: integral of (x,y) -> x*y^2 in the following figure # between y=x^2 and x=y^2
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) # The curves that delimit the graph: 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)
# I understand that the integral is much less than 1. # In the square [0,1]×[0,1] where I have to do the calculation I can simplify # the function: 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 # I can stop and take 0.05357. Otherwise: 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 # Is it a "fraction"? fraction(u6) # 3/56 more(3/56) # 0.0535714285714286 # If I want, I can use WolframAlpha and verify that 3/56 is actually the exact result # # For the graph of two variables functions see S19. Other examples of use