1623.sagews

This is an example done in class on Wed., Feb. 17.  I've added colors to the planes.  The coordinate planes are now green.  The actual surface through which we are calculating flux is contained in the red plane.

var('x y z')
p1 = parametric_plot3d((x, y, 1+2*y-x), (x,-1,2), (y,-1,2), opacity=.65, color='red')
p2 = parametric_plot3d((x,y,0), (x,-1,2), (y,-1,2), opacity=.65, color='green')
p3 = parametric_plot3d((x,0,z), (x,-1,2), (z,-2,4), opacity=.65, color='green')
p4 = parametric_plot3d((0,y,z), (y,-1,2), (z,-2,4), opacity=.65, color='green')
p5 = parametric_plot3d((x,1-x,z), (x,-1,2), (z,-2,4), opacity=.65, color='blue')
show(p1+p2+p3+p4+p5)

var('x y z u v')
F = vector([y,z,x^2])
r = vector([u,v,1+2*v-u])

integrand(u,v) = F(x=r[0], y=r[1], z=r[2]) * diff(r,u).cross_product(diff(r,v))

The first way of computing the flux across our surface mentioned in an email I sent.

numerical_integral(lambda u: numerical_integral(integrand(u,v), 0, 1-u)[0], 0, 1)

(-1.0833333333333335, 1.2155813821184732e-14)

The other method mentioned in the email, employing the Python scipy package.  Several quirks of this method: 1) it requires the parameter that appears first in integrand to be the inside integral.  2) The limits of integration for the inner integral come last (below I've called them a and b), and it seems they need to be functions---using numbers appears not to work.

import scipy
from scipy import integrate
integrand(v,u) = F(x=r[0], y=r[1], z=r[2]) * diff(r,u).cross_product(diff(r,v))
a(x)=0.0
b(x)=1-x
integrate.dblquad(integrand,0,1,a,b)

(-1.0833333333333335, 1.2159275043597867e-14)