Lecture 12¶

Desargues Theorem¶

This is modified from Homework 4.

Recall that three points are colinear if they lie on a common line. A triple of points $\{a, b, c\}$ is a triangle if the points are not colinear.

Consider two triangles $\{a_0, b_0, c_0\}$ and $\{a_1, b_1, c_1\}$.

We say the triangles are in persective centrally if the three lines $\overline{a_0 a_1}$, $\overline{b_0 b_1}$, and $\overline{c_0 c_1}$ all pass through a common point called the center of perspectivity (these lines are coincident).

We say the triangles are in persective axially if the three points $\overline{a_0 b_0} \cap \overline{a_1 b_1}$, $\overline{b_0 c_0} \cap \overline{b_1 c_1}$, and $\overline{c_0 a_0} \cap \overline{c_1 a_1}$ all lie on the same line called the axis of perspectivity (the three points are colinear).

Desargue's Theorem says that two triangles are in persective centrally if and only if they are in in persective axially.

Proof¶

We use this function to construct symbolic points in the Projective plane.

In [1]:
def sym_vector(name):
    '''name should be a string.'''
    v = vector(SR, [var(f'{name}x'), var(f'{name}y'), var(f'{name}z')])
    return v
In [2]:
a = [sym_vector('a0'), sym_vector('a1')]
a
Out[2]:
[(a0x, a0y, a0z), (a1x, a1y, a1z)]
In [3]:
b = [sym_vector('b0'), sym_vector('b1')]
c = [sym_vector('c0'), sym_vector('c1')]
In [4]:
al = a[0].cross_product(a[1])
bl = b[0].cross_product(b[1])
cl = c[0].cross_product(c[1])

The triangles are in persective centrally if and only if the following determinant is zero:

In [5]:
M1 = matrix([al, bl, cl])
M1det = M1.det().expand()
In [6]:
pab = a[0].cross_product(b[0]).cross_product(a[1].cross_product(b[1]))
pac = a[0].cross_product(c[0]).cross_product(a[1].cross_product(c[1]))
pcb = c[0].cross_product(b[0]).cross_product(c[1].cross_product(b[1]))

The triangles are in persective axially if and only if the following determinant is zero:

In [7]:
M2 = matrix([pab, pac, pcb])
M2det = M2.det().expand()
In [8]:
(M2det/M1det).rational_simplify()
Out[8]:
-(((a0z*a1z*b0y - a0y*a1z*b0z)*b1y - (a0z*a1y*b0y - a0y*a1y*b0z)*b1z)*c0x - ((a0z*a1z*b0x - a0x*a1z*b0z)*b1y - (a0z*a1y*b0x - a0x*a1y*b0z)*b1z)*c0y + ((a0y*a1z*b0x - a0x*a1z*b0y)*b1y - (a0y*a1y*b0x - a0x*a1y*b0y)*b1z)*c0z)*c1x + (((a0z*a1z*b0y - a0y*a1z*b0z)*b1x - (a0z*a1x*b0y - a0y*a1x*b0z)*b1z)*c0x - ((a0z*a1z*b0x - a0x*a1z*b0z)*b1x - (a0z*a1x*b0x - a0x*a1x*b0z)*b1z)*c0y + ((a0y*a1z*b0x - a0x*a1z*b0y)*b1x - (a0y*a1x*b0x - a0x*a1x*b0y)*b1z)*c0z)*c1y - (((a0z*a1y*b0y - a0y*a1y*b0z)*b1x - (a0z*a1x*b0y - a0y*a1x*b0z)*b1y)*c0x - ((a0z*a1y*b0x - a0x*a1y*b0z)*b1x - (a0z*a1x*b0x - a0x*a1x*b0z)*b1y)*c0y + ((a0y*a1y*b0x - a0x*a1y*b0y)*b1x - (a0y*a1x*b0x - a0x*a1x*b0y)*b1y)*c0z)*c1z
In [9]:
rs = _
rs
Out[9]:
-(((a0z*a1z*b0y - a0y*a1z*b0z)*b1y - (a0z*a1y*b0y - a0y*a1y*b0z)*b1z)*c0x - ((a0z*a1z*b0x - a0x*a1z*b0z)*b1y - (a0z*a1y*b0x - a0x*a1y*b0z)*b1z)*c0y + ((a0y*a1z*b0x - a0x*a1z*b0y)*b1y - (a0y*a1y*b0x - a0x*a1y*b0y)*b1z)*c0z)*c1x + (((a0z*a1z*b0y - a0y*a1z*b0z)*b1x - (a0z*a1x*b0y - a0y*a1x*b0z)*b1z)*c0x - ((a0z*a1z*b0x - a0x*a1z*b0z)*b1x - (a0z*a1x*b0x - a0x*a1x*b0z)*b1z)*c0y + ((a0y*a1z*b0x - a0x*a1z*b0y)*b1x - (a0y*a1x*b0x - a0x*a1x*b0y)*b1z)*c0z)*c1y - (((a0z*a1y*b0y - a0y*a1y*b0z)*b1x - (a0z*a1x*b0y - a0y*a1x*b0z)*b1y)*c0x - ((a0z*a1y*b0x - a0x*a1y*b0z)*b1x - (a0z*a1x*b0x - a0x*a1x*b0z)*b1y)*c0y + ((a0y*a1y*b0x - a0x*a1y*b0y)*b1x - (a0y*a1x*b0x - a0x*a1x*b0y)*b1y)*c0z)*c1z
In [10]:
rs.factor()
Out[10]:
-(a0z*b0y*c0x - a0y*b0z*c0x - a0z*b0x*c0y + a0x*b0z*c0y + a0y*b0x*c0z - a0x*b0y*c0z)*(a1z*b1y*c1x - a1y*b1z*c1x - a1z*b1x*c1y + a1x*b1z*c1y + a1y*b1x*c1z - a1x*b1y*c1z)

