MATH 8650 Intro to Sage

2659 days ago by luogend

MTHS 8650

Introduction to Sage

Daniel D. Warner

WWW.SAGEMATH.ORG
 
       

Mission Statement

Create a viable free open source alternative to Magma, Maple, Mathematica, and Matlab.


History of Sage

 

    • William Stein started Sage at Harvard in January 2005.

    • Sage-1.0 released February 2006 at Sage Days 1 (UC San Diego).

    • Over 50 Sage Days Workshops (!)  at UCLA, UW, Cambridge, Bristol, Austin, France, San Diego, Seattle, MSRI, ..., Barcelona, ...

    • Sage won first prize in the Trophees du Libre (November 2007)
  • Funding from Microsoft, Univ of Washington, UC San Diego, NSF, DoD, Google, Sun, private donations, etc.

 

  • Hundreds of other people subsequently got involved in Sage's development, and the scope of the project has widened to cover all mathematical computation. There is interest from all areas of mathematics, physics, engineering, etc.  (Developer mailing list has 1184 subscribers and about 50 messages/day... and sage-support has 1718 subscribers.)

 

What is Sage?

  • Sage = Python + Math
  • A unified self-contained distribution of open source free mathematical software.
  • Nearly a half million lines of new Python (and Cython) code that implements new capabilities and algorithms.
  • A "cloud" application like GMail or Google Docs: http://sagenb.org  (currently about 30,000 users);  Of course Sage also runs on your desktop and supercomputer.


 

Summary:
Sage is about building the car instead of reinventing the wheel

  1. Sage uses Python,  a mainstream programming language, instead of inventing a custom mathematics language
  2. Use straightforward method to link programs together -- C library and pseudotty's, instead of XML servers/OpenMath.  We implement all conversion routines, instead of expecting upstream to do it: we make them communicate with Sage, whether they want to or not.
  3. Give copious credit to contributors and be very developer friendly (easily build from source).
  4. Reuse, improve, and contribute to existing libraries and projects (e.g., Scipy, Numpy, R, ATLAS, CVXopt, GSL), instead of starting over and competing with them.
  5. Make the GUI using a web browser: the world of java and javascript is available and Sage integrates with the web.

 
 
       

Key take-away points: 

  • FOSS: 100% free open source software; good for clusters, sharing, research, collaboration
  • Python: very mainstream, scientific-computing friendly programming language
  • Cython: optimized Python to C/C++ compiler we develop with support for C/C++ datatypes
  • Interfaces: to control Matlab, Octave, Mathematica, etc., from Sage
  • Self contained: can have many copies of Sage at once, so no "dependency hell"
  • Peer reviewed: get your code refereed
  • Web based: notebook interface
 
       

All aboard!

We want to pause and make sure that everyone can access Sage.

  1. Access the Sage Server at  https://sage.math.clemson.edu:34567
  2. Set up an account and click on the "Published Button."  All worksheets for this course will start with MTHS 8650.
  3. You can easily download a copy of Sage from WWW.SAGEMATH.ORG, and install it on your personal computer.

 
 
       

And now a short demo...
# This is a simple introduction to Sage. # # The Sage notebook uses blocks of text, which are essentially Python statements that get executed. # Lines that start with the pound sign are comments and do not get echoed when the block is evaluated. # # A block is evaluated by clicking on the evaluate link or by typing shift-enter. # # Pressing tab will give auto-complete options, or function descriptions after a parens ( # A '?' after a command gives detailed information, e.g. 'plot?' tells you about plot # Sage can do simple calculations 3+2*6^2 
        
75
75
# and more elaborate calculations factorial(100) 
        
933262154439441526816992388562667004907159682643816214685929638952175999\
932299156089414639761565182862536979208272237582511852109168640000000000\
00000000000000
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
# or 1000 digits of the square root of pi sqpi = sqrt(pi).n(digits=1000) sqpi 
        
