PREP Tutorial, Advanced Plotting, Part 1
system:sage


<h1>Sage Tutorial for Advanced Plotting - Part 1</h1>
<p>This&nbsp;<a href="http://www.sagemath.org" target="_blank">Sage</a>&nbsp;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).</p>
<p>Thanks to Sage's integration of projects like <a href="http://matplotlib.sourceforge.net/" target="_blank">matplotlib</a>, <a href="http://jmol.sourceforge.net/" target="_blank">Jmol</a>, and <a href="http://jedi.ks.uiuc.edu/~johns/raytracer/" target="_blank">Tachyon</a>, Sage has comprehensive plotting capabilities.  This pair of worksheets aims to go more in depth with them than previous tutorials. &nbsp; This worksheet consists of the following sections:</p>
<ul>
<li><a href="#Cartesian">Cartesian Plots</a></li>
<li><a href="#Parametric">Parametric Plots</a></li>
<li><a href="#Polar">Polar Plots</a></li>
<li><a href="#Data">Plotting Data Points</a></li>
<li><a href="#Contour">Contour Plots</a></li>
<li><a href="#3D">3D Plots</a></li>
</ul>
<p>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).</p>
<h1 id="Cartesian">Cartesian Plots</h1>
<p>A simple quadratic is easy.</p>