Since $\{a_0, b_0, c_0\}$ is a triangle, the following determinant is non-zero:

In [11]:
D0 = matrix([a[0], b[0], c[0]]).det().expand()
D0
Out[11]:
-a0z*b0y*c0x + a0y*b0z*c0x + a0z*b0x*c0y - a0x*b0z*c0y - a0y*b0x*c0z + a0x*b0y*c0z
In [12]:
D1 = matrix([a[1], b[1], c[1]]).det().expand()
D1
Out[12]:
-a1z*b1y*c1x + a1y*b1z*c1x + a1z*b1x*c1y - a1x*b1z*c1y - a1y*b1x*c1z + a1x*b1y*c1z
In [13]:
(M2det + M1det*D0*D1).expand()
Out[13]:
0

3D graphics¶

References:

  • Computational Mathematics with SageMath: Chapter 4.2
  • The SageMath reference chapter on 3D Graphics

Plotting a function $f(x,y)$.¶

To make a somewhat interesting example, let us first define a triangle wave $g(t)$:

In [14]:
g(t) = 2*abs(t/2 - floor(t/2+1/2))
plot(g, (x, -3, 3), aspect_ratio=1)
Out[14]:
No description has been provided for this image
In [15]:
f(x,y) = g(x^2+y^2)
show(f)
\(\displaystyle \left( x, y \right) \ {\mapsto} \ 2 \, {\left| \frac{1}{2} \, x^{2} + \frac{1}{2} \, y^{2} - \left \lfloor \frac{1}{2} \, x^{2} + \frac{1}{2} \, y^{2} + \frac{1}{2} \right \rfloor \right|}\)
In [16]:
r = 3/2
plot3d(f, (x, -r, r), (y, -r, r))
Out[16]:

It looks like there are some defects. These can be reduced by increasing the number of points used:

In [17]:
plot3d(f, (x, -r, r), (y, -r, r), plot_points=200)
Out[17]:

The option plot_points is one of several options to plot3d. Here from the in-notebook documentation plot3d?

  • "adaptive" -- (default: False) whether to use adaptive refinement to draw the plot (slower, but may look better). This option does NOT work in conjunction with a transformation (see below).

  • "mesh" -- bool (default: False) whether to display mesh grid lines

  • "dots" -- bool (default: False) whether to display dots at mesh grid points

  • "plot_points" -- (default: "automatic") initial number of sample points in each direction; an integer or a pair of integers

