︠618cb425-0f96-4cb2-9b71-06d566528c7e︠
################################
#
# TABLE OF CONTENTS
# (You can move to any line by using the "lightning bolt" button to the right of the large 'A')
#
# 65 Distributive property for dot- and cross-product
# 112 Miscellaneous examples, done in class on R 8/28/18
# 125 Plot a plane, given three points on it
# 170 Plot a plane, given a normal vector and the distance from the plane to the origin
# 220 Miscelleaneous computations involving vectors
# 240 Interactive 3D grapher: VERY USEFUL!!
# 340 Solution to the extra credit for quiz 1
# 391 Plotting parametric curves with normal, tangent, binormal vectors,
# osculating plane, acceleration, curvature, etc.
# 542 Interactive TANGENT PLANE PLOTTER. Includes partial derivatives. Very useful!
# 645 Surface plot compared with two different styles of contour plot. Includes a "fancy" surface plot that is prettier than
# the standard one, but you cannot tilt it with the mouse. But is is designed to exactly match the
# contour plot.
# 680 Chain-rule practice. Great homework helper!
# 715 Midpoint approximation and exact computation of simple iterated integrals. Another good homework helper!
# 755 Problems 15.6.33-34
# 885 Three views of the element of volume for spherical coordinates
# 1022 Images from quiz #5 (examples of cylindrical and spherical plotting)
# 1075 Interactive (u,v) --> (x,y) transformation explorer (Section 15.9). Very instructive! Shows Jacobians, too!
# 1120 Examples of plotting vector fields
# 1158 Examples of parametric equations of surfaces
# 1175 INTERACTIVE Parametric Surface EXPLORER
#ILLUSTRATION OF DISTRIBUTIVE PROPERTY FOR DOT PRODUCT AND FOR CROSS PRODUCT
#In the diagram created, vector A is a unit vector. We draw both parallel and perpendicular projections of vectors B and C onto A.
# The main idea is that given any vectors B and A, we can resolve B into two orthogonal components: the PARALLEL component of B onto A is parallel to A, while the
# PERPENDICULAR component of B onto A is perpendicular to A and in the same plane as A and B.
# The MAGNITUDE of the PARALLEL projection of B onto A is the absolute value of the dot-product of B and A (assume A is a unit vector),
# and the PERPENDICULAR projection of B onto A (the actual VECTOR, not the magnitude) is the cross-product of B and A.
def perpProj(x,y):
#output is the perpendicular component of x wrt y, where
#it is understood for convenience that y is a unit vector
#and x, y are both vectors
return x-y*x.dot_product(y)
def parProj(x,y):
#output is the parallel component of x wrt y, where y is unit
return y*x.dot_product(y)
a=vector([0,1,0])
b=vector([-.3,.1,.7])
c=vector([-.7,.5,.3])
t=6 #thickness of some vectors
bc=b.cross_product(a)
P=plot(a, color='red', thickness=t)
P += plot(b, color='blue', thickness=t)
P += plot(c, color='black',thickness =t, start=b)
#P += plot(c,color='black',start=b)
#P += plot(b,color='blue',start=c)
P += plot(text3d("B",b*0.7-a*.05,color='blue'))
P += plot(text3d("A",a*1.05,color='red'))
P += plot(text3d("C",b*1.2+c*0.3,color='black'))
P += plot(b+c,color="magenta", thickness=t)
P += plot(text3d("B+C", (b+c)*.6, color='magenta'))
P += plot(-perpProj(b,a), color='blue',start=b)
P += plot(-perpProj(b,a), color='blue',start=b+parProj(c,a))
P += plot(parProj(b,a), color='blue',start=-0.05*perpProj(b,a))
P += plot(parProj(c,a), color='black',start = parProj(b,a)-0.08*perpProj(b,a))
P += plot(parProj(b+c,a), color='magenta',start=-0.12*perpProj(b,a))
P += plot(-perpProj(b+c,a), color='magenta',start=b+c)
P += plot(parProj(c,a),color='black',start=b)
P += plot(perpProj(c,a), color='black',start=b+parProj(c,a))
#P += plot(parProj(b+c,a), color='magenta',start=perpProj(b+c,a))
P.show(aspect_ratio=1.0,frame=false, spin=5)
︡934b104e-4b99-475f-bd73-ef0a862cc034︡{"file":{"filename":"142cacbc-a1ca-4de5-803d-4b8fe7198fab.sage3d","uuid":"142cacbc-a1ca-4de5-803d-4b8fe7198fab"}}︡{"done":true}
︠6132de36-8a02-46c5-9859-c33a18fb87c1si︠
78^52
︡98d62e8d-f427-4fad-b527-b98e6d8bbc3e︡{"stdout":"244860842756831203896798793038363593764863637715027581243662187665028113084218676375430402110652416\n"}︡{"done":true}
︠322c7941-29cf-4a87-ba7d-eaf4bd125ca6︠
plot(sin,0,pi)
︡4ee91459-f0db-432b-8582-45fcdd3256d1︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/554/tmp_UQQ26l.svg","show":true,"text":null,"uuid":"77b262f0-1d8e-43b8-ad8c-08cde07bf3a5"},"once":false}︡{"done":true}︡
︠072cf6cd-7ac8-44be-81ce-2723453f5418︠
%var x,y,z
pic=implicit_plot3d(x^2+y^2+z^2-1,(x,-1,0.5),(y,-1,1),(z,-1,1),color='red',mesh=1)
pic.show()
︡59ea9731-75e2-4161-b67d-002823487089︡{"file":{"filename":"c89cf60c-4792-45fb-b148-14d5d10b6cc5.sage3d","uuid":"c89cf60c-4792-45fb-b148-14d5d10b6cc5"}}︡{"done":true}︡
︠eea61aea-e2cf-46b0-8724-e80d2520064d︠
#Plot a plane, given three points on the plane
# try changing the points defined in line 27 of this
#declare variables
%var x, y, z
#define a function to draw axes. Very useful for understanding plots
def axes(size,xcolor='red',ycolor='green',zcolor='blue'):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= size; zSize=size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
asize = 16 #axis size
psize = 12 #plane size
#define the three points
a= (0,0,1); b= (0,1,0); c= (1,0,0)
#turn em into vectors
A = vector(a); B= vector(b); C=vector(c)
AB = B-A; AC = C-A
N = AB.cross_product(AC) #note the syntax for cross product
f = N[0]*x+N[1]*y + N[2]*z-N.dot_product(A) #set this equal to zero to get equation of plane
P = plot(axes(asize))
P += implicit_plot3d(f, (x,-psize, psize), (y,-psize, psize), (z, -psize, psize), color='orange', opacity=0.1)
P += plot(points([a,b,c], thickness=40, color='black'))
P += plot(N, start = a, width = 3, color='black')
show('The equation of the plane is $'+latex(f)+'=0$')
show('Normal vector is black.')
show(P, spin =5)
︡624c3712-c208-4e4b-b6a5-df425b7ad5db︡{"html":"
The equation of the plane is $ -x - y - z + 1 =0$
"}︡{"html":"Normal vector is black.
"}︡{"file":{"filename":"d2774532-5c4f-4000-b637-1f8d24701ace.sage3d","uuid":"d2774532-5c4f-4000-b637-1f8d24701ace"}}︡{"done":true}︡
︠b9407e1e-c750-404c-8168-499c2ba97bb6︠
#Plot a plane, given a normal vector and the distance from the origin
#declare variables
%var x, y, z
#define a function to draw axes. Very useful for understanding plots
def axes(size,xcolor='red',ycolor='green',zcolor='blue'):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= size; zSize=size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
asize = 16 #axis size
psize = 12 #plane size
#define the normal and distance to origin
N=vector((1,1,1))
d=5 #dist to origin
#equation of plane is K*=d ONLY as long as K is a unit normal vector
K=N/N.norm() #note syntax for norm of a vector: the norm of A is A.norm()
f = K[0]*x+K[1]*y + K[2]*z-d #set this equal to zero to get equation of plane
P = plot(axes(asize))
P += implicit_plot3d(f, (x,-psize, psize), (y,-psize, psize), (z, -psize, psize), color='orange', opacity=0.1)
P += plot(10*N, width = 3, color='black') #I multiplied N to make it longer. You may want to change this.
show('The equation of the plane is $'+latex(f)+'=0$')
show('Normal vector is black')
show(P, spin =5)
︡49b85d03-07ca-48bd-b9c1-340f0f29cad1︡{"html":"The equation of the plane is $ \\frac{1}{3} \\, \\sqrt{3} x + \\frac{1}{3} \\, \\sqrt{3} y + \\frac{1}{3} \\, \\sqrt{3} z - 5 =0$
"}︡{"html":"Normal vector is black
"}︡{"file":{"filename":"3e970b4f-a1cf-4ca1-be5e-909f2708298f.sage3d","uuid":"3e970b4f-a1cf-4ca1-be5e-909f2708298f"}}︡{"done":true}︡
︠45f1ac24-2169-4759-a5ad-0df5604400ad︠
︡4c165fef-faa7-4194-91c1-0af67d892b43︡
︠1f4c2f75-2e84-414e-89a5-f282cebda96a︠
a=(2,-1,5)
A=vector(a)
b=(1,1,1); c=(0,-2,4)
B=vector(b);C=vector(c)
A.dot_product(B)
︡4d33c433-975e-43c0-8fb8-c0aaf0f7f71f︡{"stdout":"6\n"}︡{"done":true}︡
︠838140c6-728d-41b4-871e-edd4deca4da3︠
A.cross_product(B)
︡e97f9d6a-fe62-4688-a75c-785ffbfb354b︡{"stdout":"(-6, 3, 3)\n"}︡{"done":true}︡
︠92fdff67-3101-43c0-9c5b-c79a3f801fa4︠
B.cross_product(A)
︡c3223560-19fb-4b0e-8ab4-f825770b4743︡{"stdout":"(6, -3, -3)\n"}︡{"done":true}︡
︠22da9dd3-e93f-4eae-ba3f-f17646962485︠
A.cross_product(B)
︡1dc872e7-ea93-493c-a192-d3a914363b07︡{"stdout":"(-6, 3, 3)\n"}︡{"done":true}︡
︠288a8f70-4268-4045-afcf-dc7b72792c5a︠
A.norm()
︡bca446be-4199-407d-a431-6deb7a6374fb︡{"stdout":"sqrt(30)\n"}︡{"done":true}︡
︠bee76581-a9dc-4abb-9200-a0b558d7cf7fs︠
################################################################################
#
# INTERACTIVE 3D GRAPHER
#
# This allows you to graph up to 3 three-dimensional plots at once, choosing whether they
# will be parametric, implicit, or z=f(x,y) style. You don't have to do any coding, but instead just fill in
# the input boxes. It **should** be self-explanatory and easy to use.
#
################################################################################
%var t,x,y,z
def axes(size,xcolor='red',ycolor='green',zcolor='blue',asp=(1,1)):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= asp[0]*size; zSize=asp[1]*size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
#return xAxis+yAxis+zAxis
###############################################################################
def plotPlane(N,P,s,clr='grey',op=0.2,meshvalue=0):
x,y,z=var('x y z')
#N=normal vec, P = point (vector), s=size
#This returns a plot of plane with the given normal vector and point, with x,y,z going plus or minus s units from the point.
#optional arguments for color and opacity (default is nearly transparent grey)
return implicit_plot3d(N.dot_product(vector((x,y,z)))-N.dot_product(P),(x,P[0]-s,P[0]+s),(y,P[1]-s,P[1]+s),(z,P[2]-s,P[2]+s),color=clr,opacity=op,mesh=meshvalue)
#default values
color1='red';color2='lightgreen';color3='lightblue'
paraf=(sin(t),cos(t),t)
implf=x^2+y^2-z^2-4
zfxy =x^2-y^2
tlims =(0,4*pi);xlims = (-3,3);ylims=(-3,3);zlims=(-3,3)
html('''f1 is RED f2 is GREEN f3 is BLUE''')
@interact(layout=dict(top=[['f1Show','f1Type','eqn1','tx'],['t1','x1','y1','z1'],['f2Show','f2Type','eqn2','tx2'],['t2','x2','y2','z2'],['f3Show','f3Type','eqn3','tx3'],['t3','x3','y3','z3']]))
def _(f1Show=checkbox(false,label= 'f1Show?'),
f1Type=selector(['parametric','implicit (mesh)','z=f(x,y) (mesh)'],width=10),
eqn1=input_box(paraf,width=50),
tx=text_control(' ',label=''),
t1=input_box(tlims,width=20,label='t1vals'),x1=input_box(xlims,width=20,label="x1vals"),y1=input_box(ylims,width=18,label="y1vals"),z1=input_box(zlims,width=18,label="z1vals"),
f2Show=checkbox(true,label= 'f2Show?'),
f2Type=selector(['parametric','implicit (mesh)','z=f(x,y) (mesh)'],default='implicit (mesh)',width=10),
eqn2=input_box(implf,width=50),
tx2=text_control(' ',label=''),
t2=input_box(tlims,width=20,label='t2vals'),x2=input_box(xlims,width=20,label="x2vals"),y2=input_box(ylims,width=18,label="y2vals"),z2=input_box(zlims,width=18,label="z2vals"),
f3Show=checkbox(false,label= 'f3Show?'),
f3Type=selector(['parametric','implicit (mesh)','z=f(x,y) (mesh)'],default='z=f(x,y) (mesh)',width=10),
eqn3=input_box(zfxy,width=50),
tx3=text_control(' ',label=''),
t3=input_box(tlims,width=20,label='t3vals'),x3=input_box(xlims,width=20,label="x3vals"),y3=input_box(ylims,width=18,label="y3vals"),z3=input_box(zlims,width=18,label="z3vals")):
P1=Graphics();P2=Graphics(); P3=Graphics()
if f1Show:
if f1Type=='parametric':
P1=parametric_plot(eqn1,(t,t1[0],t1[1]),color=color1,plot_points=1000,thickness=4)
elif f1Type=='implicit (mesh)':
P1=implicit_plot3d(eqn1,(x,x1[0],x1[1]),(y,y1[0],y1[1]),(z,z1[0],z1[1]),mesh=1,color=color1)
elif f1Type=='z=f(x,y) (mesh)':
P1=plot3d(eqn1,(x,x1[0],x1[1]),(y,y1[0],y1[1]),mesh=1,color=color1)
if f2Show:
if f2Type=='parametric':
P2=parametric_plot(eqn2,(t,t2[0],t2[1]),color=color2,plot_points=1000,thickness=4)
elif f2Type=='implicit (mesh)':
P2=implicit_plot3d(eqn2,(x,x2[0],x2[1]),(y,y2[0],y2[1]),(z,z2[0],z2[1]),mesh=1,color=color2)
elif f2Type=='z=f(x,y) (mesh)':
P2=plot3d(eqn2,(x,x2[0],x2[1]),(y,y2[0],y2[1]),mesh=1,color=color2)
if f3Show:
if f3Type=='parametric':
P3=parametric_plot(eqn3,(t,t3[0],t3[1]),color=color3,plot_points=1000,thickness=4)
elif f3Type=='implicit (mesh)':
P3=implicit_plot3d(eqn3,(x,x3[0],x3[1]),(y,y3[0],y3[1]),(z,z3[0],z2[1]),mesh=1,color=color3)
elif f3Type=='z=f(x,y) (mesh)':
P3=plot3d(eqn3,(x,x3[0],x3[1]),(y,y3[0],y3[1]),mesh=1,color=color3)
#axes set up. need size, but don't use limits if not shown
xsize=max(f1Show*x1[1],f2Show*x2[1],f3Show*x3[1])
ysize=max(f1Show*y1[1],f2Show*y2[1],f3Show*y3[1])
zsize=max(f1Show*z1[1],f2Show*z2[1],f3Show*z3[1])
ax = axes(xsize,asp=[ysize/xsize, zsize/xsize])
show(P1+P2+P3+ax, spin =2)
︡c1c9078d-babf-473d-8b68-090789bba1c0︡{"html":"f1 is RED f2 is GREEN f3 is BLUE"}︡{"interact":{"controls":[{"control_type":"checkbox","default":false,"label":"f1Show?","readonly":false,"var":"f1Show"},{"button_classes":null,"buttons":false,"control_type":"selector","default":0,"label":"f1Type","lbls":["parametric","implicit (mesh)","z=f(x,y) (mesh)"],"ncols":null,"nrows":null,"var":"f1Type","width":10},{"control_type":"input-box","default":"(sin(t), cos(t), t)","label":"eqn1","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"eqn1","width":50},{"classes":null,"control_type":"text","default":" ","label":"","var":"tx"},{"control_type":"input-box","default":"(0, 4*pi)","label":"t1vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"t1","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"x1vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"x1","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"y1vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"y1","width":18},{"control_type":"input-box","default":"(-3, 3)","label":"z1vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"z1","width":18},{"control_type":"checkbox","default":true,"label":"f2Show?","readonly":false,"var":"f2Show"},{"button_classes":null,"buttons":false,"control_type":"selector","default":1,"label":"f2Type","lbls":["parametric","implicit (mesh)","z=f(x,y) (mesh)"],"ncols":null,"nrows":null,"var":"f2Type","width":10},{"control_type":"input-box","default":"x^2 + y^2 - z^2 - 4","label":"eqn2","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"eqn2","width":50},{"classes":null,"control_type":"text","default":" ","label":"","var":"tx2"},{"control_type":"input-box","default":"(0, 4*pi)","label":"t2vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"t2","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"x2vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"x2","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"y2vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"y2","width":18},{"control_type":"input-box","default":"(-3, 3)","label":"z2vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"z2","width":18},{"control_type":"checkbox","default":false,"label":"f3Show?","readonly":false,"var":"f3Show"},{"button_classes":null,"buttons":false,"control_type":"selector","default":2,"label":"f3Type","lbls":["parametric","implicit (mesh)","z=f(x,y) (mesh)"],"ncols":null,"nrows":null,"var":"f3Type","width":10},{"control_type":"input-box","default":"x^2 - y^2","label":"eqn3","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"eqn3","width":50},{"classes":null,"control_type":"text","default":" ","label":"","var":"tx3"},{"control_type":"input-box","default":"(0, 4*pi)","label":"t3vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"t3","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"x3vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"x3","width":20},{"control_type":"input-box","default":"(-3, 3)","label":"y3vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"y3","width":18},{"control_type":"input-box","default":"(-3, 3)","label":"z3vals","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"z3","width":18}],"flicker":false,"id":"fd7c584d-f319-47f8-8a99-b0041218a59f","layout":[[["f1Show",3,null],["f1Type",3,null],["eqn1",3,null],["tx",3,null]],[["t1",3,null],["x1",3,null],["y1",3,null],["z1",3,null]],[["f2Show",3,null],["f2Type",3,null],["eqn2",3,null],["tx2",3,null]],[["t2",3,null],["x2",3,null],["y2",3,null],["z2",3,null]],[["f3Show",3,null],["f3Type",3,null],["eqn3",3,null],["tx3",3,null]],[["t3",3,null],["x3",3,null],["y3",3,null],["z3",3,null]],[["",12,null]]],"style":"None"}}︡{"done":true}
︠eeb41899-cba7-407c-bde7-d6f2bfcf02e7︠
#Solution to the extra credit for quiz 1
%var x, y, z
def axes(size,xcolor='red',ycolor='green',zcolor='blue'):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= size; zSize=size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
asize = 16
psize = 12
a= (1,0,-2); b= (3,1,-2); c= (-5,-1,0); d= (1,1,1)
A = vector(a); B= vector(b); C=vector(c); D= vector(d)
AB = B-A; AC = C-A
N = AB.cross_product(AC)
#equation of plane
f = N[0]*x+N[1]*y + N[2]*z-N.dot_product(A)
#the displacement vector which takes D to its reflection about the plane
displacement = N*2*(N.dot_product(A-D))/(N.dot_product(N))
P = plot(axes(asize))
P += implicit_plot3d(f, (x,-psize, psize), (y,-psize, psize), (z, -psize, psize), color='orange', opacity=0.1)
P += plot(points([a,b,c,d], thickness=40, color='black'))
P += plot(displacement, start = D, width = 3, color='black')
#display the coordinates of E and the amount of distance from D to E.
print D+displacement, displacement.norm()
show(P, spin =5)
︡f9826e14-9e82-4c94-97cb-9b896231e72b︡{"stdout":"(1/9, 25/9, -7/9) 8/3\n"}︡{"file":{"filename":"83f99185-6e32-45f9-9bfb-756439d8adc0.sage3d","uuid":"83f99185-6e32-45f9-9bfb-756439d8adc0"}}︡{"done":true}︡
︠356020a5-6fae-4af2-8a98-e638c3204749︠
############################################################################################
#
# An example of a parametric curve of a vector function, with tangent vectors, osculating plane, etc.
#
# ON LINES 69--72, You define a parametric curve r(t), and set the starting and ending t-values. You also include a
# value of t in the middle (the constant "c"), a display size for the osculating plane, and you will be rewarded with a parametric plot that shows the
# velocity, acceleration, unit tangent, and unit normal vectors, along with the osculating plane, and computes values for
# these vectors as well as the curvature.
#
# The acceleration vector (computed and displayed in lines 80--83, 112--116) is often too large to display nicely. Try experimenting with keeping it commented out
# (the default) and instead displaying a vector with magnitude 3 with the ACTUAL magnitude displayed, or commenting out this magnitude display (lines 113--116) and
# displaying the actual acceleration vector (line 112).
#
# THAT'S NOT ALL! You also get a graph showing speed, scalar acceleration, and curvature vs. t, which will help you to understand the important
# a = v'T + (kv^2)N equation (the decomposition of the acceleration vector into unit tangent and unit normal components)
#
#############################################################################################
# NOTE: This example uses functions that produce axes and planes functions defined above. So we will include them in our code here, starting on the next line
#########################
#
# AXES function
#
#########################
def axes(size,xcolor='red',ycolor='green',zcolor='blue',asp=(1,1)):
#returns a plot of axis vectors, labeled, with length equal to the 'size' argument
#Note the OPTIONAL color arguments. If you leave them out, the xyz axes are red, green, blue, respectively
#But, for example, the command axes(12,xcolor='black') would produce axes of length 12 with x in black and the others
#still in green and blue.
#also notice the OPTIONAL axp argument, which is two numbers to adjust the size of the y- and z- axis. For
#example, axes(12, asp=(2,.5)) will produce axes with the x-axis of length 12, the y-axis has length 25, and the z-axis is length 6.
xSize=size; ySize= asp[0]*size; zSize=asp[1]*size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to; not sure if this works)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
#return xAxis+yAxis+zAxis (if we didn't want labels)
###############################################################################
#
# PLANE drawing function
#
######################################
def plotPlane(N,P,s,clr='grey',op=0.2,meshvalue=0):
x,y,z=var('x y z')
#N=normal vec, P = point (vector), s=size
#This returns a plot of plane with the given normal vector and point, with x,y,z going plus or minus s units from the point.
#optional arguments for color and opacity and meshvalue (default is nearly transparent grey)
return implicit_plot3d(N.dot_product(vector((x,y,z)))-N.dot_product(P),(x,P[0]-s,P[0]+s),(y,P[1]-s,P[1]+s),(z,P[2]-s,P[2]+s),color=clr,opacity=op,mesh=meshvalue)
#
# end of function definitions
#
#####################################
#you need to define symbolic variables other than x, to do useful symbolic math like plotting and calculus
%var t
################# DEFINE YOUR r(t) PARAMETRIZATION BELOW, along with start and end t and intermediate t-value
#define some parametric functions for our curve
#this example uses a wiggly curve with variable curvature helix
x(t) = (4*cos(t)+0.5)*cos(t); y(t) = (4*cos(t)+0.5)*sin(t); z(t) = t
t0=0; t1= 4*pi #starting and ending t values for our curve
c= 0.3*pi #a t value in between to draw tangent and normal etc
size=1 #the size for the osculating plane. Experiment with different values
#########################################
#
# FROM THIS POINT ON, EVERYTHING IS AUTOMATED. READ THE CODE TO SEE HOW IT IS DONE
#
########################################
r(t) = (x(t),y(t),z(t)) #the position r(t) at time t
v=r.diff(t) #velocity! The "diff" suffix does symbolic differentiation; it is a function of t
a=r.diff(t,2) #acceleration (the "2" means second derivative); also a function of t
#Our convention is to use upper-case for vectors
R = vector(r); V= vector(v); A = vector(a)
#We need to turn them into vectors in order to use vector operations like cross-product, etc.
#Now compute unit tangent and unit normal vectors.
#We will compute T, N, B, and T' vectors from textbook
#recall that norm() computes magnitude of a vector
T = V/norm(V) #unit tangent vector (function of t)
Tprime = T.diff(t) #derivative of unit tangent; this is the T' vector in your textbook
#turn this into a unit normal
N = Tprime/norm(Tprime)
# Finally, the little-used binormal vector B
B=T.cross_product(N)
# So we have produced the following vectors, all functions of t:
# R, V, A, T, Tprime, N, B
#
# Let's also compute the curvature, using the formula k = |T'|/|ds/dt|=|T'|/|V|.
#curvature (two ways; should be equal)
K = norm(Tprime)/norm(V)
# This is a scalar function of t.
#create the parametric plot of curve r(t)
P = Graphics() # start with a blank slate
# Next, include the parametric plot
P += parametric_plot3d(r(t),(t,t0,t1),plot_points=400,thickness=6)
P += plot(T(t=c),start=R(t=c),color='red') #include the unit tangent at t=c; note starting point
P += plot(N(t=c),start=R(t=c),color='black') #add unit normal; note starting point
P += plot(B(t=c),start=R(t=c),color='green') #the binormal vector is normal to the osculating plane
# the next line allows you to see the acceleration vector, which, like T and N, lies on the osculating plane
# However, it is often too long, and you may want to comment out the next line or replace it with the line after it which displays a shortened vector but with the value of the magnitude labeled.
#P += plot(A(c),start=R(c),color='magenta') #the ACTUAL acceleration vector
P += plot(3.0*A(t=c)/norm(A(t=c)),start=R(t=c),color='magenta') #same direction as A, but with magnitude 3
accel_label=1.1*(3.0*A(c)/norm(A(c))+R(c)) #a location near the endpoint of this magnitude-3 vector to label the actual magnitude
P += text3d(str(float(floor(100*norm(A(t=c)))/100)),accel_label,color='magenta') #this messy code eliminates all but a few decimals
#plot the osculating plane and choose color and size and opacity
P += plotPlane(B(t=c),R(t=c),size,'grey',0.2)
#plot the point as well. Note the lower-case r, since it is a point, not a vector. we just want a dot, not an arrow.
P += point(r(t=c), thickness=30, color= 'green')
#Now compute some stuff about how fast the object is traveling
# recall that the speed is the magnitude of velocity vector
speed = V(t=c).norm()
accel = A(t=c).norm()
k = K(t=c)
#add axes, sized to be a little taller than the helix reaches
P += plot(axes(1.1*t1))
show(P,spin=1)
#Now create plots of K, A, V with color-coded legend
pK=plot(K,(t,t0,t1),legend_label='curvature',legend_color = 'green', color='green')
pA=plot(norm(A(t)),(t,t0,t1),legend_label='acceleration',legend_color = 'magenta', color='magenta')
pV=plot(norm(V),(t,t0,t1),legend_label='velocity',legend_color = 'blue', color='blue')
#print values at our t=c point.
print('At t = %f:'%c)
print('speed = %f'%speed)
print('scalar acceleration = %f'%accel)
print('curvature = %f' %k)
#finally, display the K-V-A plot
show(pA+pK+pV)
︡6cce92e8-b90c-450f-95e7-a92b76221eea︡{"file":{"filename":"7434aac8-71f4-464a-b416-80ee8f0218d5.sage3d","uuid":"7434aac8-71f4-464a-b416-80ee8f0218d5"}}︡{"stdout":"At t = 0.942478:\n"}︡{"stdout":"speed = 4.427318\n"}︡{"stdout":"scalar acceleration = 8.303751\n"}︡{"stdout":"curvature = 0.423226\n"}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/13500/tmp_unQoHR.svg","show":true,"text":null,"uuid":"6c804580-907d-426c-b8ae-38f9ffbf3072"},"once":false}︡{"done":true}︡
︠7cdc3b3a-f332-4b5d-b04e-86b801856432︠
︡6667723b-7491-46d4-bf34-0957e11dcdcd︡
︠6b017896-9856-40a4-81ec-185cf49bb75f︠
#########################################################
#
# interactive TANGENT PLANE PLOTTER. Shows a z=f(x,y) surface
# and lets you choose x=a, y=b and draws tangent vectors and tangent plane at (a,b,f(a,b))
# and as a bonus, shows you the partial derivatives and their values.
#
###########################################################
%var t,x,y,z
def axes(size,xcolor='red',ycolor='green',zcolor='blue',asp=(1,1)):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= asp[0]*size; zSize=asp[1]*size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
#return xAxis+yAxis+zAxis
def plotPlane(N,P,s=1,clr='grey',op=0.2,meshvalue=0):
x,y,z=var('x y z')
#N=normal vec, P = point (vector), s=size
#This returns a plot of plane with the given normal vector and point, with x,y,z going plus or minus s units from the point.
#optional arguments for color and opacity (default is nearly transparent grey)
return implicit_plot3d(N.dot_product(vector((x,y,z)))-N.dot_product(P),(x,P[0]-s,P[0]+s),(y,P[1]-s,P[1]+s),(z,P[2]-s,P[2]+s),color=clr,opacity=op,mesh=meshvalue)
#global values for interactive function defined below
xvals = (x,-2,2); yvals = (y,-2,2); ptvals=(1,n(1.5,12))
@interact(layout=dict(top=[['f','xrange','yrange'],['pnt','showTP','showDerivs','showAxes']]))
def _(f=input_box(sqrt(9-x^2-y^2),width=40),xrange=input_box(xvals, width=20),yrange=input_box(yvals, width = 20),pnt=input_box(ptvals,label="(a,b)",width=20),showTP=checkbox(false,label="Show tangent plane"),showDerivs=checkbox(false,label="Show derivatives"), showAxes=checkbox(true,label="Show axes")):
(a,b) = (pnt[0],pnt[1]) # we will draw tangent plane on surface at (a,b,f(a,b))
#(a,b) = (n(a,12),n(b,12))
#set options for surface plot etc
pointColor ='black'
surfaceColor = 'grey'; meshvalue = 0; opvalue= 0.4
#make the surface plot
P = plot3d(f,xrange, yrange, mesh = meshvalue, color = surfaceColor, opacity = opvalue)
# include axes if checked
if showAxes:
P += axes(1.2*max(xrange[2],yrange[2]),xcolor='grey', ycolor='grey',zcolor='grey')
#compute partial derivatives
fx(x,y) = f.diff(x)
fy(x,y) = f.diff(y)
# define point and its vector form
pt = (a,b,f(a,b))
Pt = vector(pt)
P += plot(point(pt,thickness = 10, color=pointColor))
#compute normal to tangent plane
N = vector([fx(a,b),fy(a,b),-1])
#add in the new plane if option checked
if showTP:
P += plotPlane(N,Pt, s=0.75, clr='yellow', op=0.4)
# draw normal vector
#P += plot(N/norm(N), start=Pt, color='black')
#include traces for y=b, x=a on xy-plane and on surface
#the xtrace is parallel to the x-axis, ytrace is parallel to y
xtraceColor='magenta'; ytraceColor='green'
P += parametric_plot3d((a,t,0),(t,yrange[1],yrange[2]),color=ytraceColor)
P += parametric_plot3d((t,b,0),(t,xrange[1],xrange[2]),color=xtraceColor )
P += parametric_plot3d((a,t,f(a,t)),(t,yrange[1],yrange[2]),color=ytraceColor)
P += parametric_plot3d((t,b,f(t,b)),(t,xrange[1],xrange[2]),color=xtraceColor )
#next, create the two unit tangent vectors
#for the ytrace, we take derivative of the parametric curve (a,t,f(a,t))
ytraceTangent = vector((0,1,fy(a,b)))
# for xtrace, differentiate (t, b, f(t,b))
xtraceTangent = vector((1, 0, fx(a,b)))
unitYTraceTan=ytraceTangent/norm(ytraceTangent)
unitXTraceTan=xtraceTangent/norm(xtraceTangent)
P += plot(unitYTraceTan, start = Pt, color='darkgreen')
P += plot(unitXTraceTan, start = Pt, color='darkred')
if showDerivs:
show('$f(x,y) ='+ latex(f)+ ',\quad\quad{\partial f\over\partial x} =' + latex(fx(x,y))+',\quad\quad {\partial f\over\partial y}=' + latex(fy(x,y))+'$')
show('At $(a,b)='+latex((a,b))+', \partial f/\partial x ' +'\doteq' + latex(n(fx(a,b),12))+',\quad \partial f/\partial y ' +'\doteq' + latex(n(fy(a,b),12))+'$')
#print('At (a,b) = (%s,%s), f_x=%f'%(a,b,fx(a,b)))
show(P)
︡9cb7615d-3c1e-4dae-ba34-329a4cc70319︡{"interact":{"controls":[{"control_type":"input-box","default":"sqrt(-x^2 - y^2 + 9)","label":"f","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"f","width":40},{"control_type":"input-box","default":"(x, -2, 2)","label":"xrange","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"xrange","width":20},{"control_type":"input-box","default":"(y, -2, 2)","label":"yrange","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"yrange","width":20},{"control_type":"input-box","default":"(1, 1.50)","label":"(a,b)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"pnt","width":20},{"control_type":"checkbox","default":false,"label":"Show tangent plane","readonly":false,"var":"showTP"},{"control_type":"checkbox","default":false,"label":"Show derivatives","readonly":false,"var":"showDerivs"},{"control_type":"checkbox","default":true,"label":"Show axes","readonly":false,"var":"showAxes"}],"flicker":false,"id":"1e495881-e1fe-41a2-9626-464647648748","layout":[[["f",4,null],["xrange",4,null],["yrange",4,null]],[["pnt",3,null],["showTP",3,null],["showDerivs",3,null],["showAxes",3,null]],[["",12,null]]],"style":"None"}}︡{"done":true}︡
︠7779dc57-31f7-44a4-8c59-b30a0f5fc70f︠
################################
#
# surface plot compared with contour plot
#
##################################
%var x,y,z,t
import matplotlib as mp
# pick a function and range to plot it
#f=y^2-x^2
f=6*x*y*exp(-(x^2+y^2)/2)
#f= x^3 - 3*x*y^2 #monkey saddle!!
s=3
xrange=(x,-s,s); yrange=(y,-s,s)
#set a colormap for the contour and the matplotlib plot
cmap ='coolwarm'
#cmap = mp.colors.ListedColormap([0,0,0])
#the first one is a standard 3d plot that we can rotate,etc.
plot3d(f, xrange, yrange,aspect_ratio=[1.,1.,1.], mesh=1,color='lightblue')
# next plot is fixed, but fancy and prettier
plot3d_using_matplotlib(f, xrange, yrange,azim=300,cmap=cmap)
#next plot is contour plot which matches second plot in color map
cp =contour_plot(f, xrange, yrange,colorbar=True,contours=20,labels=False,fill=true,figsize=6,linestyles='-',cmap=cmap)
#let's add in the gradient vectors
gv = plot_vector_field(f.gradient(),xrange, yrange)
show(cp+gv)
#here is another contour style that you may like: black lines with labels
import matplotlib as mp
blackCmap=mp.colors.ListedColormap([0,0,0])
contour_plot(f, xrange, yrange,colorbar=False,contours=20,labels=True,fill=false,figsize=6,linestyles='-',cmap=blackCmap)
︡292ceea9-83e9-4de1-a5b3-6e5b04610b66︡{"file":{"filename":"c7366f16-dc3f-4bfa-a357-3515d03642c6.sage3d","uuid":"c7366f16-dc3f-4bfa-a357-3515d03642c6"}}︡{"file":{"filename":"9eff2f21-b538-400b-88af-834a1f9f79d3.svg","show":true,"text":null,"uuid":"c2c9c68d-fec0-4438-8d00-c00bdf926b1f"},"once":false}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/120/tmp_PgFB3C.svg","show":true,"text":null,"uuid":"7de70e33-d96f-490b-88eb-d656049ceca5"},"once":false}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/120/tmp_rQSOV7.svg","show":true,"text":null,"uuid":"c94f2612-4eda-4208-b73d-026ce2970332"},"once":false}︡{"done":true}︡
︠6e763e56-a770-4950-ae11-7bdc3b1e9ec5︠
# chain rule practice 1
#DONT FORGET TO DECLARE VARIABLES!!
%var x,y,z,t
#First define F(x,y,z)
F(x,y,z)= x+y^2+y*z
#next, express x,y,z as functions of another variable t
#in other words, G(t) := F(f(t), g(t), h(t))
f=t^2; g = cos(t); h=3*exp(t^2)
G = F(x=f, y=g, z = h)
show('$F(x,y,z)='+latex(F(x,y,z))+'$')
show('$G(t)=F('+latex(f)+','+latex(g)+','+latex(h)+')='+latex(G(t=t))+'$')
show('$\partial G/\partial t='+latex(G.diff(t))+'$')
︡ca92e545-275c-4ebb-94f9-952ecdce8798︡{"html":"$F(x,y,z)= y^{2} + y z + x $
"}︡{"html":"$G(t)=F( t^{2} , \\cos\\left(t\\right) , 3 \\, e^{\\left(t^{2}\\right)} )= t^{2} + \\cos\\left(t\\right)^{2} + 3 \\, \\cos\\left(t\\right) e^{\\left(t^{2}\\right)} $
"}︡{"html":"$\\partial G/\\partial t= 6 \\, t \\cos\\left(t\\right) e^{\\left(t^{2}\\right)} - 2 \\, \\cos\\left(t\\right) \\sin\\left(t\\right) - 3 \\, e^{\\left(t^{2}\\right)} \\sin\\left(t\\right) + 2 \\, t $
"}︡{"done":true}︡
︠366d71d5-185f-4c5d-b649-253f29acf081︠
# chain rule practice 2
#DONT FORGET TO DECLARE VARIABLES!!
%var x,y,z,t,s
#First define F(x,y,z)
F(x,y,z)= x+y^2+y*z
#next, express x,y,z as functions of two other variables t and s
#in other words, G(t) := F(f(t,s), g(t,s), h(t,s))
f=t^2+s; g = cos(t*s); h=3*exp(t^2-s^3)
G = F(x=f, y=g, z = h)
show('$F(x,y,z)='+latex(F(x,y,z))+'$')
show('$G(t,s)=F('+latex(f)+','+latex(g)+','+latex(h)+')='+latex(G(t=t,s=s))+'$')
show('$\partial G/\partial t='+latex(G.diff(t))+'$')
show('$\partial G/\partial s='+latex(G.diff(s))+'$')
︡2d9f03ae-f2c8-41b3-b4ca-19e009dcf9ed︡{"html":"$F(x,y,z)= y^{2} + y z + x $
"}︡{"html":"$G(t,s)=F( t^{2} + s , \\cos\\left(s t\\right) , 3 \\, e^{\\left(-s^{3} + t^{2}\\right)} )= t^{2} + \\cos\\left(s t\\right)^{2} + 3 \\, \\cos\\left(s t\\right) e^{\\left(-s^{3} + t^{2}\\right)} + s $
"}︡{"html":"$\\partial G/\\partial t= 6 \\, t \\cos\\left(s t\\right) e^{\\left(-s^{3} + t^{2}\\right)} - 2 \\, s \\cos\\left(s t\\right) \\sin\\left(s t\\right) - 3 \\, s e^{\\left(-s^{3} + t^{2}\\right)} \\sin\\left(s t\\right) + 2 \\, t $
"}︡{"html":"$\\partial G/\\partial s= -9 \\, s^{2} \\cos\\left(s t\\right) e^{\\left(-s^{3} + t^{2}\\right)} - 2 \\, t \\cos\\left(s t\\right) \\sin\\left(s t\\right) - 3 \\, t e^{\\left(-s^{3} + t^{2}\\right)} \\sin\\left(s t\\right) + 1 $
"}︡{"done":true}︡
︠f740fb18-8841-4041-aa14-0f01765e2bf0︠
#MIDPOINT APPROXIMATION OF SIMPLE DOUBLE INTEGRALS
%var x,y
# the code for drawing boxes needs to be imported
from sage.plot.plot3d.shapes import Box
#define the function
f=(sin(x+y))^2
#f=13-x^2-y^2
#the function will be integrated over a rectangle centered at (cx, cy) with size sx, sy
sx,sy=2,3
cx,cy=0,0
#how much to subdivide in the x direction (m) and y direction (n)
m=30
n= 90
boxes=[]
dx=2*sx/m
dy=2*sy/n
vol=0
for i in range(m):
xi=(i+0.5)*dx-sx
for j in range(n):
yj=(j+0.5)*dy-sy
h=f(x=xi, y=yj)
boxes.append(Box([dx/2, dy/2,abs(h)/2],color=(1+(-1)^(i+j),-1-(-1)^(i+j),.5),opacity=.8,specular=0,shininess=0).translate(xi,yj,h/2))
vol += h*dx*dy
#Display the boxes
show(sum(boxes),spin=1)
#This is the code to actually do the symbolic integral
ii=integrate(integrate(f,(x,cx-sx, cx+sx)),(y,cy-sy, cy+sy))
#Display the integral value and the approximation value
show('$\int_{'+latex(cy-sy)+'}^{'+latex(cy+sy)+'}\int_{'+latex(cx-sx)+'}^{'+latex(cx+sx)+'}(' +latex(f)+')dx dy='+latex(ii)+'='+latex(ii.n())+'$')
show('Midpoint approximation '+'$ = '+latex(vol)+'$')
︡baa0a2c6-be17-4c1d-ad4d-a4f3c052dbf7︡{"file":{"filename":"1f54dbb8-c1f1-408e-a968-38c3e70090bb.sage3d","uuid":"1f54dbb8-c1f1-408e-a968-38c3e70090bb"}}︡{"html":"$\\int_{ -3 }^{ 3 }\\int_{ -2 }^{ 2 }( \\sin\\left(x + y\\right)^{2} )dx dy= \\frac{1}{4} \\, \\cos\\left(10\\right) - \\frac{1}{4} \\, \\cos\\left(2\\right) + 12 = 11.8942688268677 $
"}︡{"html":"Midpoint approximation $ = 11.8938763056596 $
"}︡{"done":true}︡
︠1d6016ab-9168-4716-b742-456d60a7146a︠
#problem 15.6.33
%var x,y,z
#reuse axes function from earlier
def axes(size,xcolor='red',ycolor='green',zcolor='blue',asp=(1,1)):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= asp[0]*size; zSize=asp[1]*size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
#return xAxis+yAxis+zAxis
###############################################################################
def plotPlane(N,P,s,clr='grey',op=0.2,meshvalue=0):
x,y,z=var('x y z')
#N=normal vec, P = point (vector), s=size
#This returns a plot of plane with the given normal vector and point, with x,y,z going plus or minus s units from the point.
#optional arguments for color and opacity (default is nearly transparent grey)
return implicit_plot3d(N.dot_product(vector((x,y,z)))-N.dot_product(P),(x,P[0]-s,P[0]+s),(y,P[1]-s,P[1]+s),(z,P[2]-s,P[2]+s),color=clr,opacity=op,mesh=meshvalue)
#default values
color1='red';color2='lightgreen';color3='lightblue'
eqn1=1-y-z
eqn2=y-sqrt(x)
xlims = (x,0,1);ylims=(y,0,1);zlims=(z,0,1)
P1=Graphics();P2=Graphics(); P3=Graphics()
P1+= implicit_plot3d(eqn1,xlims, ylims, zlims,mesh=1,color=color1)
P1+= implicit_plot3d(eqn2,xlims, ylims, zlims,mesh=1,color=color2)
#axes set up. need size, but don't use limits if not shown
ax = axes(1.25)
#make some translucent gray planes to bound the figure
yzPlane=plotPlane(vector((1,0,0)),vector((0,.5,.5)),.5)
xyPlane=plotPlane(vector((0,0,1)),vector((.5,.5,0)),.5)
print("Problem 15.6.33")
show(P1+P2+P3+ax+yzPlane+xyPlane, spin = 0)
︡c57f7925-3ab3-4d1e-a7bc-a9f408b9b19f︡{"stdout":"Problem 15.6.33\n"}︡{"file":{"filename":"06d682c2-aa3d-49a8-b01f-4c3a1e20f885.sage3d","uuid":"06d682c2-aa3d-49a8-b01f-4c3a1e20f885"}}︡{"done":true}︡
︠5a78a40d-023a-4e69-918c-946f73e8b4e3s︠
# 15.6.34
# we will compute some integrals using the function f(x,y,z)=1 and hence find the volume.
%var x,y,z
#reuse axes function from earlier
def axes(size,xcolor='red',ycolor='green',zcolor='blue',asp=(1,1)):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= asp[0]*size; zSize=asp[1]*size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
#return xAxis+yAxis+zAxis
###############################################################################
def plotPlane(N,P,s,clr='grey',op=0.2,meshvalue=0):
x,y,z=var('x y z')
#N=normal vec, P = point (vector), s=size
#This returns a plot of plane with the given normal vector and point, with x,y,z going plus or minus s units from the point.
#optional arguments for color and opacity (default is nearly transparent grey)
return implicit_plot3d(N.dot_product(vector((x,y,z)))-N.dot_product(P),(x,P[0]-s,P[0]+s),(y,P[1]-s,P[1]+s),(z,P[2]-s,P[2]+s),color=clr,opacity=op,mesh=meshvalue)
#default values
color1='red';color2='lightgreen';color3='lightblue'
eqn1=y-1+x
eqn2=z+x^2-1
xlims = (x,0,1);ylims=(y,0,1);zlims=(z,0,1)
P1=Graphics();P2=Graphics(); P3=Graphics()
P1+= implicit_plot3d(eqn1,xlims, ylims, zlims,mesh=1,color=color1)
P1+= implicit_plot3d(eqn2,xlims, ylims, zlims,mesh=1,color=color2)
#axes set up. need size, but don't use limits if not shown
ax = axes(1.25)
#make some translucent gray planes to bound the figure
yzPlane=plotPlane(vector((1,0,0)),vector((0,.5,.5)),.5)
xyPlane=plotPlane(vector((0,0,1)),vector((.5,.5,0)),.5)
#also create planes parallel to the three coordinate planes to visualize different ways to set up integrals
#we make these planes yellow and meshed with high opacity so they are very visible
xyPlane2= plotPlane(vector((0,0,1)),vector((.5,.5,.5)),.5,clr='yellow',meshvalue=1, op=.7)
yzPlane2= plotPlane(vector((1,0,0)),vector((.5,.5,.5)),.5,clr='yellow',meshvalue=1, op=.7)
xzPlane2= plotPlane(vector((0,1,0)),vector((.5,.5,.5)),.5,clr='yellow',meshvalue=1, op=.7)
print("Problem 15.6.34")
# this is the "dy dz dx" version in the problem statement
Iyzx =integrate(integrate(integrate(1,y,0,1-x),z,0,1-x^2),x,0,1)
# this is the "dy dx dz" version in the problem statement
Iyxz =integrate(integrate(integrate(1,y,0,1-x),x,0,sqrt(1-z)),z,0,1)
# this is the "dx dy dz" version
Ixyz=integrate(integrate(integrate(1,x,0,1-y),y,1-sqrt(1-z),1),z,0,1)+integrate(integrate(integrate(1,x,0,sqrt(1-z)),y,0,1-sqrt(1-z)),z,0,1)
#they should be the same value!
print("Iyzx = %s"%Iyzx)
print("Ixyz = %s"%Ixyz)
print("Iyxz = %s"%Iyxz)
show(P1+P2+P3+ax+xyPlane+yzPlane+xyPlane2, spin = 0)
︡e7928069-bd83-4464-9694-712f317487a9︡{"stdout":"Problem 15.6.34\n"}︡{"stdout":"Iyzx = 5/12\n"}︡{"stdout":"Ixyz = 5/12\n"}︡{"stdout":"Iyxz = 5/12\n"}︡{"file":{"filename":"412dbdb5-de52-4e4f-8954-66e690fb49bd.sage3d","uuid":"412dbdb5-de52-4e4f-8954-66e690fb49bd"}}︡{"done":true}
︠9291405e-2500-4fa7-bbdf-caeb045a173fs︠
###########################
# These three cells illustrate ways to see how the element of volume depends on dRho, dTheta, and dPhi
#
# the first two use the spherical_plot3d function which automatically draws a surface of the form rho=f(theta, phi).
# the third cell just uses the regular plot3d function, but explicitly specifies the transformation that takes the spherical variables and converts them
# into cartesian coordinates (i.e., takes (rho, theta, phi)--->(x,y,z)). This method is very flexible but also pretty sophisticated.
#
#############################################################################
# this first picture reminds you that longitudes get close to one another at the poles!
%var r,theta, phi
f(theta,phi)=1
s=0.1
thetaRange=(theta,0,2*pi)
phiRange=(phi,0,pi/2)
thetaRange1=(theta,0,pi*s)
phiRange1=(phi,0,pi/2)
S0 = spherical_plot3d(f,thetaRange, phiRange,mesh=0,opacity=0.2, color='yellow')
S1=spherical_plot3d(f,thetaRange1, phiRange1,mesh=0, color='lightblue',opacity=0.9)
S2 =spherical_plot3d(f*1.1,thetaRange1, phiRange1,color='red',mes0=1,opacity=0.5)
show(S0+S1+S2)
︡84079866-2cd8-430e-abce-898a76d2e8e8︡{"file":{"filename":"bce5fba2-3a91-4ea7-91e8-5fe43903dab6.sage3d","uuid":"bce5fba2-3a91-4ea7-91e8-5fe43903dab6"}}︡{"done":true}︡
︠be4a4a93-9951-4b6f-a71d-55b138f56880s︠
# In this second picture, try altering the values for theta and dtheta, etc. to see how the element of volume changes
%var theta, phi
########## try altering these values
th =0 #starting value of theta
deltaTh =0.1*pi
ph=pi/4 #starting phi value
deltaPh = 0.1*pi
r=1 #starting rho
deltaR=0.1
##############################
def onSphere(rho,theta, phi):
#coords on a sphere
return((rho*sin(phi)*cos(theta),rho*sin(phi)*sin(theta), rho*cos(phi)))
f(theta,phi)=r #equation of sphere rho=1
thetaRange=(theta,th,th+ deltaTh)
phiRange=(phi,ph,ph+ deltaPh)
#for a hemisphere
thetaRange0=(theta,0,pi*2)
phiRange0=(phi,0,pi/2)
L1= line([(0,0,1),(0,0,0),(cos(th),sin(th),0),(0,0,0),(cos(th+deltaTh),sin(th+deltaTh),0)])
L2=line([(0,0,0),onSphere(r,th,ph),(0,0,0),onSphere(r, th, ph+deltaPh),(0,0,0),onSphere(r,th+deltaTh,ph),(0,0,0),onSphere(r, th+deltaTh, ph+deltaPh)])
L3 = line([(sin(ph)*cos(th),sin(ph)*sin(th),0),onSphere(r,th,ph)],color='red',linestyle=':')
L4 = line([(sin(ph)*cos(th+deltaTh),sin(ph)*sin(th+deltaTh),0),onSphere(r,th+deltaTh,ph)],color='red',linestyle=':')
Arc1= parametric_plot3d((sin(ph)*cos(theta),sin(ph)*sin(theta),0),thetaRange)
Arc2= parametric_plot3d((sin(ph+ deltaPh)*cos(theta),sin(ph+ deltaPh)*sin(theta),0),thetaRange)
S0 = spherical_plot3d(r,thetaRange0, phiRange0,mesh=0,opacity=0.2, color='yellow') #yellow hemisphere
S1=spherical_plot3d(r,thetaRange, phiRange, color='lightblue',opacity=0.7) #inner surface
S2 =spherical_plot3d(r+deltaR,thetaRange, phiRange,color='red',mesh=0,opacity=0.7) #outer surface
show(S0+S1+S2+L1+L2+L3+L4+Arc1+Arc2,spin=5)
︡cf484d6c-12d7-4ef7-89b3-5a8ec25ed04c︡{"file":{"filename":"92bb6e60-f7e0-4630-9a81-deda26b6601a.sage3d","uuid":"92bb6e60-f7e0-4630-9a81-deda26b6601a"}}︡{"done":true}︡
︠ecc80e43-c8ba-4f24-8f7a-eba91ce52646︠
︡0260df3e-9e1f-43e6-aa9b-432476ae70ff︡
︠8669e62b-bba6-4e34-b7b2-9b662d41d9b5s︠
%var rho,theta, phi,x,y,z
############
#
# This final picture shows two different elements of volume, sitting on a "taco" of "theta wedges"
# It also illustrates very sophisticated ways to use the spherical-to-cartesian transformation to draw interesting shapes, where we
# declare which variables are independent (the undecleared one, that is not in the final brackets, is the dependent variable)
#
######################
##################################
#define three types of spherical trans, depending on the dependent and indep variable
#rho is dep, theta and phi are indep
SphereRho = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [theta,phi])
#phi is dep, theta and rho are indep
SpherePhi = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [theta,rho])
#theta is dep, phi and rho are indep
SphereTheta = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [phi,rho])
##################################
#specify starting theta and phi and delta theta and delta phi.
th =0 #starting value of theta
deltaTh =0.1*pi
ph=.1*pi #starting phi value
ph2=.4*pi #second starting phi value
deltaPh = 0.05*pi
#start rho and delta rho
rh=1
deltaRh=0.1
thetaRange=(theta,th,th+ deltaTh)
phiRange=(phi,ph,ph+ deltaPh)
phiRange2=(phi,ph2,ph2+ deltaPh)
rhoRange=(rho, rh,rh+deltaRh)
###################################
#yellow hemisphere, radius = 1
P0=implicit_plot3d(x^2+y^2+z^2-1,(x,-1,1),(y,-1,1),(z,0,1),color='yellow',opacity=0.1)
#"theta wedge"
P1=plot3d(th, (phi,0,pi/2),(rho,0,1),transformation=SphereTheta, color='grey', opacity =0.5)
P2=plot3d(th+deltaTh, (phi,0,pi/2),(rho,0,1),transformation=SphereTheta, color='brown', opacity =0.3)
#first "phi wedge"
P3=plot3d(rh, thetaRange,phiRange,transformation=SphereRho, color='green',opacity=0.3)
P6=plot3d(rh+deltaRh, thetaRange,phiRange,transformation=SphereRho, color='green',opacity=0.3)
P4=plot3d(th,phiRange,rhoRange, transformation=SphereTheta, color='red', opacity =0.3)
P5=plot3d(th+deltaTh,phiRange,rhoRange, transformation=SphereTheta, color='red', opacity =0.3)
PhiWedge=P3+P6+P4+P5
#second "phi wedge"
P7=plot3d(rh, thetaRange,phiRange2,transformation=SphereRho, color='green',opacity=0.3)
P8=plot3d(rh+deltaRh, thetaRange,phiRange2,transformation=SphereRho, color='green',opacity=0.3)
P9=plot3d(th,phiRange2,rhoRange, transformation=SphereTheta, color='red', opacity =0.3)
P10=plot3d(th+deltaTh,phiRange2,rhoRange, transformation=SphereTheta, color='red', opacity =0.3)
PhiWedge2=P7+P8+P9+P10
show(P1+P2+P0+PhiWedge+PhiWedge2,spin=0)
︡2a162ae4-dafb-4a13-9416-b4e530b38d44︡{"file":{"filename":"d189cc48-be26-4a1a-8ec2-a7ff7e2100be.sage3d","uuid":"d189cc48-be26-4a1a-8ec2-a7ff7e2100be"}}︡{"done":true}︡
︠ed38a467-2fb5-453e-8299-b41acbc1ee26s︠
#pictures from quiz 5
#Problem 1 (cylindrical)
%var r,theta,z
zRange=(z,0,2)
rRange=(r,0,1)
thetaRange=(theta,0,2*pi)
##IMPORTANT NOTE: cylindrical_plot3d allows one ONLY to plot a function of the form r=f(theta, z). So
#in this problem, z=2r, so we write r=z/2.
Cone=cylindrical_plot3d(z/2,thetaRange,zRange,color='blue', mesh=1)
Cyl =cylindrical_plot3d(1,thetaRange,zRange,color='red', mesh=0, opacity=0.4)
show(Cone+Cyl)
︡3b90a739-871a-4536-af21-ff98c57d0fb7︡{"file":{"filename":"69b24813-68a0-44a1-9f9a-8a50e9116208.sage3d","uuid":"69b24813-68a0-44a1-9f9a-8a50e9116208"}}︡{"done":true}︡
︠9d00a7cd-195b-44a9-a72e-487338eb3c0a︠
#quiz 5, problem 2 (spherical)
%var rho, theta, phi
##################################
#define three types of spherical trans, depending on the dependent and indep variable
#rho is dep, theta and phi are indep
SphereRho = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [theta,phi])
#phi is dep, theta and rho are indep
SpherePhi = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [theta,rho])
#theta is dep, phi and rho are indep
SphereTheta = (rho*cos(theta)*sin(phi), rho*sin(theta)*sin(phi), rho*cos(phi), [phi,rho])
##################################
thetaRange=(theta,0,1.9*pi) #really should be 0 to 2pi but this way we can peek inside!
rhoRange=(rho, 1,2)
phiRange=(phi,0,pi/6)
Cone= plot3d(pi/6, thetaRange, rhoRange,transformation=SpherePhi,color='blue', mesh=1,opacity=0.3)
S1=plot3d(1,thetaRange,phiRange, transformation=SphereRho, color='red',mesh=1)
S2=plot3d(2,thetaRange,phiRange, transformation=SphereRho, color='green',mesh=1)
#'bottoms' of objects
S1b=plot3d(1,thetaRange,(phi,pi/6,pi/2), transformation=SphereRho, color='red',opacity=0.1)
S2b=plot3d(2,thetaRange,(phi,pi/6,pi/2), transformation=SphereRho, color='green',opacity=0.1)
Coneb=plot3d(pi/6, thetaRange, (rho,0,1),transformation=SpherePhi,color='blue', opacity=0.1)
show(Cone+Coneb+S1+S2+S1b+S2b)
︡13af5d1c-ef40-4272-9ec2-574a8df16cc6︡{"file":{"filename":"11f2ce0c-59b2-47d9-93e8-79f2507198d6.sage3d","uuid":"11f2ce0c-59b2-47d9-93e8-79f2507198d6"}}︡{"done":true}︡
︠0f3101a9-51e0-4c43-9337-de3dbcb97f1d︠
#INTERACTIVE U,V--->X,Y TRANSFORMATION EXPLORER
#INPUT FORMULA FOR X=X(U,V), Y=Y(U,V), A STARTING POINT (U0,Y0) AND THE SIZE AND NUMBER OF DIVISIONS
#THE TWO GRAPHS ARE DISPLAYED, ALONG WITH THE JACOBIAN WHICH MEASURES THE AREA SCALING FACTOR.
%var u,v,s
xdef=u^2-v^2
ydef = 2*u*v
uvdef= (2,3)
@interact(layout=dict(top=[['x','y','D'],['n','pt','ss']]))
def _(x=input_box(xdef, label='x(u,v)'), y=input_box(ydef, label='y(u,v)'), D=input_box(1,label='delta',width=20), n=input_box(5,label='nDivisions',width=20),pt=input_box(uvdef,label='(u0,v0)',width=20),ss=checkbox(false,label='same scale')):
u0=pt[0]; v0=pt[1]
d=D/n
UVPlot=Graphics()
XYPlot=Graphics()
#add dots to each plot
UVPlot += point((u0+d/2,v0+d/2))
XYPlot += point((x(u=u0+d/2,v=v0+d/2),y(u=u0+d/2,v=v0+d/2)))
for t in srange(u0-D,u0+D+d,d):
UVPlot += line([(t, v0-D),(t,v0+D)])
XYPlot += parametric_plot((x(u=t,v=s),y(u=t,v=s)),(s,pt[1]-D, pt[1]+D))
UVPlot += line([(u0,v0-D),(u0,v0+D)],color='red',thickness=3)
UVPlot += line([(u0-D,v0),(u0+D,v0)],color='green',thickness=3)
for t in srange(v0-D,v0+D+d,d):
UVPlot += line([(u0-D,t),(u0+D,t)])
XYPlot += parametric_plot((x(u=s,v=t),y(u=s,v=t)),(s,pt[0]-D, pt[0]+D))
XYPlot += parametric_plot((x(u=u0,v=s),y(u=u0,v=s)),(s,pt[1]-D, pt[1]+D),color='red',thickness=3)
XYPlot += parametric_plot((x(u=s,v=v0),y(u=s,v=v0)),(s,pt[0]-D, pt[0]+D),color='green',thickness=3)
L = line([(-3,2), (3,2)],thickness=0)
L.axes(false)
if ss:#if same scale give uniform xmins, etc
xmn=min(XYPlot.xmin(),UVPlot.xmin())
XYPlot.xmin(xmn)
UVPlot.xmin(xmn)
xmx=max(XYPlot.xmax(),UVPlot.xmax())
XYPlot.xmax(xmx)
UVPlot.xmax(xmx)
graphics_array([UVPlot,L,XYPlot]).show(aspect_ratio=1)
jacobian=x.diff(u)*y.diff(v)-x.diff(v)*y.diff(u)
print('UV is left graph, XY is right graph. Jacobian scale factor is %f' %abs(jacobian(u=u0,v=v0)))
show(' $\partial x/\partial u='+latex(x.diff(u))+',\quad\partial x/\partial v ='+latex(x.diff(v))+'$')
show(' $\partial y/\partial u='+latex(y.diff(u))+',\quad\partial y/\partial v ='+latex(y.diff(v))+'$')
show(' $\partial(x,y)/\partial (u,v)='+latex(jacobian)+'$')
︡b4bd7d65-525b-4249-8cac-83fd03fafb8a︡{"interact":{"controls":[{"control_type":"input-box","default":"u^2 - v^2","label":"x(u,v)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"x","width":null},{"control_type":"input-box","default":"2*u*v","label":"y(u,v)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"y","width":null},{"control_type":"input-box","default":1,"label":"delta","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"D","width":20},{"control_type":"input-box","default":5,"label":"nDivisions","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"n","width":20},{"control_type":"input-box","default":"(2, 3)","label":"(u0,v0)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"pt","width":20},{"control_type":"checkbox","default":false,"label":"same scale","readonly":false,"var":"ss"}],"flicker":false,"id":"53067e16-cd0d-4d6e-bb11-bcf6b8a0cbe3","layout":[[["x",4,null],["y",4,null],["D",4,null]],[["n",4,null],["pt",4,null],["ss",4,null]],[["",12,null]]],"style":"None"}}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/149/tmp_AYR79_.svg","show":true,"text":null,"uuid":"3e226350-55b3-4058-af44-17ea6c3c7023"},"once":false}︡{"html":"gradient field =
"}︡{"html":"$\\displaystyle \\left(2 \\, x,\\,6 \\, y\\right)$
"}︡{"file":{"filename":"ce07b6c2-6cce-4460-83a5-867a5e5a4b2d.sage3d","uuid":"ce07b6c2-6cce-4460-83a5-867a5e5a4b2d"}}︡{"html":"vector field =
"}︡{"html":"($\\displaystyle x$, $\\displaystyle y$, $\\displaystyle 0$)
"}︡{"done":true}︡
︠109801e8-98dd-4aab-b62e-cec4ea52ef3fs︠
# PLOTTING VECTOR FIELDS
#
%var x,y,z,u,v,t
f = x^2+3*y^2 #ellipse
g = f.gradient() #duh
a = 3 # scale
cp = contour_plot(f,(x,-a,a),(y,-a,a),fill=false)
vf = plot_vector_field(g,(x,-a,a),(y,-a,a),color='blue')
show(cp+vf)
show("gradient field =");
show(g)
︡0203b30f-3f8b-45a6-9d49-c4c4a45a121f︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/149/tmp_i8Itib.svg","show":true,"text":null,"uuid":"cc1d2175-dba2-4bf9-807a-4b32fc54fd58"},"once":false}︡{"html":"gradient field =
"}︡{"html":"$\\displaystyle \\left(2 \\, x,\\,6 \\, y\\right)$
"}︡{"done":true}︡
︠01b00799-7b05-4cbd-bdeb-1747e45acf8ds︠
# a 3-d example
%var x,y,z
P=x+y
Q=y
R=z*y
a = 3 # scale
xrange=(x,-a,a)
yrange=(y,-a,a)
zrange=(z,-a,a)
plot_vector_field3d((P,Q,R),xrange,yrange,zrange)
show("vector field =");
show((P,Q,R))
︡333c11b4-f6c3-400c-b0d7-d2e872eaa84c︡{"file":{"filename":"645e93db-4a0c-4c65-ae3a-a988a5b5bb06.sage3d","uuid":"645e93db-4a0c-4c65-ae3a-a988a5b5bb06"}}︡{"html":"vector field =
"}︡{"html":"($\\displaystyle x + y$, $\\displaystyle y$, $\\displaystyle y z$)
"}︡{"done":true}
︠4abff4b1-9f97-4a43-8b24-0a98b79fc8e1︠
%var x,y
P=1+sin(y); Q=1+x*cos(y)
f=y+x+x*sin(y)
show(P.diff(y))
show(Q.diff(x))
vp=plot_vector_field((P,Q), (x,0,3), (y, 0,3),aspect_ratio=1)
cp = contour_plot(f,(x,0,3), (y, 0,3),fill=false)
show(vp+cp)
︡0cacfe45-8ae7-4d08-8ea0-2840efdbd71e︡{"html":"$\\displaystyle \\cos\\left(y\\right)$
"}︡{"html":"$\\displaystyle \\cos\\left(y\\right)$
"}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/149/tmp_m4kzNT.svg","show":true,"text":null,"uuid":"9e06d680-33b5-4f39-be5e-9270b59d4bda"},"once":false}︡{"done":true}︡
︠d2af7142-7348-4cf4-8115-83c2229738bas︠
%var x,y,t
P=2*x*y; Q=x^2
f=y+x+x*sin(y)
show(P.diff(y))
show(Q.diff(x))
vp=plot_vector_field((P,Q), (x,0,3.5), (y, 0,3.5),aspect_ratio=1)
cp = contour_plot(f,(x,0,3), (y, 0,3),fill=false)
circ1=parametric_plot((cos(t)+2,sin(t)+2),(t,0,pi),thickness=3)
circ2=parametric_plot((cos(t)+2,1.5*sin(t)+2),(t,pi,2*pi),thickness=3)
wiggle=parametric_plot((2*t/(2*pi)+1,2+.1*sin(6*t)),(t,0,2*pi),thickness=3)
pts=point([(1,2),(3,2)],size=40)
ta=text("A",(2,3.1),fontsize=16)
tb=text("B",(2,1.6),fontsize=16)
tc=text("C",(1,1),fontsize=16)
show(vp+circ1+wiggle+pts+circ2+ta+tb+tc)
︡043bacd9-359a-48ba-a9ba-f8e6a15f5562︡{"html":"$\\displaystyle 2 \\, x$
"}︡{"html":"$\\displaystyle 2 \\, x$
"}︡{"file":{"filename":"/home/user/.sage/temp/project-2ad65871-00da-4e99-b1a6-6cbf035f8bc5/149/tmp_Moqrth.svg","show":true,"text":null,"uuid":"9a9b327c-8bd1-4adf-b065-651e6235009d"},"once":false}︡{"done":true}︡
︠4cad4bc1-968b-4a22-bd5f-6451cde583cb︠
#Examples of parametric equations of surfaces
%var u,v
r=1;R=4; h=.5
nu=10; nv =10
u0=0;u1=4*pi
v0=.5; v1=1
urange=(u,u0,u1); vrange=(v,v0,v1)
r(u,v) = (R*v*cos(u),R*v*sin(u),u)
parametric_plot3d(r(u,v),urange,vrange,mesh=0,opacity=0.2,plotpoints=10)
︡ca7faac9-961d-4923-855a-3562eb86aca4︡{"file":{"filename":"8024ceb7-6870-43f9-bffe-bf52e5198f25.sage3d","uuid":"8024ceb7-6870-43f9-bffe-bf52e5198f25"}}︡{"done":true}︡
︠710bd2dc-e9e9-4796-95d5-f5e722e1a78bs︠
#INTERACTIVE Parametric Surface EXPLORER
#INPUT FORMULA FOR X=X(U,V), Y=Y(U,V), Z(U,V); also input a "central" point (U0,V0) AND THE number of gridlines
#TWO GRAPHS ARE DISPLAYED: the rectangular UV "domain space" ALONG WITH THE 3D surface in xyz space.
%var u,v,s
#defaults : change them if you wish!
xdef=cos(u)*v
ydef = sin(u)*v
zdef = u
urangedef=(u,0,2*pi)
vrangedef=(v,1,2)
uvdef= (pi/4,1.5)
#define a function to draw axes. Very useful for understanding plots
def axes(size,xcolor='red',ycolor='green',zcolor='blue'):
#returns a plot of axis vectors, labeled, with size
xSize=size; ySize= size; zSize=size
xAxis = plot(vector((xSize,0,0)),color=xcolor,width=3)
yAxis = plot(vector((0,ySize,0)),color=ycolor,width=3)
zAxis = plot(vector((0,0,zSize)),color=zcolor,width=3)
#Labels for axes
#set fontsize (change if you need to)
fsize=20
xT=text3d("X",(xSize*1.1,0,0),color=xcolor,fontsize=fsize)
yT=text3d("Y",(0,ySize*1.1,0),color=ycolor,fontsize=fsize)
zT=text3d("Z",(0,0,zSize*1.1),color=zcolor,fontsize=fsize)
return xAxis+xT+yAxis+yT+zAxis+zT
@interact(layout=dict(top=[['x','y','z','n'],['pt','urange','vrange','showax']]))
def _(x=input_box(xdef, label='x(u,v)'), y=input_box(ydef, label='y(u,v)'),z=input_box(zdef, label='z(u,v)'), n=input_box(6,label='nDivisions',width=20),pt=input_box(uvdef,label='(u0,v0)',width=40),urange=input_box(urangedef,label='u range',width=36),vrange=input_box(vrangedef,label='v range',width=36),showax=checkbox(false,label='show axes')):
u0=pt[0]; v0=pt[1]
Du=urange[2]-urange[1]
du=Du/n
Dv=vrange[2]-vrange[1]
dv=Dv/n
UVPlot=Graphics()
SPlot= Graphics()
#add dots to each plot
UVPlot += point((u0+du/2,v0+dv/2),color='black')
SPlot += point((x(u=u0+du/2,v=v0+dv/2),y(u=u0+du/2,v=v0+dv/2),z(u=u0+du/2,v=v0+dv/2)),color='black',thickness=5)
for t in srange(urange[1],urange[2]+du,du):
UVPlot += line([(t, vrange[1]),(t,vrange[2])],color='red')
SPlot += parametric_plot((x(u=t,v=s),y(u=t,v=s),z(u=t,v=s)),(s,vrange[1], vrange[2]),color='red')
UVPlot += line([(u0,vrange[1]),(u0,vrange[2])],color='red',thickness=3)
SPlot += parametric_plot((x(u=u0,v=s),y(u=u0,v=s),z(u=u0,v=s)),(s,vrange[1], vrange[2]),color='red',thickness=5)
UVPlot += line([(urange[1],v0),(urange[2],v0)],color='green',thickness=3)
SPlot += parametric_plot((x(u=s,v=v0),y(u=s,v=v0),z(u=s,v=v0)),(s,urange[1], urange[2]),color='green',thickness=5)
for t in srange(vrange[1],vrange[2]+dv,dv):
UVPlot += line([(urange[1],t),(urange[2],t)],color='green')
SPlot += parametric_plot((x(u=s,v=t),y(u=s,v=t),z(u=s,v=t)),(s,urange[1], urange[2]),color='green')
show('U axis is green, V axis is red. Axes cross at (u0,v0)')
SPlot += parametric_plot3d((x,y,z),urange, vrange,color='grey',opacity=0.1)
UVPlot.show(aspect_ratio=1)
if showax:
SPlot += axes(4)
show(SPlot)
︡afe96a69-d9c5-4d7d-959e-03b4de33499d︡{"interact":{"controls":[{"control_type":"input-box","default":"v*cos(u)","label":"x(u,v)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"x","width":null},{"control_type":"input-box","default":"v*sin(u)","label":"y(u,v)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"y","width":null},{"control_type":"input-box","default":"u","label":"z(u,v)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"z","width":null},{"control_type":"input-box","default":6,"label":"nDivisions","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"n","width":20},{"control_type":"input-box","default":"(1/4*pi, 1.50000000000000)","label":"(u0,v0)","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"pt","width":40},{"control_type":"input-box","default":"(u, 0, 2*pi)","label":"u range","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"urange","width":36},{"control_type":"input-box","default":"(v, 1, 2)","label":"v range","nrows":1,"readonly":false,"submit_button":null,"type":null,"var":"vrange","width":36},{"control_type":"checkbox","default":false,"label":"show axes","readonly":false,"var":"showax"}],"flicker":false,"id":"e335a3f9-be16-424f-a55f-0841b792eed4","layout":[[["x",3,null],["y",3,null],["z",3,null],["n",3,null]],[["pt",3,null],["urange",3,null],["vrange",3,null],["showax",3,null]],[["",12,null]]],"style":"None"}}︡{"done":true}
︠c7ae4f27-d3b6-420d-b8b6-1178b26352a7s︠
ic=icosahedron(size=4, frame=false,frame_size=3, frame_color='red',opacity=0.5)
show(ic)
︡334e09de-b01e-4571-85d0-38aced3d56ef︡{"file":{"filename":"8a3a0aec-59a8-4293-ada8-cb830b0e1d53.sage3d","uuid":"8a3a0aec-59a8-4293-ada8-cb830b0e1d53"}}︡{"done":true}
︠efdd11f7-e45b-4856-bd13-59c08bcbdb7ds︠
ic
︡d8fe6d99-a556-4ce4-9ea3-54dda58ea8a5︡{"file":{"filename":"8a3a0aec-59a8-4293-ada8-cb830b0e1d53.sage3d","uuid":"8a3a0aec-59a8-4293-ada8-cb830b0e1d53"}}︡{"done":true}
︠9d91def3-0c08-4852-be0c-f9b4dc0e78f4︠
︡f59e2314-7dd8-4af0-a51a-d07b4b071745︡
︠bc915724-eb0a-40a7-997c-7f6ec142b65d︠