Sage Quickstart for Multivariable Calculus#

This Sage quickstart tutorial was developed for the MAA PREP Workshop “Sage: Using Open-Source Mathematics Software with Undergraduates” (funding provided by NSF DUE 0817071). It is licensed under the Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA).

Because much related material was covered in the calculus tutorial, this quickstart is designed to:

  • Give concrete examples of some things not covered there, without huge amounts of commentary, and

  • Remind the reader of a few of the things in that tutorial.

Vector Calculus#

In Sage, vectors are primarily linear algebra objects, but they are slowly becoming simultaneously analytic continuous functions.

Dot Product, Cross Product#

sage: v = vector(RR, [1.2, 3.5, 4.6])
sage: w = vector(RR, [1.7,-2.3,5.2])
sage: v*w
17.9100000000000

>>> from sage.all import *
>>> v = vector(RR, [RealNumber('1.2'), RealNumber('3.5'), RealNumber('4.6')])
>>> w = vector(RR, [RealNumber('1.7'),-RealNumber('2.3'),RealNumber('5.2')])
>>> v*w
17.9100000000000

sage: v.cross_product(w)
(28.7800000000000, 1.58000000000000, -8.71000000000000)

>>> from sage.all import *
>>> v.cross_product(w)
(28.7800000000000, 1.58000000000000, -8.71000000000000)

Lines and Planes#

Intersect \(x-2y-7z=6\) with \(\frac{x-3}{2}=\frac{y+4}{-3}=\frac{z-1}{1}\). One way to plot the equation of a line in three dimensions is with parametric_plot3d.

sage: # designed with intersection at t = 2, i.e. (7, -10, 3)
sage: var('t, x, y')
(t, x, y)
sage: line = parametric_plot3d([2*t+3, -3*t-4, t+1], (t, 0, 4),color='red')
sage: plane = plot3d((1/5)*(-12+x-2*y), (x, 4, 10), (y, -13,-7), opacity=0.5)
sage: intersect=point3d([7,-10,3],color='black',size=30)
sage: line+plane+intersect
Graphics3d Object

>>> from sage.all import *
>>> # designed with intersection at t = 2, i.e. (7, -10, 3)
>>> var('t, x, y')
(t, x, y)
>>> line = parametric_plot3d([Integer(2)*t+Integer(3), -Integer(3)*t-Integer(4), t+Integer(1)], (t, Integer(0), Integer(4)),color='red')
>>> plane = plot3d((Integer(1)/Integer(5))*(-Integer(12)+x-Integer(2)*y), (x, Integer(4), Integer(10)), (y, -Integer(13),-Integer(7)), opacity=RealNumber('0.5'))
>>> intersect=point3d([Integer(7),-Integer(10),Integer(3)],color='black',size=Integer(30))
>>> line+plane+intersect
Graphics3d Object

Vector-Valued Functions#

We can make vector-valued functions and do the usual analysis with them.

sage: var('t')
t
sage: r=vector((2*t-4, t^2, (1/4)*t^3))
sage: r
(2*t - 4, t^2, 1/4*t^3)

>>> from sage.all import *
>>> var('t')
t
>>> r=vector((Integer(2)*t-Integer(4), t**Integer(2), (Integer(1)/Integer(4))*t**Integer(3)))
>>> r
(2*t - 4, t^2, 1/4*t^3)

sage: r(t=5)
(6, 25, 125/4)

>>> from sage.all import *
>>> r(t=Integer(5))
(6, 25, 125/4)

The following makes the derivative also a vector-valued expression.

sage: velocity = r.diff(t) # velocity(t) = list(r.diff(t)) also would work
sage: velocity
(2, 2*t, 3/4*t^2)

>>> from sage.all import *
>>> velocity = r.diff(t) # velocity(t) = list(r.diff(t)) also would work
>>> velocity
(2, 2*t, 3/4*t^2)

Currently, this expression does not function as a function, so we need to substitute explicitly.

sage: velocity(t=1)
(2, 2, 3/4)

>>> from sage.all import *
>>> velocity(t=Integer(1))
(2, 2, 3/4)

sage: T=velocity/velocity.norm()

>>> from sage.all import *
>>> T=velocity/velocity.norm()

sage: T(t=1).n()
(0.683486126173409, 0.683486126173409, 0.256307297315028)

>>> from sage.all import *
>>> T(t=Integer(1)).n()
(0.683486126173409, 0.683486126173409, 0.256307297315028)

Here we compute the arclength between \(t=0\) and \(t=1\) by integrating the normalized derivative. As pointed out in the calculus tutorial, the syntax for numerical_integral is slightly nonstandard – we just put in the endpoints, not the variable.

sage: arc_length = numerical_integral(velocity.norm(), 0,1)
sage: arc_length
(2.3169847271197814, 2.572369791753217e-14)