{{{id=3|
plot(x^2, (x,-2,2))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>You can combine "plot objects" by adding them.</p>

{{{id=6|
regular = plot(x^2, (x,-2,2), color= 'purple')
skinny = plot(4*x^2, (x,-2,2), color = 'green')
regular + skinny
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p><strong>Problem</strong>:&nbsp; Plot a green <span style="color: #008000;">$y=\sin(x)$</span> together with a red <span style="color: #ff0000;">$y=2\,\cos(x)$</span>.&nbsp; (Hint: you can use "pi" as part of your range.)</p>

{{{id=12|

///
}}}

<p>Boundaries of a plot can be specified, in addition to the overall size.</p>

{{{id=9|
plot(1+e^(-x^2), xmin=-2, xmax=2, ymin=0, ymax=2.5, figsize=10)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p><strong>Problem</strong>:&nbsp; Plot $y=5+3\,\sin(4x)$ with suitable boundaries.</p>

{{{id=52|

///
}}}

<p>You can add lots of extra information.</p>

{{{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
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>You can fill regions with transparent color, and thicken the curve. &nbsp;This example uses several options to fine-tune our graphic.</p>

{{{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
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

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

{{{id=19|

///
}}}

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

{{{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)
///
}}}

<p><strong>Problem</strong>: These parametric equations will create a hypocycloid,</p>
<p>$$x(t)=17\cos(t)+3\cos(17t/3)$$</p>
<p>$$y(t)=17\sin(t)-3\sin(17t/3)$$</p>
<p>Create this as a parametric plot.</p>

{{{id=34|

///
}}}

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

{{{id=89|
var('t')
parametric_plot((t^2,sin(t)), (t,0,pi))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{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))
///
}}}

<h1 id="Polar">Polar Plots</h1>
<p>Sage can also do polar plots.</p>

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

<p>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.</p>
<p>More than one polar curve can be specified in a list (square brackets).&nbsp; Notice the automatic graded shading of the fill color.</p>

{{{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)
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>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)$.</p>

{{{id=31|

///
}}}

{{{id=55|

///
}}}

{{{id=54|

///
}}}

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

{{{id=47|
%hide
%auto
html('<h2>Sine and unit circle (by Jurgis Pralgauskis)</h2> inspired by <a href="http://www.youtube.com/watch?v=Ohp6Okk_tww&feature=related">this video</a>' )
# 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)
///
<html><font color='black'><h2>Sine and unit circle (by Jurgis Pralgauskis)</h2> inspired by <a href="http://www.youtube.com/watch?v=Ohp6Okk_tww&feature=related">this video</a></font></html>
<html><!--notruncate--><div padding=6 id="div-interact-47"> <table width=800px height=20px bgcolor="#c5c5c5"
                 cellpadding=15><tr><td bgcolor="#f9f9f9" valign=top align=left><table>
<tr><td colspan=3><table><tr><td align=right><font color="black">Pasirinkite kampą:    &nbsp;</font></td><td><table><tr><td>
        <div id="slider-angle-47" style="margin:0px; margin-left: 1.0em; margin-right: 1.0em; width: 15.0em;"></div>
        </td><td><font color="black" id="slider-angle-47-lbl"></font></td></tr></table><script>(function(){ var values = ["0","1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30","31","32","33","34","35","36","37","38","39","40","41","42","43","44","45","46","47","48","49","50","51","52","53","54","55","56","57","58","59","60","61","62","63","64","65","66","67","68","69","70","71","72","73","74","75","76","77","78","79","80","81","82","83","84","85","86","87","88","89","90","91","92","93","94","95","96","97","98","99","100","101","102","103","104","105","106","107","108","109","110","111","112","113","114","115","116","117","118","119","120","121","122","123","124","125","126","127","128","129","130","131","132","133","134","135","136","137","138","139","140","141","142","143","144","145","146","147","148","149","150","151","152","153","154","155","156","157","158","159","160","161","162","163","164","165","166","167","168","169","170","171","172","173","174","175","176","177","178","179","180","181","182","183","184","185","186","187","188","189","190","191","192","193","194","195","196","197","198","199","200","201","202","203","204","205","206","207","208","209","210","211","212","213","214","215","216","217","218","219","220","221","222","223","224","225","226","227","228","229","230","231","232","233","234","235","236","237","238","239","240","241","242","243","244","245","246","247","248","249","250","251","252","253","254","255","256","257","258","259","260","261","262","263","264","265","266","267","268","269","270","271","272","273","274","275","276","277","278","279","280","281","282","283","284","285","286","287","288","289","290","291","292","293","294","295","296","297","298","299","300","301","302","303","304","305","306","307","308","309","310","311","312","313","314","315","316","317","318","319","320","321","322","323","324","325","326","327","328","329","330","331","332","333","334","335","336","337","338","339","340","341","342","343","344","345","346","347","348","349","350","351","352","353","354","355","356","357","358","359","360"]; setTimeout(function() {
    $('#slider-angle-47').slider({
        step: 1, min: 0, max: 360, value: 30,
        change: function (e,ui) { var position = ui.value; if(values!=null) $('#slider-angle-47-lbl').text(values[position]); interact(47, '_interact_.update(\'47\', \'angle\', 1, _interact_.standard_b64decode(\''+encode64(position)+'\'), globals()); _interact_.recompute(\'47\');'); },
        slide: function(e,ui) { if(values!=null) $('#slider-angle-47-lbl').text(values[ui.value]); }
    });
    if(values != null) $('#slider-angle-47-lbl').text(values[$('#slider-angle-47').slider('value')]);
    }, 1); })();</script></td>
</tr></table></td></tr>
<tr><td></td><td style='width: 100%;'><div id="cell-interact-47"><?__SAGE__START>
        <table border=0 bgcolor="white" width=100%>
        <tr><td bgcolor="white" align=left valign=top><pre><?__SAGE__TEXT></pre></td></tr>
        <tr><td  align=left valign=top><?__SAGE__HTML></td></tr>
        </table><?__SAGE__END></div></td><td></td></tr>
<tr><td colspan=3></td></tr>
</table></td>
                 </tr></table></div>
                 </html>
}}}

<h1 id="Data">Plotting Data Points</h1>
<p>Sometimes one wishes to simply plot data. &nbsp;Here, we demonstrate several ways of plotting points and data via the simple approximation to the Fibonacci numbers given by</p>
<p>$$F_n=\frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n\; ,$$</p>
<p>which is quite good after about $n=5$.</p>
<p>First, we notice that the Fibonacci numbers are built in.</p>

{{{id=61|
fibonacci_sequence(6)
///
<generator object fibonacci_sequence at 0x634c460>
}}}

{{{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)]
}}}

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

{{{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))]
}}}

<p>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.</p>

{{{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|

///
}}}

<h1 id="Contour">Contour Plots</h1>
<p>Contour plotting can be very useful when trying to get a handle on multivariable functions, as well as modeling. &nbsp;The basic syntax is essentially the same as for 3D plotting - simply an extension of the 2D plotting syntax.</p>

{{{id=118|
f(x,y)=y^2+1-x^3-x
contour_plot(f, (x,-pi,pi), (y,-pi,pi))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

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

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

<p>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). &nbsp;The strange-looking construction corresponding to 'label_fmt' is a Sage/Python data type called a <em>dictionary</em>, and turns out to be useful for more advanced Sage use;&nbsp;it consists of pairs connected by a colon, all inside curly braces.</p>

{{{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')
///
}}}

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

{{{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|

///
}}}

<p>A density plot is like a contour plot, but without discrete levels.</p>
</div>

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

<p>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.</p>

{{{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')
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

<p>Note: mention something about the options</p>

{{{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))
///
<html><font color='black'><img src='cell://sage0.png'></font></html>
}}}

{{{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|

///
}}}

<p>3d contour plots are given by the implicit_plot3d function.</p>

{{{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)
///
}}}

<p><strong style="font-weight: bold;">Problem</strong>:&nbsp;</p>
<ol>
<li>Plot the surface $z=f(x,y)=x^2y$ on a reasonable domain.</li>
<li>In a contrasting color, plot the tangent plane to this surface at $x=1$, $y=3$.</li>
<li>Adjust the opacity of each to make viewing easier.</li>
<li>Place a "point" graphics object in a third color at the location of the tangency.</li>
</ol></div>

{{{id=125|

///
}}}

<h1 id="3D">3D Plotting</h1>
<p>Plotting surfaces works as you might expect.&nbsp; The JMOL Java applet has lots of built-in features:&nbsp; zooming, rotating, stereoscopic viewing, automatic spin.&nbsp; Note the need to specify the variables.</p>

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

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

{{{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])
///
}}}

<h2>Plotting with transformations</h2>
<p>You can easily apply a transformation to your 3d function. &nbsp;Built-in transformations include spherical and cylindrical transformations.</p>

{{{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))
///
}}}

<p>You can make your own transformations and apply them easily. &nbsp;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$).</p>

{{{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|

///
}}}