1.7724538509055160272981674833411451827975494561223871282138077898529112\
845910321813749506567385446654162268236242825706662361528657244226025250\
937096027870684620376986531051228499251730289508262289320953792679628001\
746390153514797205167001901852340185854469744949126403139217755259062164\
054193325009063984076137334774751534336679897893658518364087954511651617\
387600590673934317913328098548462481849020546548521956132515616474675150\
427387610561079961271072100603720444836723652966137080943234988316684242\
138457096091204204277857780686947665700052183056851254133966369446541815\
107166938833219429293570622688652244205421499480499207564863988748385059\
306402182140292858112330649789452036211490789622873894032459781985131348\
712665125062932600446563821096750268124969305954204615607619522173915250\
702077927580990543329006622230676144696612481887430699788352050614644438\
541853079735742571791856359597499599522638492422038891039664064472939728\
41345043002140564233433039261756134176336320017037654163476320669
1.772453850905516027298167483341145182797549456122387128213807789852911284591032181374950656738544665416226823624282570666236152865724422602525093709602787068462037698653105122849925173028950826228932095379267962800174639015351479720516700190185234018585446974494912640313921775525906216405419332500906398407613733477475153433667989789365851836408795451165161738760059067393431791332809854846248184902054654852195613251561647467515042738761056107996127107210060372044483672365296613708094323498831668424213845709609120420427785778068694766570005218305685125413396636944654181510716693883321942929357062268865224420542149948049920756486398874838505930640218214029285811233064978945203621149078962287389403245978198513134871266512506293260044656382109675026812496930595420461560761952217391525070207792758099054332900662223067614469661248188743069978835205061464443854185307973574257179185635959749959952263849242203889103966406447293972841345043002140564233433039261756134176336320017037654163476320669
sqpi^2 
        
3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815209209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623799627495673518857527248912279381830119491298\
336733624406566430860213949463952247371907021798609437027705392171762931\
767523846748184676694051320005681271452635608277857713427577896091736371\
787214684409012249534301465495853710507922796892589235420199561121290219\
608640344181598136297747713099605187072113499999983729780499510597317328\
160963185950244594553469083026425223082533446850352619311881710100031378\
387528865875332083814206171776691473035982534904287554687311595628638823\
53787593751957781857780532171226806613001927876611195909216420199
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
pi.n(digits=1000) 
        
3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446095505822317253\
594081284811174502841027019385211055596446229489549303819644288109756659\
334461284756482337867831652712019091456485669234603486104543266482133936\
072602491412737245870066063155881748815209209628292540917153643678925903\
600113305305488204665213841469519415116094330572703657595919530921861173\
819326117931051185480744623799627495673518857527248912279381830119491298\
336733624406566430860213949463952247371907021798609437027705392171762931\
767523846748184676694051320005681271452635608277857713427577896091736371\
787214684409012249534301465495853710507922796892589235420199561121290219\
608640344181598136297747713099605187072113499999983729780499510597317328\
160963185950244594553469083026425223082533446850352619311881710100031378\
387528865875332083814206171776691473035982534904287554687311595628638823\
53787593751957781857780532171226806613001927876611195909216420199
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
# Sage can factor integers factor(798) 
        
2 * 3 * 7 * 19
2 * 3 * 7 * 19
# The same operations can be applied to algebraic expressions p = x^5-1 p 
        
x^5 - 1
x^5 - 1
factor(p) 
        
(x^4 + x^3 + x^2 + x + 1)*(x - 1)
(x^4 + x^3 + x^2 + x + 1)*(x - 1)
p^3; factor(p^3); 
        
(x^5 - 1)^3
(x^4 + x^3 + x^2 + x + 1)^3*(x - 1)^3
(x^5 - 1)^3
(x^4 + x^3 + x^2 + x + 1)^3*(x - 1)^3
# Algebraic expressions can be plotted P1 = plot( p^2 , -0.5, 0.5, rgbcolor=(1,0,0) ) P2 = plot( p^3, -0.5, 0.5, rgbcolor=(0,1,0) ) show(P1+P2) 
       
# Sage handles rational functions with ease rf = x / (x^2 - 1) rf 
        
x/(x^2 - 1)
x/(x^2 - 1)
P3 = plot( rf , -0.75, 0.75, rgbcolor=(1,0,0) ) show(P3) 
       
# and partial fraction expansions pf = rf.partial_fraction(x) pf 
        
1/2/(x + 1) + 1/2/(x - 1)
1/2/(x + 1) + 1/2/(x - 1)
# Sage can display results in the more familiar mathematical format show(pf) 
       



# Sage uses LaTeX for the mathematical formatting latex(pf) 
        
\frac{1}{2 \, {\left(x + 1\right)}} + \frac{1}{2 \, {\left(x -
1\right)}}
\frac{1}{2 \, {\left(x + 1\right)}} + \frac{1}{2 \, {\left(x - 1\right)}}

Solving Algebraic Equations
# A quadratic example eq1 = x^2 - 2 == x show(eq1) sol = solve([eq1],x) sol 
        