The mesh option is useful for seeing why the initial plot looked bumpy:

In [18]:
plot3d(f, (x, -r, r), (y, -r, r), adaptive=True)
Out[18]:
In [19]:
plot3d(f, (x, -r, r), (y, -r, r), mesh=True)
Out[19]:

There are also different viewers.

The default viewer seems to be Three.js JavaScript WebGL Renderer. That link provides a lot of further options to customize your plot:

  • aspect_ratio -- (default: [1,1,1]) list or tuple of three numeric values

  • auto_scaling -- (default: [False, False, False]) list or tuple of three booleans; set to True to automatically scale down the corresponding direction if it is too large

  • axes -- (default: False) Boolean determining whether coordinate axes are drawn

  • axes_labels -- (default: ['x','y','z']) list or tuple of three strings; set to False to remove all labels

  • axes_labels_style -- (default: None) list of three dicts, one per axis, or a single dict controlling all three axes; supports the same styling options as :func:~sage.plot.plot3d.shapes2.text3d such as color, opacity, fontsize, fontweight, fontstyle, and fontfamily

  • color -- (default: 'blue') color of the 3D object

  • decimals -- (default: 2) integer determining decimals displayed in labels

  • depth_write -- (default: True for opaque surfaces, False for transparent surfaces) whether to write the surface's depth into the depth buffer for the purpose of occluding objects behind it

  • frame -- (default: True) Boolean determining whether frame is drawn

  • online -- (default: False) Boolean determining whether the local standard package files are replaced by links to an online content delivery network

  • opacity -- (default: 1) numeric value for transparency of lines and surfaces

  • page_title -- (default: None) string containing the title of the generated HTML page; often displayed in the browser window's title bar, its tab list, and/or the operating system's task bar

  • projection -- (default: 'perspective') the type of camera projection to use; 'perspective' or 'orthographic'

  • radius -- (default: None) numeric value for radius of lines; use to render lines thicker than available using thickness or on Windows platforms where thickness is ignored

  • render_order -- (default: 0) numeric value for rendering order of transparent surfaces; objects render from lowest to highest value ensuring that lower-valued objects render completely

  • single_side -- (default: False) Boolean determining whether both sides of a surface material are rendered; set to True to reduce rendering artifacts for closed transparent surfaces

  • theme -- (default: 'light') the color scheme to use for the scene and user interface; 'light' or 'dark'

  • thickness -- (default: 1) numeric value for thickness of lines

  • viewpoint -- (default: None) list or tuple of the form [[x,y,z],angle] setting the initial viewpoint of the scene, where angle is in degrees; can be determined using the 'Get Viewpoint' option of the information menu

In addition, there are some animation options.

In [20]:
plot3d(f, (x,-r,r), (y,-r,r), viewer='threejs', 
       plot_points=200, frame=False, opacity=0.9)
Out[20]:
In [21]:
plot3d(f, (x,-r,r), (y,-r,r), viewer='threejs', 
       plot_points=200, frame=False, opacity=0.9)
Out[21]:

The Tachyon 3D Ray Tracer is a non-interactive but can generate high quality images. Requires some playing with options to get a nice plot.

Here we demonstrate the non-interactive tachyon viewer. It presumably looks a bit better on a smooth surface.

In [22]:
plot3d(f, (x,-r,r), (y,-r,r), 
       antialiasing=True, raydepth=3, opacity = 1,
       camera_position = (3,0.6,1.5),
       plot_points=300, color='purple', 
       frame=False, viewer='tachyon', specular=0.5,
       viewpoint=[[-0.3644,5,-0.8207],111.59])
Out[22]:
No description has been provided for this image

Points in space¶

Lists of points can be dispayed in 3D with point3d:

In [23]:
point3d([(random(),random(),random()) for i in range(50)])
Out[23]:

Spheres¶

A sphere can be constructed with a center and a radius:

In [24]:
sphere((0,0,0), 1)
Out[24]:

A higher quality version of the plot above.

In [25]:
plt = None
for i in range(50):
    plt2 = sphere([random(), random(), random()], 0.03)
    if plt is None:
        plt = plt2
    else:
        plt += plt2
