#! /usr/bin/env python # -*- coding: utf-8 -*- """ Plots and computes the volume delimited by two surfaces (in 3D). @date: Fri Feb 06 11:10:46 2015. @author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015. @licence: GNU Public License version 3. """ # 3D plotting of the two surfaces # We import the pylab environment. It is like calling %pylab in IPython from pylab import * zMax = 36 # WARNING cannot be modified, the R is dependant on zMax R = 3 n = 25 x, y = linspace(-R, R, n), linspace(-R, R, n) X, Y = meshgrid(x, y) print "Shape of X:", X.shape print "Shape of Y:", Y.shape # %% Plot the intersection circle # A 3D plot for two volumes # Comes from Exercise 24.a) of the Tutorial Sheet 2 for MA102 fig = figure(14) from mpl_toolkits.mplot3d import Axes3D ax = Axes3D(fig) title("Volume limited by two paraboloids.") # We first draw the intersection circle nt = 10*n t = linspace(0, 2*pi, nt) xt, yt, zt = R*cos(t), R*sin(t), array([9]*nt) ax.plot(xt, yt, zt, ':', color='black', linewidth=2, alpha=0.5) ax.legend("Circle of intersection ($R{}, z={}$).".format(R, R**2)) # find a way to get these as 2D indexes circleBool = X**2 + Y**2 <= R**2 print "Shape of that boolean array for being inside or not of the circle C(O, 3):", circleBool.shape indexLimited_find = find(circleBool) print "Shape of that boolean array with find() function:", indexLimited_find.shape indexes = [(k % n, k / n) for k in indexLimited_find] X_Limited = array([X[i, j] for i, j in indexes]) Y_Limited = array([Y[i, j] for i, j in indexes]) # %% Next step # Lower "egg" volume Z1 = X**2 + Y**2 Z1[X**2 + Y**2 > R**2] = R**2 # FIXME Z1_Limited = X_Limited**2 + Y_Limited**2 #ax.plot(X_Limited, Y_Limited, Z1_Limited, alpha=0.8) #ax.plot(X_Limited, Y_Limited, Z1_Limited, '+', alpha=0.8) #surf1 = ax.plot_surface(X, Y, -Z1, rstride=1, cstride=1, cmap='coolwarm', alpha=1.0, linewidth=0) trisurf1 = ax.plot_trisurf(X_Limited, Y_Limited, Z1_Limited, cmap='coolwarm', alpha=1.0, linewidth=0) ax.legend("Lower paraboloid ($z = x^2 + y^2$).") #fig.colorbar(surf1, shrink=0.5, aspect=5) fig.colorbar(trisurf1, shrink=0.5, aspect=5) #ax.contourf(X, Y, Z1, zdir='z', offset=0, cmap='winter') # Upper "egg" volume Z2 = zMax - 3*(X**2 + Y**2) Z2[X**2 + Y**2 > R**2] = R**2 # FIXME Z2_Limited = zMax - 3*(X_Limited**2 + Y_Limited**2) #ax.plot(X_Limited, Y_Limited, Z2_Limited, alpha=0.8) #ax.plot(X_Limited, Y_Limited, Z2_Limited, '*', alpha=0.8) #surf2 = ax.plot_surface(X, Y, -Z2, rstride=1, cstride=1, cmap='hot', alpha=0.8, linewidth=0) trisurf2 = ax.plot_trisurf(X_Limited, Y_Limited, Z2_Limited, cmap='hot', alpha=0.8, linewidth=0) ax.legend("Upper paraboloid ($z = {} - {}(x^2 + y^2$)).".format(zMax, R)) #fig.colorbar(surf2, shrink=0.5, aspect=5) fig.colorbar(trisurf2, shrink=0.5, aspect=5) #ax.contourf(X, Y, Z2, zdir='z', offset=0, cmap='hot') # Conclude ax.set_zlim(0, zMax) #ax.set_zlim(-zMax, 0) #savefig('DemoTwoVolumes__figure1.svg') savefig('Demo_Two_Volumes.png') # %% Asking SymPy some help to compute the volume raw_input("Continue ? (SymPy will be used to formally compute the volume).") from sympy import * from sympy.abc import x, y, z init_printing(use_unicode=False, wrap_line=False, no_global=True) #f1 = lambdify((x, y), x**2 + y**2) #f2 = lambdify((x, y), zMax - 3*(x**2 + y**2)) #yBound = lambdify(x, sqrt(9 - x**2)) f1 = lambda x, y: x**2 + y**2 f2 = lambda x, y: zMax - 3*(x**2 + y**2) yBound = lambda x: sqrt(9 - x**2) # Compute the triple-integral, first on x, then on y, then on z print "\nThat triple-integral is first on $x$, then on $y$, then on $z$:" print r"\iiint_{\mathcal{D}} 1 dx~dy~dz $\int_{x=-R}^{x=R} \int_{y=-\sqrt{9 - x^2}}^{y=\sqrt{9 - x^2}} \int_{z=x^2+y^2}^{z=36-3(x^2+y^2)} 1 dx~dy~dz$" # Formal input s_xyz = Integral(1, (z, f1(x, y), f2(x, y)), (y, -yBound(x), yBound(x)), (x, -R, R)) print "SymPy represents formally this as:", s_xyz, "\nAnd it can be sexy:" pprint(s_xyz) # %% Formal computation: let integrate it! print "But SymPy can also compute it formally." yBound = lambda x: (9 - x**2)**0.5 #intf_x_y_z = integrate( # integrate( # integrate( # 1.0, # (z, f1(x, y), f2(x, y)) # ), # (y, -yBound(x), yBound(x)) # ), # (x, -R, R) #) # #print "And that integral can be printed (but might be ugly):\n", intf_x_y_z # Let us help SymPy a little bit by computing step by step, from the inner first to the last one f_z = lambdify((x, y), integrate(1.0, (z, f1(x, y), f2(x, y)))) f_y = lambdify(x, integrate(f_z(x, y), (y, -yBound(x), yBound(x)))) intf_x = integrate(f_y(x), (x, -R, R)) print "And that integral can be printed (but might be ugly):\n", intf_x try: # value = float(intf_x_y_z) value = float(intf_x) print "And we can compute its value: {:g}.".format(value) print "(More precisely: {:.20g})".format(value) except Exception as e: print "Error:", e print "And it is too ulgy to compute its approximate value!" # %% Final test print "Finally, that volume is :" pprint(s_xyz) print "Which can be computed and is almost equal to:" pprint(value)