>>> from sage.all import *
>>> arc_length = numerical_integral(velocity.norm(), Integer(0),Integer(1))
>>> arc_length
(2.3169847271197814, 2.572369791753217e-14)

We can also plot vector fields, even in three dimensions.

sage: x,y,z=var('x y z')
sage: plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors=['red','green','blue'])
Graphics3d Object

>>> from sage.all import *
>>> x,y,z=var('x y z')
>>> plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,Integer(0),pi), (y,Integer(0),pi), (z,Integer(0),pi),colors=['red','green','blue'])
Graphics3d Object

If we know a little vector calculus, we can also do line integrals. Here, based on an example by Ben Woodruff of BYU-Idaho, we compute a number of quantities in a physical three-dimensional setting.

Try to read through the entire code. We make an auxiliary function that will calculate \(\int (\text{integrand})\, dt\). We’ll use our formulas relating various differentials to \(dt\) to easily specify the integrands.

sage: density(x,y,z)=x^2+z
sage: r=vector([3*cos(t), 3*sin(t),4*t])
sage: tstart=0
sage: tend=2*pi
sage: t = var('t')
sage: t_range=(t,tstart,tend)
sage: def line_integral(integrand):
....:     return RR(numerical_integral((integrand).subs(x=r[0], y=r[1],z=r[2]), tstart, tend)[0])
sage: _ = var('x,y,z,t')
sage: r_prime = diff(r,t)
sage: ds=diff(r,t).norm()
sage: s=line_integral(ds)
sage: centroid_x=line_integral(x*ds)/s
sage: centroid_y=line_integral(y*ds)/s
sage: centroid_z=line_integral(z*ds)/s
sage: dm=density(x,y,z)*ds
sage: m=line_integral(dm)
sage: avg_density = m/s
sage: moment_about_yz_plane=line_integral(x*dm)
sage: moment_about_xz_plane=line_integral(y*dm)
sage: moment_about_xy_plane=line_integral(z*dm)
sage: center_mass_x = moment_about_yz_plane/m
sage: center_mass_y = moment_about_xz_plane/m
sage: center_mass_z = moment_about_xy_plane/m
sage: Ix=line_integral((y^2+z^2)*dm)
sage: Iy=line_integral((x^2+z^2)*dm)
sage: Iz=line_integral((x^2+y^2)*dm)
sage: Rx = sqrt(Ix/m)
sage: Ry = sqrt(Iy/m)
sage: Rz = sqrt(Iz/m)

>>> from sage.all import *
>>> __tmp__=var("x,y,z"); density = symbolic_expression(x**Integer(2)+z).function(x,y,z)
>>> r=vector([Integer(3)*cos(t), Integer(3)*sin(t),Integer(4)*t])
>>> tstart=Integer(0)
>>> tend=Integer(2)*pi
>>> t = var('t')
>>> t_range=(t,tstart,tend)
>>> def line_integral(integrand):
...     return RR(numerical_integral((integrand).subs(x=r[Integer(0)], y=r[Integer(1)],z=r[Integer(2)]), tstart, tend)[Integer(0)])
>>> _ = var('x,y,z,t')
>>> r_prime = diff(r,t)
>>> ds=diff(r,t).norm()
>>> s=line_integral(ds)
>>> centroid_x=line_integral(x*ds)/s
>>> centroid_y=line_integral(y*ds)/s
>>> centroid_z=line_integral(z*ds)/s
>>> dm=density(x,y,z)*ds
>>> m=line_integral(dm)
>>> avg_density = m/s
>>> moment_about_yz_plane=line_integral(x*dm)
>>> moment_about_xz_plane=line_integral(y*dm)
>>> moment_about_xy_plane=line_integral(z*dm)
>>> center_mass_x = moment_about_yz_plane/m
>>> center_mass_y = moment_about_xz_plane/m
>>> center_mass_z = moment_about_xy_plane/m
>>> Ix=line_integral((y**Integer(2)+z**Integer(2))*dm)
>>> Iy=line_integral((x**Integer(2)+z**Integer(2))*dm)
>>> Iz=line_integral((x**Integer(2)+y**Integer(2))*dm)
>>> Rx = sqrt(Ix/m)
>>> Ry = sqrt(Iy/m)
>>> Rz = sqrt(Iz/m)

Finally, we can display everything in a nice table. Recall that we use the r"stuff" syntax to indicate “raw” strings so that backslashes from LaTeX won’t cause trouble.