plt
Out[25]:

Boxes¶

Some other shapes are obtainable by running an import. For example Boxes.

In [26]:
from sage.plot.plot3d.shapes import Box
In [27]:
Box([1,2,3])
Out[27]:
In [28]:
Box([1,2,3]).translate([1,2,3])
Out[28]:

Block stacking¶

Say a block has width one. How far can we push slide it off the table until it falls? Position the table so that the surface has negative $x$-coordinate with the edge at $x=0$.

In [29]:
var('x')
block_sticks_out = x
center_of_mass = (x-1/2)
equation = center_of_mass==0
solve(equation, x)
Out[29]:
[x == (1/2)]

Now suppose we have two blocks. The top one can be pushed halfway off the bottom block. How far can we shift the bottom block before the stack falls?

In [30]:
var('x')
bottom_block = x
top_block = x + 1/2
center_of_mass = 1/2 * ((bottom_block-1/2)+(top_block-1/2))
equation = center_of_mass==0
solve(equation, x)
Out[30]:
[x == (1/4)]

Now if we have three blocks:

In [31]:
var('x')
bottom_block = x
middle_block = bottom_block + 1/4
top_block = middle_block + 1/2
center_of_mass = 1/2 * ((bottom_block-1/2)+(middle_block-1/2)+(top_block-1/2))
equation = center_of_mass==0
solve(equation, x)
Out[31]:
[x == (1/6)]
In [32]:
from sage.plot.plot3d.shapes import Box
In [33]:
def block_stack(n, dims=[8, 4, 1] ):
    dims = 1/2*vector(QQ, dims)
    def box(i):
        if i%2 == 0:
            return Box(list(dims), color='yellow').translate([-dims[0], 0, dims[2]])
        else:
            return Box(list(dims), color='blue').translate([-dims[0], 0, dims[2]])
    plt = Box([dims[0], dims[0], 1], color='brown').translate([-dims[0], 0, -1])
    delta_x = 0
    shifts = []
    for i in range(n):
        shifts.append( 1/(2*(n-i)))
        delta_x += 2*dims[0]/(2*(n-i))
        delta_z = i*2*dims[2]
        b = box(i).translate((delta_x, 0, delta_z))
        if plt is None:
            plt = b
        else:
            plt += b
    show(shifts)
    return plt
block_stack(3)
\(\displaystyle \left[\frac{1}{6}, \frac{1}{4}, \frac{1}{2}\right]\)
Out[33]:
In [34]:
@interact
def interactive_block_stack(n='1'):
    n = ZZ(n)
    show(block_stack(n), frame=False, viewpoint=[[-0.085,-0.6331,-0.7694],168.2])

Interactive function <function interactive_block_stack at 0x7f6a23963ec0> with 1 widget
  n: Text(value='1', description='n')

Wikipedia has an article on the block stacking problem.

Lines in space¶

In [35]:
line3d([(0,0,0),(1,1,1)])
Out[35]:
In [36]:
line3d([(0,0,0),(1,1,1)], thickness=5)
Out[36]:
In [37]:
line3d([(0,0,0),(1,1,1)], radius=0.1)
Out[37]:
In [38]:
random()
Out[38]:
0.03725204456474551
In [39]:
def random_walk(n):
    steps = [
        vector([1,0,0]),
        vector([-1,0,0]),
        vector([0,1,0]),
        vector([0,-1,0]),
        vector([0,0,1]),
        vector([0,0,-1]),
    ]
    pos = vector([0,0,0])
    path = [pos]
    for i in range(n):
        pos += steps[floor(6*random())]
        path.append(pos)
    return path
In [40]:
random_walk(5)
Out[40]:
[(0, 0, 0), (0, 0, -1), (-1, 0, -1), (-1, 0, -2), (0, 0, -2), (0, -1, -2)]
In [41]:
line3d(random_walk(50), radius=0.05)
Out[41]:
In [42]:
line3d(random_walk(5000), thickness=3)
Out[42]:

Parametric curves¶

