Sage Tutorial for Advanced Plotting - Part 1

This Sage worksheet is one of the tutorials developed for the MAA PREP Workshop "Sage: Using Open-Source Mathematics Software with Undergraduates" (funding provided by NSF DUE 0817071).

Thanks to Sage's integration of projects like matplotlib, Jmol, and Tachyon, Sage has comprehensive plotting capabilities. This pair of worksheets aims to go more in depth with them than previous tutorials.   This worksheet consists of the following sections:

This tutorial assumes that one is familiar with the basics of Sage, such as evaluating a cell by clicking the "evaluate" link, or by pressing Shift-Enter (hold down Shift while pressing the Enter key).

Cartesian Plots

A simple quadratic is easy.

{{{id=3| plot(x^2, (x,-2,2)) /// }}}

You can combine "plot objects" by adding them.

{{{id=6| regular = plot(x^2, (x,-2,2), color= 'purple') skinny = plot(4*x^2, (x,-2,2), color = 'green') regular + skinny /// }}}

Problem:  Plot a green $y=\sin(x)$ together with a red $y=2\,\cos(x)$.  (Hint: you can use "pi" as part of your range.)

{{{id=12| /// }}}

Boundaries of a plot can be specified, in addition to the overall size.

{{{id=9| plot(1+e^(-x^2), xmin=-2, xmax=2, ymin=0, ymax=2.5, figsize=10) /// }}}

Problem:  Plot $y=5+3\,\sin(4x)$ with suitable boundaries.

{{{id=52| /// }}}

You can add lots of extra information.

{{{id=11| exponential = plot(1+e^(-x^2), xmin=-2, xmax=2, ymin=0, ymax=2.5) max_line = plot(2, xmin=-2, xmax=2, linestyle='-.', color = 'red') min_line = plot(1, xmin=-2, xmax=2, linestyle=':', color = 'red') exponential + max_line + min_line /// }}}

You can fill regions with transparent color, and thicken the curve.  This example uses several options to fine-tune our graphic.

{{{id=15| exponential = plot(1+e^(-x^2), xmin=-2, xmax=2, ymin=0, ymax=2.5, fill=0.5, fillcolor='grey', fillalpha=0.3) min_line = plot(1, xmin=-2, xmax=2, linestyle='-', thickness= 6, color = 'red') exponential + min_line /// }}}

Problem:  Create a plot showing the cross-section area for the following solid of revolution problem:  Consider the area bounded by $y=x^2-3x+6$ and the line $y=4$.  Find the volume created by rotating this area around the line $y=1$.

{{{id=19| /// }}}

Parametric Plots

A parametric plot needs a list of two functions of the parameter; in Sage, we use square brackets to delimit the list.  Notice also that we must declare $t$ as a variable first.  Because the graphic is slightly wider than it is tall, we use the 'aspect_ratio' option (such options are called keywords) to ensure the axes are correct for how we want to view this object.

{{{id=17| var('t') parametric_plot([cos(t) + 3 * cos(t/9), sin(t) - 3 * sin(t/9)], (t, 0, 18*pi), fill = True, aspect_ratio=1) /// }}}

Problem: These parametric equations will create a hypocycloid,

$$x(t)=17\cos(t)+3\cos(17t/3)$$

$$y(t)=17\sin(t)-3\sin(17t/3)$$

Create this as a parametric plot.

{{{id=34| /// }}}

Sage automatically plots a 2d or 3d plot, and a curve or a surface, depending on how many variables and coordinates you specify.

{{{id=89| var('t') parametric_plot((t^2,sin(t)), (t,0,pi)) /// }}} {{{id=80| var('t') parametric_plot((t^2,sin(t),cos(t)), (t,0,pi)) /// }}} {{{id=51| var('t,r') parametric_plot((t^2,sin(r*t),cos(r*t)), (t,0,pi),(r,-1,1)) /// }}}

Polar Plots

Sage can also do polar plots.

{{{id=28| polar_plot(2 + 2*cos(x), (x, 0, 2*pi), color=hue(0.5), thickness=4, aspect_ratio=1) /// }}}

Although they aren't essential, many of these examples try to demonstrate things like coloring, fills, and shading to give you a sense of the possibilities.

More than one polar curve can be specified in a list (square brackets).  Notice the automatic graded shading of the fill color.

{{{id=29| var('t') polar_plot([cos(4*t) + 1.5, 0.5 * cos(4*t) + 2.5], (t, 0, 2*pi), color='black', thickness=2, fill=True, fillcolor='orange', aspect_ratio=1) /// }}}

Problem: Create a plot for the following problem. Find the area that is inside the circle $r=2$, but outside the cardiod $2+2\cos(\theta)$.

{{{id=31| /// }}} {{{id=55| /// }}} {{{id=54| /// }}}

Interactive Graphics

It may be of interest to see all these things put together in a very nice pedagogical graphic.   Even though this is fairly advanced (and so we hide the code) it is not as difficult as you might think to put together.

{{{id=47| %hide %auto html('

Sine and unit circle (by Jurgis Pralgauskis)

inspired by this video' ) # http://www.sagemath.org/doc/reference/sage/plot/plot.html radius = 100 # scale for radius of "unit" circle graph_params = dict(xmin = -2*radius, xmax = 360, ymin = -(radius+30), ymax = radius+30, aspect_ratio=1, axes = False ) def sine_and_unit_circle( angle=30, instant_show = True, show_pi=True ): ccenter_x, ccenter_y = -radius, 0 # center of cirlce on real coords sine_x = angle # the big magic to sync both graphs :) current_y = circle_y = sine_y = radius * sin(angle*pi/180) circle_x = ccenter_x + radius * cos(angle*pi/180) graph = Graphics() # we'll put unit circle and sine function on the same graph # so there will be some coordinate mangling ;) # CIRCLE unit_circle = circle((ccenter_x, ccenter_y), radius, color="#ccc") # SINE x = var('x') sine = plot( radius * sin(x*pi/180) , (x, 0, 360), color="#ccc" ) graph += unit_circle + sine # CIRCLE axis # x axis graph += arrow( [-2*radius, 0], [0, 0], color = "#666" ) graph += text("$(1, 0)$", [-16, 16], color = "#666") # circle y axis graph += arrow( [ccenter_x,-radius], [ccenter_x, radius], color = "#666" ) graph += text("$(0, 1)$", [ccenter_x, radius+15], color = "#666") # circle center graph += text("$(0, 0)$", [ccenter_x, 0], color = "#666") # SINE x axis graph += arrow( [0,0], [360, 0], color = "#000" ) # let's set tics # or http://aghitza.org/posts/tweak_labels_and_ticks_in_2d_plots_using_matplotlib/ # or wayt for http://trac.sagemath.org/sage_trac/ticket/1431 # ['$-\pi/3$', '$2\pi/3$', '$5\pi/3$'] for x in range(0, 361, 30): graph += point( [x, 0] ) angle_label = ". $%3d^{\circ}$ " % x if show_pi: angle_label += " $(%s\pi) $"% x/180 graph += text(angle_label, [x, 0], rotation=-90, vertical_alignment='top', fontsize=8, color="#000" ) # CURRENT VALUES # SINE -- y graph += arrow( [sine_x,0], [sine_x, sine_y], width=1, arrowsize=3) graph += arrow( [circle_x,0], [circle_x, circle_y], width=1, arrowsize=3) graph += line(([circle_x, current_y], [sine_x, current_y]), rgbcolor="#0F0", linestyle = "--", alpha=0.5) # LABEL on sine graph += text("$(%d^{\circ}, %3.2f)$"%(sine_x, float(current_y)/radius), [sine_x+30, current_y], color = "#0A0") # ANGLE -- x # on sine graph += arrow( [0,0], [sine_x, 0], width=1, arrowsize=1, color='red') # on circle graph += disk( (ccenter_x, ccenter_y), float(radius)/4, (0, angle*pi/180), color='red', fill=False, thickness=1) graph += arrow( [ccenter_x, ccenter_y], [circle_x, circle_y], rgbcolor="#cccccc", width=1, arrowsize=1) if instant_show: show (graph, **graph_params) return graph ##################### # make Interaction ###################### @interact def _( angle = slider([0..360], default=30, step_size=5, label="Pasirinkite kampÄ…: ", display_value=True) ): sine_and_unit_circle(angle, show_pi = False) ///

Sine and unit circle (by Jurgis Pralgauskis)

inspired by this video
Pasirinkite kampÄ…:  
}}}

Plotting Data Points

Sometimes one wishes to simply plot data.  Here, we demonstrate several ways of plotting points and data via the simple approximation to the Fibonacci numbers given by

$$F_n=\frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n\; ,$$

which is quite good after about $n=5$.

First, we notice that the Fibonacci numbers are built in.

{{{id=61| fibonacci_sequence(6) /// }}} {{{id=62| list(fibonacci_sequence(6)) /// [0, 1, 1, 2, 3, 5] }}} {{{id=60| list(enumerate(fibonacci_sequence(6))) /// [(0, 0), (1, 1), (2, 1), (3, 2), (4, 3), (5, 5)] }}}

In this cell, we just define the numbers and coordinate pairs we are about to plot.

{{{id=56| fibonnaci = list(enumerate(fibonacci_sequence(6))) f(n)=(1/sqrt(5))*((1+sqrt(5))/2)^n asymptotic = [(i, f(i)) for i in range(6)] /// }}} {{{id=137| fibonnaci /// [(0, 0), (1, 1), (2, 1), (3, 2), (4, 3), (5, 5)] }}} {{{id=138| asymptotic /// [(0, 1/5*sqrt(5)), (1, 1/10*(sqrt(5) + 1)*sqrt(5)), (2, 1/20*(sqrt(5) + 1)^2*sqrt(5)), (3, 1/40*(sqrt(5) + 1)^3*sqrt(5)), (4, 1/80*(sqrt(5) + 1)^4*sqrt(5)), (5, 1/160*(sqrt(5) + 1)^5*sqrt(5))] }}}

And here we plot not just the two sets of points, but also use several of the documented options for plotting points; the choice of markers for different data is particularly helpful for generating publishable graphics.

{{{id=36| fib_plot=scatter_plot(fibonnaci, facecolor='red', marker='o',markersize=40) asy_plot = line(asymptotic, marker='D',color='black',thickness=2) show(fib_plot+asy_plot, aspect_ratio=1) /// }}} {{{id=45| /// }}}

Contour Plots

Contour plotting can be very useful when trying to get a handle on multivariable functions, as well as modeling.  The basic syntax is essentially the same as for 3D plotting - simply an extension of the 2D plotting syntax.

{{{id=118| f(x,y)=y^2+1-x^3-x contour_plot(f, (x,-pi,pi), (y,-pi,pi)) /// }}}

We can change colors, specify contours, label curves, and many other things.  When there are many levels, the 'colorbar' keyword becomes quite useful for keeping track of them.  Notice that, as opposed to many other options, it can only be 'True' or 'False' (corresponding to whether it appears or does not appear).

{{{id=121| contour_plot(f, (x,-pi,pi), (y,-pi,pi),colorbar=True,labels=True) /// }}}

This example is fairly self-explanatory, but demonstrates the power of formatting, labeling, and the wide variety of built-in color gradations (colormaps or 'cmap's).  The strange-looking construction corresponding to 'label_fmt' is a Sage/Python data type called a dictionary, and turns out to be useful for more advanced Sage use; it consists of pairs connected by a colon, all inside curly braces.

{{{id=117| contour_plot(f, (x,-pi,pi), (y,-pi,pi), contours=[-4,0,4], fill=False, cmap='cool', labels=True, label_inline=True, label_fmt={-4:"low", 0:"medium", 4: "hi"}, label_colors='black') /// }}}

Implicit plots are a special type of contour plot (they just plot the zero contour).

{{{id=139| f(x,y) /// -x^3 + y^2 - x + 1 }}} {{{id=116| implicit_plot(f(x,y)==0,(x,-pi,pi),(y,-pi,pi)) /// }}} {{{id=133| /// }}}

A density plot is like a contour plot, but without discrete levels.

{{{id=132| x,y = var('x,y') density_plot(f, (x, -2, 2), (y, -2, 2)) /// }}}

In the last session, Ron brought up that contour plots can be a little misleading. Here we combine a density plot and contour plot to show even better what is happening with the function.

{{{id=134| density_plot(f,(x,-2,2),(y,-2,2))+contour_plot(f,(x,-2,2),(y,-2,2),fill=False,labels=True,cmap='jet') /// }}}

Note: mention something about the options

{{{id=142| contour_plot.options /// {'labels': False, 'linestyles': None, 'frame': True, 'axes': False, 'plot_points': 100, 'linewidths': None, 'colorbar': False, 'contours': None, 'fill': True} }}} {{{id=141| contour_plot.options["fill"]=False /// }}} {{{id=143| contour_plot(f,(x,-2,2),(y,-2,2)) /// }}} {{{id=145| contour_plot.reset() /// }}} {{{id=144| contour_plot.options /// {'labels': False, 'linestyles': None, 'frame': True, 'axes': False, 'plot_points': 100, 'linewidths': None, 'colorbar': False, 'contours': None, 'fill': False} }}} {{{id=140| /// }}}

3d contour plots are given by the implicit_plot3d function.

{{{id=115| var('x y z') implicit_plot3d(-(cos(x) + cos(y) + cos(z)), (x, -4, 4), (y, -4, 4), (z, -4, 4),contour = [0,-2],opacity=0.4) /// }}}

Problem

  1. Plot the surface $z=f(x,y)=x^2y$ on a reasonable domain.
  2. In a contrasting color, plot the tangent plane to this surface at $x=1$, $y=3$.
  3. Adjust the opacity of each to make viewing easier.
  4. Place a "point" graphics object in a third color at the location of the tangency.
{{{id=125| /// }}}

3D Plotting

Plotting surfaces works as you might expect.  The JMOL Java applet has lots of built-in features:  zooming, rotating, stereoscopic viewing, automatic spin.  Note the need to specify the variables.

{{{id=38| var('x y') plot3d(x^2+3*y^2, (x,-2,2),(y,-2,2)) /// }}}

Here is an example illustrating multiple integration, built up from two surfaces and various graphics primitives.

{{{id=40| from sage.plot.plot3d.shapes2 import frame3d var('x y') up = plot3d(x^2+y^2, (x, -2, 2), (y, -2, 2), color='red', opacity = 0.5) down = plot3d(8-x^2-y^2, (x, -2, 2), (y, -2, 2), color='purple', opacity = 0.5) intersect = circle((0,0), 2, thickness=4, color='black').plot3d(z=4) shadow = circle((0,0), 2, thickness=4, color='black').plot3d(z=0) base = point((1,1,0), color='green') epsilon=0.1 column = frame3d((1-epsilon,1-epsilon,2),(1+epsilon,1+epsilon,6), thickness=3, color='green') cube = frame3d((1-epsilon,1-epsilon,5-2*epsilon),(1+epsilon,1+epsilon,5+2*epsilon), thickness=3, color='black') show(up+down+intersect+shadow+base+column+cube, aspect_ratio=[2,2,1]) /// }}}

Plotting with transformations

You can easily apply a transformation to your 3d function.  Built-in transformations include spherical and cylindrical transformations.

{{{id=114| var('phi,theta') spherical_plot3d((2+cos(2*theta))*(phi+1),(theta,0,2*pi),(phi,0,pi),color='red') /// }}} {{{id=112| cylindrical_plot3d(theta*cosh(z),(theta,0,2*pi),(z,-2,2)) /// }}}

You can make your own transformations and apply them easily.  Here we define the parabolic cylindrical transformation by $$x=wv,\quad y=(v^2-w^2)/2,\quad z=u$$ and plot the function $w=4+u$ (Sage knows that we are plotting $w=2$ since we specified ranges for $u$ and $v$).

{{{id=130| var('u,v,w') parabolic_cylindrical=(w*v,(v^2-w^2)/2,u) #(x,y,z) coordinates p=plot3d(4+u,(u,-pi,pi),(v,-pi,pi),transformation=parabolic_cylindrical) /// }}} {{{id=129| p /// }}} {{{id=146| /// }}}