[x == 2, x == -1]
[x == 2, x == -1]
# Cubic polynomials are just as easy eq2 = x^3 + 2*x^2 - 5*x - 6 == 0 show(eq2) sol = solve([eq2],x) sol 
        

[x == -3, x == -1, x == 2]
[x == -3, x == -1, x == 2]
# Sage can even find the ugly answers eq3 = x^3 - 2 == x show(eq3) sol = solve([eq3],x) for s in sol: show(s) 
       
# A linear systems of equations var('x,y') eq4 = x-1 == y+1 eq5 = x+1 == 2*(y-1) eq4;eq5 
        
x - 1 == y + 1
x + 1 == 2*y - 2
x - 1 == y + 1
x + 1 == 2*y - 2
sol = solve([eq4,eq5],x,y) sol 
        
[[x == 7, y == 5]]
[[x == 7, y == 5]]
sol = solve([eq4,eq5],x,y,solution_dict=True) sol 
        
[{x: 7, y: 5}]
[{x: 7, y: 5}]
[[s[x].n(digits=2),s[y].n(digits=5)] for s in sol] 
        
[[7.0, 5.0000]]
[[7.0, 5.0000]]
# A nonlinear system of equations var('x y') eq6 = (x-1)^2 + (y-1)^2 == 4; show(eq6) eq7 = 4*y - 3*x == -1; show(eq7) sol3 = solve([eq6,eq7],x,y,solution_dict=True) show(sol3) 
       
# Use the dictionary result and display each solution to a different precision for s in sol3: print(s[x].n(digits=5),s[y].n(digits=25)) 
        
(-0.32767, -0.4957550765359254871346963)
(2.8077, 1.855755076535925487134696)
(-0.32767, -0.4957550765359254871346963)
(2.8077, 1.855755076535925487134696)
# Display the nonlinear example and the solutions x0 = sol3[0][x] y0 = sol3[0][y] s0 = '(%5.2f,%5.2f)' % (x0,y0) p0 = (x0-0.47,y0) x1 = sol3[1][x] y1 = sol3[1][y] s1 = '(%5.2f,%5.2f)' % (x1,y1) p1 = (x1-0.47,y1) g = Graphics() g += circle((1,1),2) g += line([(-1,-1),(3,2)]) g += points([(s[x],s[y]) for s in sol3],rgbcolor=(1,0,0)) g += text(s0,p0) g += text(s1,p1) show(g,aspect_ratio=1) 
       
# The following example by Jason Grout uses Sage to solve a system of three non-linear equations. # The fourth variable, u, is treated as a parameter # First, we solve the system symbolically: var('x y u v') eq8 = u+v==9; show(eq8) eq9 = v*y+u*x==-6; show(eq9) eq10 = v*y^2+u*x^2==24; show(eq10) show(solve([eq8,eq9,eq10,u==1],u,v,x,y)) 
       
# Second, we can solve it numerically solns = solve([eq8,eq9,eq10,u==1],u,v,x,y, solution_dict=True) for s in solns: print(s[u].n(digits=10), s[v].n(digits=10), s[x].n(digits=10), s[y].n(digits=10)) 
        
(1.000000000, 8.000000000, -4.883036880, -0.1396203900)
(1.000000000, 8.000000000, 3.549703547, -1.193712943)
(1.000000000, 8.000000000, -4.883036880, -0.1396203900)
(1.000000000, 8.000000000, 3.549703547, -1.193712943)

Matrix Calculations

Try a random 3 by 3 matrix with integer entries.  Then calculate its determinant.
a = matrix.random(ZZ,3,3) show(a) print type(a) time a.determinant() 
        

<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
39
Time: CPU 0.00 s, Wall: 0.00 s
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
39
Time: CPU 0.00 s, Wall: 0.00 s

Try a random 200 by 200 matrix with integer entries.  Then time the calculation of its determinant.  It could be a big number.
a = matrix.random(ZZ,200,200) print a[0,0] print a[199,199] print type(a) time a.determinant() 
        
