Thank you to the organizers for the invitation to come and introduce you all to Sage.
Let's get started. Hopefully by now you and I don't need help with the following computations:
{{{id=0| 2+2 /// 4 }}} {{{id=2| integrate(x^3,x) /// 1/4*x^4 }}} {{{id=48| 3*vector([1,2,3]) /// (3, 6, 9) }}} On the other hand, if I needed one of the following as an intermediate step, I might appreciate some assistance: {{{id=33| factorial(100) /// 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 }}} {{{id=49| matrix([[1,2,3],[4,5,7],[6,8,9]])^-1 /// [-11/7 6/7 -1/7] [ 6/7 -9/7 5/7] [ 2/7 4/7 -3/7] }}}What if someone sends me a possible counterexample to FLT?
{{{id=57| 27784^3 + 35385^3 - 40362^3 /// 1 }}}Good thing we checked. Now let's do an integral.
{{{id=9| integral(3*x^2/sqrt(25*x^2-3), x) /// 3/50*sqrt(25*x^2 - 3)*x + 9/250*log(50*x + 10*sqrt(25*x^2 - 3)) }}} Hmm, the antiderivative looks kind of ugly. Let me try to simplify it and typeset it nicely. {{{id=7| show(expand(integral(3*x^2/sqrt(25*x^2-3),x))) /// }}} That's better. Well, maybe I should just type it in my notes, so I don't forget. $$\frac{3}{50}\sqrt{25x^2-3x}+\frac{9}{250}\log\left(50x+10\sqrt{25x^2-3}\right)$$ Of course, maybe I needed a definite integral. {{{id=51| show(integral(3*x^2/sqrt(25*x^2-3),x,1,2)) /// }}} And maybe I needed it numerically approximated, with a built-in error bound computed. {{{id=50| numerical_integral(3*x^2/sqrt(25*x^2-3),1,2) /// (0.9262503194124578, 1.0283444311781162e-14) }}} And maybe I need to visualize that as well... {{{id=53| show(plot(3*x^2/sqrt(25*x^2-3),(x,1,2),fill=True)+plot(3*x^2/sqrt(25*x^2-3),(x,.5,3)),figsize=3) html("This shows $\int_1^2 \\frac{3x^2}{\sqrt{25x^2-3}}\; dx$") ///
This shows
}}}
Even if you haven't done these sorts of things yourself, you know you can do so with computer help, and possibly have something in a lab or on your laptop which can do them - or you use some web app or other.
But creating something like this takes just a bit more time.
{{{id=35|
%hide
def Newton_Iterates(a_function, guess, iterations = 5):
'''
Returns a list of the iterations of the Newton-Raphson method approximations to a zero of a_function.
INPUT:
a_function: a function of one variable
start: the starting value of the iteration
iterations: (optional) the number of iterations to draw
'''
f(x)=a_function(x)
deriv=f.derivative()
approx(x)=x-f(x)/deriv(x)
input = guess
iter_list = [input]
for i in range(iterations):
input = approx(input).n()
iter_list.append(input)
return iter_list
def Newton_Graph(a_function, guess, iterations = 5):
'''
Returns a graphics object of a plot of a_function and Newton-Raphson method tangent lines.
INPUT:
a_function: a function of one variable
start: the starting value of the iteration
iterations: (optional) the number of iterations to draw
xmin: (optional) the lower end of the plotted interval
xmax: (optional) the upper end of the plotted interval
EXAMPLES:
sage: f = lambda x: x^3-3*x^2+2*x
sage: show(Newton_Graph(f,.5))
Note: This is very slow with symbolic functions.
'''
iter_list = [[guess,0]]
f(x)=a_function(x)
deriv=f.derivative()
approx(x)=x-f(x)/deriv(x)
input = guess
for i in range(iterations):
iter_list.append([input,f(input)])
input=approx(input).n()
iter_list.append([input,0])
approx_tangents = line(iter_list,rgbcolor=(1,0,0))
xmin=min([point[0] for point in iter_list])-1
xmax=max([point[0] for point in iter_list])+1
ymin=min([point[1] for point in iter_list])-.5
ymax=max([point[1] for point in iter_list])+.5
basic_plot = plot(a_function, xmin, xmax, color='blue')
P=basic_plot + approx_tangents
P.show(ymin=ymin,ymax=ymax)
def cubic_approx(guess, intercept=0, iterations = 5):
f(x)=x^3-3*x^2+2*x+intercept
return Newton_Graph(f,guess,iterations=iterations)
def cubic_iterate(guess, intercept=0, iterations = 5):
f(x)=x^3-3*x^2+2*x+intercept
return Newton_Iterates(f,guess,iterations=iterations)
@interact
def _(guess=RR(.5),intercept=slider(-2,2,.5,0),iterations=slider(1,20,1,5)):
root = cubic_iterate(guess,intercept=intercept,iterations=iterations)[-1]
print 'start guess=', guess
print 'intercept=', intercept
print 'iterations =', iterations
print 'root approximation =', root
cubic_approx(guess,intercept=intercept,iterations=iterations)
///
|
|||||||||||||||||||||
It's so great to have this computational and visualization power. Unfortunately, with some such software, you may have run into one of the following three barriers to using it more widely.
Sage, the software I am currently using in this talk, is a free, powerful option which avoids all these issues. You can download it freely, OR you can run it off a server for interface in your web browser- just as I am doing right now! Indeed, there are four interfaces:
Let me demonstrate each of these for you briefly.
(Demo Interlude)
However, we'll stick with the notebook here.
Sage handles other topics in undergraduate mathematics with aplomb, such as vector calculus, number theory, applied linear algebra, differential equations, statistics...
{{{id=36| x,y,z=var('x y z') plot_vector_field3d((x*cos(z),-y*cos(z),sin(z)), (x,0,pi), (y,0,pi), (z,0,pi),colors=['red','green','blue']) /// }}} This one is best viewed using 3D glasses, of course. The following one shows many theorems about numbers, color-encoded! {{{id=37| @interact def power_table_plot(p=(7,prime_range(50))): P=matrix_plot(matrix(p-1,[mod(a,p)^b for a in range(1,p) for b in srange(p)]),cmap='jet') show(P) ///
|
|||||||||||||
If you like numerical methods, these come built-in as well. Note the interactivity Sage can provide.
{{{id=56| %hide import pylab A_image = pylab.mean(pylab.imread(DATA + 'san_thome.png'), 2) B_image = pylab.mean(pylab.imread(DATA + 'imsc_night.png'), 2) @interact def svd_image(i=(5,(1..50))): u,s,v = pylab.linalg.svd(A_image) A = sum(s[j]*pylab.outer(u[0:,j],v[j,0:]) for j in range(i)) u1,s1,v1 = pylab.linalg.svd(B_image) B = sum(s1[j]*pylab.outer(u1[0:,j],v1[j,0:]) for j in range(i)) g = graphics_array([[matrix_plot(A),matrix_plot(B)],[matrix_plot(A_image),matrix_plot(B_image)]]) show(g,axes=False, figsize=6) html('
|
|||||||||||||||
|
|||||||||||||||||||
Statistics are possible via some very basic native Sage functionality, or via several very powerful components of Sage that are primarily statistical/numerical. By far the most powerful of these is R.
{{{id=44| %r library("MASS") data(Cars93) mean(Cars93$MPG.city) xtabs( ~ Origin + MPG.city, data=Cars93) /// [1] 22.36559 MPG.city Origin 15 16 17 18 19 20 21 22 23 24 25 26 28 29 30 31 32 33 39 42 46 USA 2 3 4 5 8 3 3 4 7 2 2 0 1 2 0 2 0 0 0 0 0 non-USA 0 0 4 7 2 5 3 3 1 3 4 2 1 4 1 0 1 1 1 1 1 }}}Of course, one reason I like Sage is not all these 'useful' things, but producing eye candy:
{{{id=47| plot([x^i for i in [1..12]],(x,0,1),fill = dict((i,[i+1]) for i in [0..11])) ///
}}}
Note that this means any computation you want to do that Sage can do with graphs, you can in principle do by just "drawing" the graph.
Sage consists of many separate components, each of which is an open-source standard in its field. Let's see a nice example of all these things working together.
{{{id=31| H = graphs.TetrahedralGraph() # Takes NetworkX graph /// }}} {{{id=38| H = H.cartesian_product(H) # Uses Sage functionality /// }}} {{{id=40| show(H) # uses matplotlib ///
}}}
{{{id=24|
Sym = H.automorphism_group() # Native Sage Code using C backend for graphs
///
}}}
{{{id=25|
len( Sym.list() ) # Needs GAP to calculate
///
1152
}}}
{{{id=26|
Sym._gap_().ComputedPCentralSeriess() # Directly within GAP
///
[ ]
}}}
{{{id=60|
poly = H.chromatic_polynomial() # Uses compiled C code from Cython in Sage
///
}}}
{{{id=61|
deriv = derivative(poly,x) # Uses Ginac/Pynac for symbolics
///
}}}
{{{id=63|
deriv.subs(x=0)
///
-188580960
}}}
And yes, this graph is connected.
{{{id=65| integral(poly,x) # Uses Maxima for integration, symbolic summation, etc. /// 1/17*x^17 - 3*x^16 + 1096/15*x^15 - 7902/7*x^14 + 161032/13*x^13 - 102358*x^12 + 7242454/11*x^11 - 16789056/5*x^10 + 123182360/9*x^9 - 89140779/2*x^8 + 114946568*x^7 - 230543668*x^6 + 1742815741/5*x^5 - 376891986*x^4 + 264255876*x^3 - 94290480*x^2 }}}The book "Quo Vadis, Graph Theory?: A Source Book for Challenges and Directions" asks for interpretations of such things.
{{{id=30| /// }}}There are many other great aspects of Sage I haven't mentioned yet. Let's see how each of these works.
But the best thing to do is to try it out or download it!
Here is a gallery of other research-level things. But to really explore, one needs to peruse the gargantuan reference manual.
I like this example, as it shows that Riemann's exact formula for $\pi(x)$, the prime counting function, is completely different in nature than other approximations.
{{{id=29| import mpmath var('y') L = lcalc.zeros_in_interval(10,50,0.1) @interact def _(n=(100,(60,10^3))): P = plot(prime_pi,n-50,n,color='black',legend_label='$\pi(x)$') P += plot(Li,n-50,n,color='green',legend_label='$Li(x)$') G = lambda x: sum([mpmath.li(x^(1/j))*moebius(j)/j for j in [1..3]]) P += plot(G,n-50,n,color='red',legend_label='$\sum_{j=1}^{%s}\\frac{\mu(j)}{j}Li(x^{1/j})$'%3) F = lambda x: sum([(mpmath.li(x^(1/j))-log(2)+numerical_integral(1/(y*(y^2-1)*log(y)),x^(1/j),oo)[0])*moebius(j)/j for j in [1..3]])-sum([(mpmath.ei(log(x)*((0.5+l[0]*i)/j))+mpmath.ei(log(x)*((0.5-l[0]*i)/j))).real for l in L for j in [1..3]]) P += plot(F,n-50,n,color='blue',legend_label='Really good estimate',plot_points=50) show(P) ///
|
|||||||||||||||
Here's a result from the last couple years connected to work the founder of Sage, William Stein, is doing, as well as to Ramanujan (some of you will know the connection).
Theorem:
Let $r_{12}(n)$ denote the number of ways to write $n$ as a sum of twelve squares (where order and sign both matter, so $(1,2)$ and $(2,1)$ and $(-2,1)$ are all different). As we let $p$ go through the set of all prime numbers, the distribution of the fraction $$\frac{r_{12}(p)-8(p^5+1)}{32p^{5/2}}$$ is asymptotic to $$\frac{2}{\pi}\sqrt{1-t^2}\; .$$
{{{id=66| def dist(v, b, left=float(0), right=float(pi)): """ We divide the interval between left (default: 0) and right (default: pi) up into b bins. For each number in v (which must left and right), we find which bin it lies in and add this to a counter. This function then returns the bins and the number of elements of v that lie in each one. ALGORITHM: To find the index of the bin that a given number x lies in, we multiply x by b/length and take the floor. """ length = right - left normalize = float(b/length) vals = {} d = dict([(i,0) for i in range(b)]) for x in v: n = int(normalize*(float(x)-left)) d[n] += 1 return d, len(v) def graph(d, b, num=5000, left=float(0), right=float(pi)): s = Graphics() left = float(left); right = float(right) length = right - left w = length/b k = 0 for i, n in d.iteritems(): k += n # ith bin has n objects in it. s += polygon([(w*i+left,0), (w*(i+1)+left,0), \ (w*(i+1)+left, n/(num*w)), (w*i+left, n/(num*w))],\ rgbcolor=(0,0,0.5)) return s def sqrt2(): PI = float(pi) return plot(lambda x: (2/PI)*math.sqrt(1-x^2), -1,1, plot_points=200, rgbcolor=(0.3,0.1,0.1), thickness=2) delta = delta_qexp(10^5) @interact def delta_dist(b=(20,(10..150)), number = (500,1000,..,delta.prec())): D = delta[:number] w = [float(D[p])/(2*float(p)^(5.5)) for p in prime_range(number + 1)] d, total_number_of_points = dist(w,b,float(-1),float(1)) show(graph(d, b, total_number_of_points,-1,1) + sqrt2(), frame=True, gridlines=True) ///
|
|||||||||||||||||||
The Sage-combinat crew has a huge number of good tutorials that are not (yet) in the official Sage documentation. While we're on the subject of prime numbers, here is one which shows how Cython can be used to speed up calculations of interest.
{{{id=67| sage: def first_primes_python(m): ... primes_list = [] ... n = 2 ... while len(primes_list) < m: ... n_is_prime = True ... for p in primes_list: ... if n % p == 0: ... n_is_prime = False ... break ... if n_is_prime == True: ... primes_list.append(n) ... n = n + 1 ... return primes_list /// }}} {{{id=69| %cython def first_primes_cython_v1(m): primes_list = [] n = 2 while len(primes_list) < m: n_is_prime = True for p in primes_list: if n % p == 0: n_is_prime = False break if n_is_prime == True: primes_list.append(n) n = n + 1 return primes_list /// }}} {{{id=70| time p = first_primes_python(5000) /// Time: CPU 1.82 s, Wall: 1.82 s }}} {{{id=71| time p = first_primes_cython_v1(5000) /// Time: CPU 0.45 s, Wall: 0.45 s }}}It turns out that one can do quite complicated things combining combinatorics and power series within Sage. This example is too long for me to include directly, so let's go to the link! (A world-readable link to a non-live version is here.)
As it says, "Therefore, instead of solving the equation, we look for the equation describing the object which is best suited to the problem we want to solve."
{{{id=72| /// }}}I don't have a heading for this, but it just shows that you can do a lot in Sage. This actually has applications in social choice theory, usually considered as a subbranch of economics.
{{{id=74| d5 = {12345:[12354,12435,13245,15234,13452],12354:[12534,13254,14235,13542],12435:[15243,12453,14352,14235],13245:[15324,12453,13254,13425],15234:[14523,15243,12534,15324],13452:[13425,13542,14523,14352],15432:[15423,15342,14532,14325,12543],12543:[13254,15243,12534,12453],14325:[14352,14235,13254,13425],14532:[15324,14352,12453,14523],15342:[13425,15324,12534,13542],15423:[15243,14235,14523,13542],13524:[13542,15243,13254,14352,15324],14253:[12534,12453,14523,14235,13425]} G5 = Graph(d5) /// }}} {{{id=75| plot(G5,figsize=4) ///
}}}
{{{id=76|
Gp5 = G5.automorphism_group()
print Gp5.order()
///
480
}}}
{{{id=77|
Gp5.character_table()
///
[ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[ 1 -1 1 1 -1 1 -1 -1 1 1 1 1 -1 -1 1 1 1 1 1]
[ 1 -1 1 1 1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1]
[ 1 1 1 1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 -1 -1 1]
[ 2 0 -2 2 0 0 0 0 0 0 2 -2 0 0 -2 2 0 0 -2]
[ 4 0 0 -1 -2 0 1 0 -1 1 1 1 2 -1 -1 0 1 -4 4]
[ 4 0 0 -1 -2 0 1 0 1 -1 1 1 -2 1 -1 0 -1 4 4]
[ 4 0 0 -1 2 0 -1 0 -1 1 1 1 -2 1 -1 0 1 -4 4]
[ 4 0 0 -1 2 0 -1 0 1 -1 1 1 2 -1 -1 0 -1 4 4]
[ 5 -1 1 0 -1 -1 -1 1 1 0 -1 -1 1 1 0 1 0 -5 5]
[ 5 -1 1 0 1 1 1 -1 -1 0 -1 -1 1 1 0 1 0 5 5]
[ 5 1 1 0 -1 1 -1 1 -1 0 -1 -1 -1 -1 0 1 0 5 5]
[ 5 1 1 0 1 -1 1 -1 1 0 -1 -1 -1 -1 0 1 0 -5 5]
[ 6 0 -2 1 0 -2 0 0 0 1 0 0 0 0 1 -2 1 6 6]
[ 6 0 -2 1 0 2 0 0 0 -1 0 0 0 0 1 -2 -1 -6 6]
[ 6 0 2 1 0 0 0 0 0 2*zeta5^3 + 2*zeta5^2 + 1 0 0 0 0 -1 -2 -2*zeta5^3 - 2*zeta5^2 - 1 0 -6]
[ 6 0 2 1 0 0 0 0 0 -2*zeta5^3 - 2*zeta5^2 - 1 0 0 0 0 -1 -2 2*zeta5^3 + 2*zeta5^2 + 1 0 -6]
[ 8 0 0 -2 0 0 0 0 0 0 2 -2 0 0 2 0 0 0 -8]
[ 10 0 -2 0 0 0 0 0 0 0 -2 2 0 0 0 2 0 0 -10]
}}}
{{{id=78|
Chars5 = Gp5.irreducible_characters()
CCRs5 = Gp5.conjugacy_classes_representatives()
P5 = [sum([len(z) for z in x.cycle_tuples(True) if len(z)==1]) for x in CCRs5]
for ch in Chars5:
mysum = 0
for i in range(19):
mysum += P5[i]*ch.values()[i] * 480 / Gp5.centralizer(CCRs5[i]).order()
print mysum/480
///
1
0
0
1
0
0
0
0
0
0
0
1
1
0
0
1
1
0
0
}}}
To me, this is basically magic. I hope you will enjoy similar magical experiences with Sage!
{{{id=79| /// }}}