In [43]:
p(t) = (t*cos(t), t*sin(t), t)
In [44]:
parametric_plot3d(p, (t,0, 10*pi), thickness=3, plot_points=200)
Out[44]:

Parametric surfaces¶

We plot a helicoid, which in cylinderical coordinates is given by $\theta=z$ (or here $\theta=6z$)

In [45]:
p(r, z) = (r*cos(6*z), r*sin(6*z), z)
show(p)
\(\displaystyle \left( r, z \right) \ {\mapsto} \ \left(r \cos\left(6 \, z\right),\,r \sin\left(6 \, z\right),\,z\right)\)
In [46]:
parametric_plot3d(p, (r, -1, 1), (z, 0, pi/3))
Out[46]:

There are more examples in the Sage documentation for parametric_plot_3d of the use of color maps.

In [47]:
def coloring_function(r,z):
    return abs(sin(6*z))
cm = colormaps.rainbow
In [48]:
parametric_plot3d(p, (r, -1, 1), (z, 0, pi/3), color=(coloring_function, cm), plot_points=200)
Out[48]:

Least squares in 2D:¶

Consider the function $f:\mathbb R^2 \to \mathbb R$ given by $$f(x,y) = \cos(x) \cos(y).$$ Find the best approximation $g(x,y)$ from the span of the collection of functions: $$\{1, x, y, x^2, y^2, xy, \sin(x), \sin(y)\}$$ in the sense that it minimizes the integral: $$\iint_S \big(f(x,y)-g(x,y)\big)^2~dx~dy$$ over the square $S=[-r,r] \times [-r, r]$ with $r=\frac{3}{2}$.

Remark: In class, I made the mistake of typing $f$ incorrectly, see below. So, the solution is worked out in this case, which is much simpler. See the lecture notes for a more challenging example.