1
3
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
-14382230777275826589566212167887142284022202947000473201750073971367920\
653260467820935889968458550212479781015365015156203347597476573364560383\
574201803886840423741839879833901820985502312201902589867086430411262481\
120475209803413262175401159033540547216117117017782120211302785040379159\
124538202527695693944994729888152435993127745476801520782655482046465499\
178358461181604148574107498068001571168294012935483809321263021696306732\
000731690393250585869217175998640
Time: CPU 0.54 s, Wall: 0.54 s
1
3
<type 'sage.matrix.matrix_integer_dense.Matrix_integer_dense'>
-14382230777275826589566212167887142284022202947000473201750073971367920653260467820935889968458550212479781015365015156203347597476573364560383574201803886840423741839879833901820985502312201902589867086430411262481120475209803413262175401159033540547216117117017782120211302785040379159124538202527695693944994729888152435993127745476801520782655482046465499178358461181604148574107498068001571168294012935483809321263021696306732000731690393250585869217175998640
Time: CPU 0.54 s, Wall: 0.54 s
 
       

Calculus

Sage does Calculus:
p = (x^5-1) diff(p^3,x) 
        
15*(x^5 - 1)^2*x^4
15*(x^5 - 1)^2*x^4
integral(p^3,x) 
        
1/16*x^16 - 3/11*x^11 + 1/2*x^6 - x
1/16*x^16 - 3/11*x^11 + 1/2*x^6 - x
# Slightly more complicated formula with trigonometric functions with parameters var('a b c') q = (x-a)*(cos(b*x)+sin(c*x)) show(q) 
       
# Differentiate # Illustrates the product rule, the chain rule, and the derivatives of sin and cos. show(diff(q,x)) 
       
# Integrate using integration by parts r = integral(q,x) show(r) 
       
show(expand(r)) 
       
# Define and work with functions f(x) = sin(x)^2*cos(x)*exp(x) show(f) z=4 print f(z) plot(f,0,3) 
        

cos(4)*e^4*sin(4)^2

cos(4)*e^4*sin(4)^2
g = integral(f, x) show(g) 
       
 
       

Illustrating Partial Derivatives
var('x y u v') 
        
(x, y, u, v)
(x, y, u, v)
plot3d((x^2 - y^2) / (x^2 + y^2), (x, -1, 1), (y, -1, 1)) 
       
f(x, y) = x^3 + x^2 * y^3 - 2 * y^2 P = plot3d(f(x, y), (x, 1, 3), (y, 0, 2), opacity = 0.2) P += parametric_plot3d([2 + u, 1, f(2, 1) + diff(f(x, y), x)(x = 2, y = 1) * u], (u, -1/2, 1/2), color = "red", thickness = 2) P += parametric_plot3d([2, 1 + u, f(2, 1) + diff(f(x, y), y)(x = 2, y = 1) * u], (u, -1/2, 1/2), color = "green", thickness = 2) P += point((2, 1, f(2, 1)), color = "black", size = 10) show(P) 
       
f(x, y) = 4 - x^2 - 2 * y^2 P = parametric_plot3d([u * cos(v), 1 / sqrt(2) * u * sin(v), 4 - u^2], (u, 0, 2), (v, 0, 2 * pi), aspect_ratio = [1, 1, 1], opacity = 0.2) P += parametric_plot3d([1 + u, 1, f(1, 1) + diff(f(x, y), x)(x = 1, y = 1) * u], (u, -1/2, 1/2), color = "red", thickness = 2) P += parametric_plot3d([1, 1 + u, f(1, 1) + diff(f(x, y), y)(x = 1, y = 1) * u], (u, -1/2, 1/2), color = "green", thickness = 2) P += point((1, 1, f(1, 1)), color = "black", size = 10) show(P) 
       

3-Dimensional Plots

Sage can draw 3d plots:
# Surface plot example var('x,y') plot3d(cos(x^2 + y^2), (x, -2, 2), (y, -2, 2)) 
       
# Two intersecting surfaces var('x,y') b = 2.2 (plot3d(sin(x^2-y^2),(x,-b,b),(y,-b,b), opacity=.9) + plot3d(0, (x,-b,b), (y,-b,b), color='red')) 
       
# Implicit 3D plot var('x,y,z') T = golden_ratio p = 2 - (cos(x + T*y) + cos(x - T*y) + cos(y + T*z) + cos(y - T*z) + cos(z - T*x) + cos(z + T*x)) r = 4.78 implicit_plot3d(p, (x, -r, r), (y, -r, r), (z, -r, r), plot_points=50) 
       
var('t') 
        
t
t
parametric_plot3d([t^3, ln(3 - t), sqrt(t)], (t, 0.01, 2.99)) 
       
parametric_plot3d([cos(t), sin(t), t], (t, 0, 6 * pi)) 
       
var('s') parametric_plot3d([cos(t), sin(t), s], (t, 0, 2 * pi), (s, 0, 6 * pi)) \ + parametric_plot3d([cos(t), sin(t), t], (t, 0, 6 * pi), color = 'red') 
       
