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?