sage: html.table([
....:     [r"Density $\delta(x,y)$", density],
....:     [r"Curve $\vec r(t)$",r],
....:     [r"$t$ range", t_range],
....:     [r"$\vec r'(t)$", r_prime],
....:     [r"$ds$, a little bit of arclength", ds],
....:     [r"$s$ - arclength", s],
....:     [r"Centroid (constant density) $\left(\frac{1}{m}\int x\,ds,\frac{1}{m}\int y\,ds, \frac{1}{m}\int z\,ds\right)$", (centroid_x,centroid_y,centroid_z)],
....:     [r"$dm=\delta ds$ - a little bit of mass", dm],
....:     [r"$m=\int \delta ds$ - mass", m],
....:     [r"average density $\frac{1}{m}\int ds$" , avg_density.n()],
....:     [r"$M_{yz}=\int x dm$ - moment about $yz$ plane", moment_about_yz_plane],
....:     [r"$M_{xz}=\int y dm$ - moment about $xz$ plane", moment_about_xz_plane],
....:     [r"$M_{xy}=\int z dm$ - moment about $xy$ plane", moment_about_xy_plane],
....:     [r"Center of mass $\left(\frac1m \int xdm, \frac1m \int ydm, \frac1m \int z dm\right)$", (center_mass_x, center_mass_y, center_mass_z)],
....:     [r"$I_x = \int (y^2+z^2) dm$", Ix],[r"$I_y=\int (x^2+z^2) dm$", Iy],[mp(r"$I_z=\int (x^2+y^2)dm$"), Iz],
....:     [r"$R_x=\sqrt{I_x/m}$", Rx],[mp(r"$R_y=\sqrt{I_y/m}"), Ry],[mp(r"$R_z=\sqrt{I_z/m}"),Rz]
....:     ])

>>> from sage.all import *
>>> html.table([
...     [r"Density $\delta(x,y)$", density],
...     [r"Curve $\vec r(t)$",r],
...     [r"$t$ range", t_range],
...     [r"$\vec r'(t)$", r_prime],
...     [r"$ds$, a little bit of arclength", ds],
...     [r"$s$ - arclength", s],
...     [r"Centroid (constant density) $\left(\frac{1}{m}\int x\,ds,\frac{1}{m}\int y\,ds, \frac{1}{m}\int z\,ds\right)$", (centroid_x,centroid_y,centroid_z)],
...     [r"$dm=\delta ds$ - a little bit of mass", dm],
...     [r"$m=\int \delta ds$ - mass", m],
...     [r"average density $\frac{1}{m}\int ds$" , avg_density.n()],
...     [r"$M_{yz}=\int x dm$ - moment about $yz$ plane", moment_about_yz_plane],
...     [r"$M_{xz}=\int y dm$ - moment about $xz$ plane", moment_about_xz_plane],
...     [r"$M_{xy}=\int z dm$ - moment about $xy$ plane", moment_about_xy_plane],
...     [r"Center of mass $\left(\frac1m \int xdm, \frac1m \int ydm, \frac1m \int z dm\right)$", (center_mass_x, center_mass_y, center_mass_z)],
...     [r"$I_x = \int (y^2+z^2) dm$", Ix],[r"$I_y=\int (x^2+z^2) dm$", Iy],[mp(r"$I_z=\int (x^2+y^2)dm$"), Iz],
...     [r"$R_x=\sqrt{I_x/m}$", Rx],[mp(r"$R_y=\sqrt{I_y/m}"), Ry],[mp(r"$R_z=\sqrt{I_z/m}"),Rz]
...     ])

Functions of Several Variables#

This connects directly to other issues of multivariable functions.

How to view these was mostly addressed in the various plotting tutorials. Here is a reminder of what can be done.

sage: # import matplotlib.cm; matplotlib.cm.datad.keys()
sage: # 'Spectral', 'summer', 'blues'
sage: g(x,y)=e^-x*sin(y)
sage: contour_plot(g, (x, -2, 2), (y, -4*pi, 4*pi), cmap = 'Blues', contours=10, colorbar=True)
Graphics object consisting of 1 graphics primitive

>>> from sage.all import *
>>> # import matplotlib.cm; matplotlib.cm.datad.keys()
>>> # 'Spectral', 'summer', 'blues'
>>> __tmp__=var("x,y"); g = symbolic_expression(e**-x*sin(y)).function(x,y)
>>> contour_plot(g, (x, -Integer(2), Integer(2)), (y, -Integer(4)*pi, Integer(4)*pi), cmap = 'Blues', contours=Integer(10), colorbar=True)
Graphics object consisting of 1 graphics primitive

Partial Differentiation#

The following exercise is from Hass, Weir, and Thomas, University Calculus, Exercise 12.7.35. This function has a local minimum at \((4,-2)\).

sage: f(x, y) = x^2 + x*y + y^2 - 6*x + 2

>>> from sage.all import *
>>> __tmp__=var("x,y"); f = symbolic_expression(x**Integer(2) + x*y + y**Integer(2) - Integer(6)*x + Integer(2)).function(x,y)

Quiz: Why did we not need to declare the variables in this case?

sage: fx(x,y)= f.diff(x)
sage: fy(x,y) = f.diff(y)
sage: fx; fy
(x, y) |--> 2*x + y - 6
(x, y) |--> x + 2*y

>>> from sage.all import *
>>> __tmp__=var("x,y"); fx = symbolic_expression(f.diff(x)).function(x,y)
>>> __tmp__=var("x,y"); fy = symbolic_expression(f.diff(y)).function(x,y)
>>> fx; fy
(x, y) |--> 2*x + y - 6
(x, y) |--> x + 2*y

sage: f.gradient()
(x, y) |--> (2*x + y - 6, x + 2*y)

>>> from sage.all import *
>>> f.gradient()
(x, y) |--> (2*x + y - 6, x + 2*y)

sage: solve([fx==0, fy==0], (x, y))
[[x == 4, y == -2]]

>>> from sage.all import *
>>> solve([fx==Integer(0), fy==Integer(0)], (x, y))
[[x == 4, y == -2]]

sage: H = f.hessian()
sage: H(x,y)
[2 1]
[1 2]

>>> from sage.all import *
>>> H = f.hessian()
>>> H(x,y)
[2 1]
[1 2]

And of course if the Hessian has positive determinant and \(f_{xx}\) is positive, we have a local minimum.

sage: html("$f_{xx}=%s$"%H(4,-2)[0,0])
sage: html("$D=%s$"%H(4,-2).det())

>>> from sage.all import *
>>> html("$f_{xx}=%s$"%H(Integer(4),-Integer(2))[Integer(0),Integer(0)])
>>> html("$D=%s$"%H(Integer(4),-Integer(2)).det())

Notice how we were able to use many things we’ve done up to now to solve this.

  • Matrices

  • Symbolic functions

  • Solving

  • Differential calculus

  • Special formatting commands

  • And, below, plotting!

sage: plot3d(f,(x,-5,5),(y,-5,5))+point((4,-2,f(4,-2)),color='red',size=20)
Graphics3d Object

>>> from sage.all import *
>>> plot3d(f,(x,-Integer(5),Integer(5)),(y,-Integer(5),Integer(5)))+point((Integer(4),-Integer(2),f(Integer(4),-Integer(2))),color='red',size=Integer(20))
Graphics3d Object

Multiple Integrals and More#

Naturally, there is lots more that one can do.

sage: f(x,y)=x^2*y
sage: # integrate in the order dy dx
sage: f(x,y).integrate(y,0,4*x).integrate(x,0,3)
1944/5
sage: # another way to integrate, and in the opposite order too
sage: integrate( integrate(f(x,y), (x, y/4, 3)), (y, 0, 12) )
1944/5

>>> from sage.all import *
>>> __tmp__=var("x,y"); f = symbolic_expression(x**Integer(2)*y).function(x,y)
>>> # integrate in the order dy dx
>>> f(x,y).integrate(y,Integer(0),Integer(4)*x).integrate(x,Integer(0),Integer(3))
1944/5
>>> # another way to integrate, and in the opposite order too
>>> integrate( integrate(f(x,y), (x, y/Integer(4), Integer(3))), (y, Integer(0), Integer(12)) )
1944/5

sage: var('u v')
(u, v)
sage: surface = plot3d(f(x,y), (x, 0, 3.2), (y, 0, 12.3), color = 'blue', opacity=0.3)
sage: domain = parametric_plot3d([3*u, 4*(3*u)*v,0], (u, 0, 1), (v, 0,1), color = 'green', opacity = 0.75)
sage: image = parametric_plot3d([3*u, 4*(3*u)*v, f(3*u, 12*u*v)], (u, 0, 1), (v, 0,1), color = 'green', opacity = 1.00)
sage: surface+domain+image
Graphics3d Object

>>> from sage.all import *
>>> var('u v')
(u, v)
>>> surface = plot3d(f(x,y), (x, Integer(0), RealNumber('3.2')), (y, Integer(0), RealNumber('12.3')), color = 'blue', opacity=RealNumber('0.3'))
>>> domain = parametric_plot3d([Integer(3)*u, Integer(4)*(Integer(3)*u)*v,Integer(0)], (u, Integer(0), Integer(1)), (v, Integer(0),Integer(1)), color = 'green', opacity = RealNumber('0.75'))
>>> image = parametric_plot3d([Integer(3)*u, Integer(4)*(Integer(3)*u)*v, f(Integer(3)*u, Integer(12)*u*v)], (u, Integer(0), Integer(1)), (v, Integer(0),Integer(1)), color = 'green', opacity = RealNumber('1.00'))
>>> surface+domain+image
Graphics3d Object

Quiz: why did we need to declare variables this time?