In [49]:
f(x,y) = cos(x)*sin(y)
In [50]:
basis = [1, x, y, x^2, y^2, x*y, sin(x), sin(y)]
basis
Out[50]:
[1, x, y, x^2, y^2, x*y, sin(x), sin(y)]
In [51]:
c = [var(f'c{i}') for i in range(8)]
c
Out[51]:
[c0, c1, c2, c3, c4, c5, c6, c7]
In [52]:
g(x,y) = sum([c[i]*basis[i] for i in range(8)])
g
Out[52]:
(x, y) |--> c3*x^2 + c5*x*y + c4*y^2 + c1*x + c2*y + c6*sin(x) + c7*sin(y) + c0
In [53]:
sq = ( f(x,y)-g(x,y) )^2
sq
Out[53]:
(c3*x^2 + c5*x*y + c4*y^2 + c1*x + c2*y + c6*sin(x) + c7*sin(y) - cos(x)*sin(y) + c0)^2
In [54]:
r = 3/2
In [55]:
I = sq.integral(x, -r, r).integral(y, -r, r)
I
Out[55]:
-6*c1*c6*(3*cos(3/2) - 2*sin(3/2)) - 3/2*c6^2*(sin(3) - 3) - 3/2*c7^2*(sin(3) - 3) + 9*c0^2 + 27/4*c1^2 + 4*(3*cos(3/2)*sin(3/2) - 2*sin(3/2)^2)*c2 + 27/4*c2^2 + 27/2*c0*c3 + 729/80*c3^2 + 27/8*(4*c0 + 3*c3)*c4 + 729/80*c4^2 + 81/16*c5^2 - 2*(3*c2*(3*cos(3/2) - 2*sin(3/2)) - sin(3)*sin(3/2) + 3*sin(3/2))*c7 - 1/4*sin(3)^2 + 9/4
In [56]:
I.coefficient(c[1]*c[6])
Out[56]:
-18*cos(3/2) + 12*sin(3/2)
In [57]:
M = matrix([
    [ I.coefficient(c[i]*c[j]) if i==j else 1/2*I.coefficient(c[i]*c[j]) for i in range(8)]
    for j in range(8)
])
show(M)
\(\displaystyle \left(\begin{array}{rrrrrrrr} 9 & 0 & 0 & \frac{27}{4} & \frac{27}{4} & 0 & 0 & 0 \\ 0 & \frac{27}{4} & 0 & 0 & 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 \\ 0 & 0 & \frac{27}{4} & 0 & 0 & 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) \\ \frac{27}{4} & 0 & 0 & \frac{729}{80} & \frac{81}{16} & 0 & 0 & 0 \\ \frac{27}{4} & 0 & 0 & \frac{81}{16} & \frac{729}{80} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{81}{16} & 0 & 0 \\ 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 & 0 & 0 & 0 & -\frac{3}{2} \, \sin\left(3\right) + \frac{9}{2} & 0 \\ 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 & 0 & 0 & 0 & -\frac{3}{2} \, \sin\left(3\right) + \frac{9}{2} \end{array}\right)\)
In [58]:
cv = vector(c)
cv
Out[58]:
(c0, c1, c2, c3, c4, c5, c6, c7)
In [59]:
I2 = (I - cv*M*cv).expand()
I2
Out[59]:
12*c2*cos(3/2)*sin(3/2) + 2*c7*sin(3)*sin(3/2) - 8*c2*sin(3/2)^2 - 1/4*sin(3)^2 - 6*c7*sin(3/2) + 9/4
In [60]:
lin = vector([I2.coefficient(c[i]) for i in range(8)])
lin
Out[60]:
(0, 0, 12*cos(3/2)*sin(3/2) - 8*sin(3/2)^2, 0, 0, 0, 0, 2*sin(3)*sin(3/2) - 6*sin(3/2))
In [61]:
const = (I2 - lin*cv).expand()
const
Out[61]:
-1/4*sin(3)^2 + 9/4
In [62]:
( I - (cv*M*cv + lin*cv + const) ).expand()
Out[62]:
0
In [63]:
show(M)
\(\displaystyle \left(\begin{array}{rrrrrrrr} 9 & 0 & 0 & \frac{27}{4} & \frac{27}{4} & 0 & 0 & 0 \\ 0 & \frac{27}{4} & 0 & 0 & 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 \\ 0 & 0 & \frac{27}{4} & 0 & 0 & 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) \\ \frac{27}{4} & 0 & 0 & \frac{729}{80} & \frac{81}{16} & 0 & 0 & 0 \\ \frac{27}{4} & 0 & 0 & \frac{81}{16} & \frac{729}{80} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{81}{16} & 0 & 0 \\ 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 & 0 & 0 & 0 & -\frac{3}{2} \, \sin\left(3\right) + \frac{9}{2} & 0 \\ 0 & 0 & -9 \, \cos\left(\frac{3}{2}\right) + 6 \, \sin\left(\frac{3}{2}\right) & 0 & 0 & 0 & 0 & -\frac{3}{2} \, \sin\left(3\right) + \frac{9}{2} \end{array}\right)\)
In [64]:
gradient_eq = [I.derivative(c[i])==0 for i in range(8)]
gradient_eq
Out[64]:
[18*c0 + 27/2*c3 + 27/2*c4 == 0,
 -6*c6*(3*cos(3/2) - 2*sin(3/2)) + 27/2*c1 == 0,
 -6*c7*(3*cos(3/2) - 2*sin(3/2)) + 12*cos(3/2)*sin(3/2) - 8*sin(3/2)^2 + 27/2*c2 == 0,
 27/2*c0 + 729/40*c3 + 81/8*c4 == 0,
 27/2*c0 + 81/8*c3 + 729/40*c4 == 0,
 81/8*c5 == 0,
 -6*c1*(3*cos(3/2) - 2*sin(3/2)) - 3*c6*(sin(3) - 3) == 0,
 -6*c2*(3*cos(3/2) - 2*sin(3/2)) - 3*c7*(sin(3) - 3) + 2*sin(3)*sin(3/2) - 6*sin(3/2) == 0]
In [65]:
sol = solve(gradient_eq, c, solution_dict=True)[0]
In [66]:
g0 = g.substitute(sol)
g0
Out[66]:
(x, y) |--> 2/3*sin(3/2)*sin(y)
In [67]:
plot3d(f, (x, -r, r), (y, -r, r)) + plot3d(g0, (x, -r, r), (y, -r, r), color='red')
Out[67]: