source("http://macosa.dima.unige.it/r.R") # If I have not already loaded the library ---------- ---------- ---------- ---------- ---------- ---------- ---------- ---------- Volumes of solids of revolution Rotation of a plane curve around a straight line that lies on the same plane [disc integration which integrates along the axis of revolution] An example (ellipsoid): BF=3; HF=2.5 f = function(x) sqrt(1-(x/2)^2)*1.5 PLANE(-2,2, -1.5, 1.5) g = function(x) -f(x) graph(f, -2,2, "brown"); graph1(g, -2,2, "brown") segm(0.6, f(0.6), 0.6,0, "seagreen") line(0.6,-f(0.6), 0.6,0, "seagreen") H = function(x) f(x)^2*pi integral(H, -2,2) # 18.84956 integral(H,-2,2)/pi # 6 The volume is 6*π Shell integration, as if the solid was composed of many small cylindrical shells An example: volume of a kind of "glass": BF=2.5; HF=2.5 PLANE(-2.5,2.5, 0,5) f = function(x) x^2; graph2(f, -2,2, "brown") F <- function(x,y) x^2+y^2-1; multrarot(F, 2,1/3, 0,4, 0, "brown") line(0,0, -2,4, "brown"); line(0,0, 2,4, "brown") # cylindrical shells of radius x, height 2x-x^2, thickness dx H = function(x) 2*pi*x*(2*x-x^2) integral(H, 0,2) # 8.37758 fraction(integral(H, 0,2)/pi) # 8/3 The volume is 8*π/3 The volume of a barrel (with parabolic section) D = 7.9; d = 6.65; h = 9.5 # dm # y=a*x^2+c fromx=-h/2, y=-d/2 to x=h/2,y=d/2 passing through x=0,y=D/2. # x=0, y=D/2 → c=D/2; x=h/2, y=d/2 → a*(h/2)^2+D/2 = d/2 → a=(d/2-D/2)/(h/2)^2 # y = (d/2-D/2)/(h/2)^2*x^2 + D/2 f = function(x) (d/2-D/2)/(h/2)^2*x^2+D/2 g = function(x) f(x)^2*pi; integral(g, -h/2,h/2) # 418.8702 (419 liters) # It can be find that: V = function(h,D,d) pi*h*(3*d^2 + 4*d*D + 8*D^2)/60 # verification: V(9.5, 7.9, 6.65) # 418.8702 OK # # The figure: BF=4; HF=4; BOXW(0,0,0,0) # I use BOXW to reduce white space x = c(-7.9/2,7.9/2); y = c(-7.9/2,7.9/2) ; w = c(0,9.5) # the box z = array(c(w[1],w[1],w[1],w[1]), dim=c(2,2)) th = 120; ph = 40 F = persp(x,y,z, theta=th,phi=ph, scale=FALSE, xlim=x,ylim=y,zlim=w,d=3,border="brown") # border="white", box=FALSE and co="white", if I don't want the box co="red" axes = function(F) { lines(trans3d(c(0,0),c(0,0),c(0,w[2]),pmat=F),col=co) lines(trans3d(c(0,x[2]),c(0,0),c(0,0),pmat=F),col=co) lines(trans3d(c(0,0),c(0,y[2]),c(0,0),pmat=F),col=co) lines(trans3d(c(0,0),c(0,0),c(w[1],0),pmat=F),col=co,lty=3) lines(trans3d(c(x[1],0),c(0,0),c(0,0),pmat=F),col=co,lty=3) lines(trans3d(c(0,0),c(y[1],0),c(0,0),pmat=F),col=co,lty=3) } axes(F) t=seq(0,9.5,0.1) lines(trans3d(0,f(t-9.5/2),t,pmat=F),col="black") for(a in seq(0,360,360/20)) lines(trans3d(cos(a*pi/180)*f(t-9.5/2),sin(a*pi/180)*f(t-9.5/2),t,pmat=F),col="black") a = (0:50)*2*pi/50 lines(trans3d(cos(a)*d/2,sin(a)*d/2,0,pmat=F),col="blue") lines(trans3d(cos(a)*d/2,sin(a)*d/2,9.5,pmat=F),col="blue") lines(trans3d(cos(a)*D/2,sin(a)*D/2,9.5/2,pmat=F),col="blue") Other examples of use