parametric_plot3d([cos(t), sin(t), s], (t, 0, 2 * pi), (s, 0, 4)) \ + parametric_plot3d([s, t, 2 - t], (s, -2, 2), (t, -2, 2), color = 'red') \ + parametric_plot3d([cos(t), sin(t), 2 - sin(t)], (t, 0, 2 * pi), color = 'green') 
       
parametric_plot3d([t^2, 7 * t - 12, t^2], (t, 0, 5)) \ + parametric_plot3d([4 * t - 3, t^2, 5 * t - 6], (t, 0, 5), color = 'red') 
       
 
       

Dynamical Systems Modeling
var('t') y = function('y', t) show(y) type(y) 
        

<type 'sage.symbolic.expression.Expression'>
<type 'sage.symbolic.expression.Expression'>

A Simple Growth Model
var('alpha beta') DE1 = diff(y,t) == alpha*y - beta show(DE1) u1 = desolve(DE1, y, ivar=t) show(u1) 
       
f(t) = expand(u1) show(f) 
       
f(0);f(1);f(1/alpha);f(2) 
        
c + beta/alpha
c*e^alpha + beta/alpha
c*e + beta/alpha
c*e^(2*alpha) + beta/alpha
c + beta/alpha
c*e^alpha + beta/alpha
c*e + beta/alpha
c*e^(2*alpha) + beta/alpha
g(t) = f(c=20,alpha=1/2,beta=40) g 
        
t |--> 20*e^(1/2*t) + 80
t |--> 20*e^(1/2*t) + 80
plot(g,0,5) 
       

A Spring Mass Damper Model
var('t k b') x = function('x',t) DE2 = diff(x,t,t) == -2*x - 1/2*diff(x,t) u2 = desolve(DE2, x, ivar=t) u2 
        
(k2*cos(1/4*sqrt(31)*t) + k1*sin(1/4*sqrt(31)*t))*e^(-1/4*t)
(k2*cos(1/4*sqrt(31)*t) + k1*sin(1/4*sqrt(31)*t))*e^(-1/4*t)
f(t) = u2(k1=1/2,k2=1) f 
        
t |--> 1/2*(2*cos(1/4*sqrt(31)*t) + sin(1/4*sqrt(31)*t))*e^(-1/4*t)
t |--> 1/2*(2*cos(1/4*sqrt(31)*t) + sin(1/4*sqrt(31)*t))*e^(-1/4*t)
plot(f,0,10) 
       

Rumor Model
var('t alpha P') r = function('r',t) model = diff(r,t) == alpha*r*(P-r) u3 = desolve(model,r,ivar=t) u3 
        
-(log(-P + r(t)) - log(r(t)))/(P*alpha) == c + t
-(log(-P + r(t)) - log(r(t)))/(P*alpha) == c + t
var('K') eq = (r-P)/r == K*exp(-alpha*P*t) eq 
        
-(P - r(t))/r(t) == K*e^(-P*alpha*t)
-(P - r(t))/r(t) == K*e^(-P*alpha*t)
solve(eq,r) 
        
[r(t) == -P*e^(P*alpha*t)/(K - e^(P*alpha*t))]
[r(t) == -P*e^(P*alpha*t)/(K - e^(P*alpha*t))]
r=function('r',t) rp = desolve_rk4(diff(r,t) == (1/100)*r*(100-r), r, ics=[0,1],step=0.05, end_points=10) 
       
list_plot(rp) 
       
 
       

Graph Theory

Sage can do graph theory:
g = graphs.FlowerSnark(); g 
        
Flower Snark: Graph on 20 vertices
Flower Snark: Graph on 20 vertices
g 
        
Flower Snark: Graph on 20 vertices
Flower Snark: Graph on 20 vertices

Sage contains many unique and deep algorithms:
g.automorphism_group() 
        
Permutation Group with generators
[(2,11)(3,12)(4,13)(5,8)(6,9)(7,10)(16,19)(17,18),
(1,14)(2,4)(5,7)(8,10)(11,13),
(0,3,6,9,12)(1,4,7,10,13)(2,5,8,11,14)(15,16,17,18,19)]
Permutation Group with generators [(2,11)(3,12)(4,13)(5,8)(6,9)(7,10)(16,19)(17,18), (1,14)(2,4)(5,7)(8,10)(11,13), (0,3,6,9,12)(1,4,7,10,13)(2,5,8,11,14)(15,16,17,18,19)]