CoCalc -- Ewha_Sage_Demo_March

Ewha Sage Demo, 3/18

Path: Ewha_Sage_Demo_March_2018.sagews

Views: ¹³⁴

A Walk Among the SageMath
Travis Scrimshaw (The University of Queensland)
Ewha Womans University
March 20, 2018

SageMath (http://www.sagemath.org/), or Sage for short, is an open-source mathematical software that is built upon numerous open-source packages and uses Python as a common language. It began in 2005 by William Stein (who now runs cocalc, https://wstein.org) as a way to build a better car rather than reinvent the wheel by building on these previously disjoint packages.

In this demo, we will explore some of its features ranging from being a simple calculator to how to use Sage in research.

To install an old version of Sage (7.4) on Windows, you can use the .exe from https://github.com/embray/sage-windows/releases/. You can also build it from source using Cygwin64.

Like every good Sage tutorial, we begin by verifying Sage can do basic arithmetic operations.

2+2

5 - 2

3**2

5 % 3

4/3 + 3/4

If this was plain Python, this would be considered int (floor) division, and would return $1 + 0 = 1$ . However, these are a special Sage integer, so we get rational numbers and can do extra operations on integers.

256.factorial()

12365123412.factor()

6^2

x = 5
x + 1/2 * x - x^4

In order to see what operations are available on integers (or really any object), all one needs to do is x.<tab>, where <tab> means hitting the tab character on your keyboard. This is called tab-completion (it can also fill out partially typed commands).

Don't forget the dot .!

x.binomial(2)

However, it is not always clear what some of these methods (functions acting on a particular type of object, yes you can think of . that way) do. I do not even know all of them. However, Sage has a very useful way to find out: just add ? to the end and it will display the documentation.

Tab-completion and ? will get you very far in exploring what Sage can do and how to use Sage. If you take anything away from this introduction, remember these two things.

x.is_pseudoprime?

Unfortunatly, we cannot write every real number exactly using a computer as that will take an infinite amount of time and space. $\pi$ anyone? So we need to consider real numbers up to a fixed precision.

RR

RR(2)

exp(RR(2))

R3 = RealField(3)
R3

R3(1.2314124312)

R3(5) / R3(3) - R3(1.8)

RR(5) / RR(3) - RR(1.8)

As we can see with the above example, this can cause a loss of precision (remember significant figures?). To get a better sense about the error that can be introduced, we can use interval arithematic.

(RIF, RBF)

RBF3 = RealBallField(3)
RBF3

RBF3(5) / RBF3(3) - RBF3(1.8)

RIF3 = RealIntervalField(3)
RIF3

RIF3(5) / RIF3(3) - RIF3(1.8)

We have the same thing for the complex numbers.

CC, CBF, CIF

ComplexField(4), ComplexBallField(5), ComplexIntervalField(3)

However, if we want to work with known constants, Sage has them.

pi + e

log(2)

log(8,2)

a = (pi+I)^2
a

sin(3)

b = sin(3)^2 + cos(3)^2
b

Since simplifying expressions in general can require a lot of computations, Sage does not always simplify things. Instead, you must tell Sage to do so.

a.expand()

b.simplify_trig()

Subsets of the complex numbers are good, and we can do a lot by just applying basic functions. However, Sage allows you to construct much more complicated functions and do things like calculus.

For the following examples, we will turn Sage's output for CoCalc into latex by default.

%typeset_mode True   # Turn on display output as latex
x,t,z,y = var('x,t,z,y')
t,x,y,z

f = (x^2 - y + cos(z + t^2)) / exp(2+log(x/y + 3))
f

f(x=5, y=pi, z=log(3), t=4)

f(x=5, y=pi, z=log(3), t=4).n()  # numerical approximation

f.derivative(x)

derivative(f, x)

derivative(f, y)

derivative(f, x, t, y)

f.derivative(z,y)

g = f(x=t, y=t, z=-t^2+t)
g

g.simplify_full()

g.integrate(t)

plot(g, [t, -7, 7])

plot(g, [t, -7, 8]) + plot(exp(-2)/4 * (x^2 - x + 1), [x, -7, 8], color='red')

While we are talking about functions and plotting, let's do some other types of plots.

%typeset_mode False     # back to normal text output

parametric_plot((cos(x) + 2 * cos(x/4), sin(x) - 2 * sin(x/4)), [x, 0, 8*pi])

parametric_plot([(3 + cos(3*t))*cos(4*t), (3 + cos(3*t))*sin(4*t), sin(3*t)], [t, -pi, pi], plot_points=200)

parametric_plot( (cos(x) - sin(y), sin(x)^2 + cos(y), cos(x*y)), (x, -pi,pi), (y,-pi,pi), mesh=True)

plot_slope_field(sin(x+y) + cos(x+y), (x,-3,3), (y,-3,3))

p = plot_vector_field((x^2 - y^2*cos(x+y), -3*x/(y^2+1)), (x, -5, 5), (y, -5, 5), aspect_ratio=1)
show(p, figsize=[7,7])

(p + streamline_plot((x^2 - y^2*cos(x+y), -3*x/(y^2+1)), (x, -5, 5), (y, -5, 5), density=1.8)).show(figsize=[7,7])

animate([circle((i,0), i, hue=sin(i/10)) for i in srange(10,0,-1/10)], xmax=20, ymin=-10, ymax=10, aspect_ratio=1).show(delay=10)

k   = 1
f   = (2*x^2-x)*exp(-x^2)
p   = plot(f,0,2*pi, thickness=1, linestyle='--')
c   = 1/pi
order = 15
alpha = [(n, c*numerical_integral(f*sin(x*n/2),0,2*pi)[0] ) for n in range(1,order)]

def ft(time):
    return sum( a*sin(x*n/2)*exp(-k*(n/2)^2*time) for n,a in alpha)
A = animate([p + plot(ft(0.003*j), 0, 2*pi, color='red', thickness=2) for j in range(50)], ymin=-.2)
A.show(delay=12)

x_coords = [cos(a)^3 for a in srange(0, 2*pi, 0.02)]
y_coords = [sin(a)^3 for a in srange(0, 2*pi, 0.02)]
list_plot(list(zip(x_coords, y_coords)))

colors = rainbow(30)
p = sum([plot(x^n, [x, -0.1,1.1], color=colors[n-1]) for n in range(1,31)])
p.show(ymax=1.2)

plot([x^n for n in range(1,31)], [x, -0.1,1.1], ymax=1.2)

plot(1.13*log(x), 1, 100, fill=lambda x: nth_prime(x)/floor(x), fillcolor='red')

polar_plot(2 + 3*cos(x), (x, 0, 2*pi), color=hue(0.3))

One of the powerful features of using CoCalc (or the Sage notebook) is through interacts, where one can dynamically update inputs in the same cell. These are not just limited to plots, but the examples we give below are all for plots.

For more information and examples, see:

# Example with torus knots
@interact
def _(m=slider(1,10,1, default=3),n=slider(1,10,1,default=5)):
    T = parametric_plot([(3 + cos(x))*cos(y), (3 + cos(x))*sin(y), sin(x)], [x, -pi, pi], [y,-pi,pi], opacity=0.7)
    K = parametric_plot([(3.05 + 1.1*cos(m*t))*cos(n*t), (3.05 + 1.1*cos(m*t))*sin(n*t), 1.1*sin(m*t)], [t, -pi, pi], plot_points=200, thickness=3, color='red')
    show(T + K)

# Heat equation (by Fourier series)
k   = 1
f   = (2*x^2-x)*exp(-x^2)
p   = plot(f,0,2*pi, thickness=1, linestyle='--')
c   = 1/pi
order = 20
alpha = [(n, c*numerical_integral(f*sin(x*n/2),0,2*pi)[0] ) for n in range(1,order)]

@interact
def _(time=(0.01*j for j in range(20)) ):
    ft = sum( a*sin(x*n/2)*exp(-k*(n/2)^2*time) for n,a in alpha)
    pt = plot(ft, 0, 2*pi, color='red', thickness=2)
    show( p + pt, ymin = -.2)

x0  = 0
f   = sin(x) * exp(-x)
p   = plot(f, -1,5, thickness=2)
dot = point((x0,f(x=x0)), pointsize=80, rgbcolor=(1,0,0))

@interact
def _(order=slider(1, 12, 1, default=2, label='order')):
  ft = f.taylor(x,x0,order)
  pt = plot(ft,-1, 5, color='green', thickness=2)
  html('$f(x)\;=\;%s$'%latex(f))
  html('$\hat{f}(x;%s)\;=\;%s+\mathcal{O}(x^{%s})$'%(x0,latex(ft),order+1))
  show(dot + p + pt, ymin = -.5, ymax = 1)

Linear algebra is one of the core components of any computer algebra system. Sage provides support for doing linear algebra over a variety of fields and some support over rings.

m = matrix(QQ, [[1,2],[-3,2]])
m

m.det()

~m

m^2

m.charpoly()

m.charpoly()(m)

m.eigenvectors_right()

view(m.exp())

M = matrix([[497,100,-60],[-2400,-483,291],[125,25,-13]])
M

view( M.jordan_form(transformation=True) )

view(M.LU())

Mp = (M - 2*matrix.identity(3))
Mp

Mp^2

(Mp^2).jordan_form()

Mp.right_kernel()

Mp.column_space()

m = random_matrix(RDF, 500)
e = m.eigenvalues()
w = [(i, abs(e[i])) for i in range(len(e))]
p = parametric_plot((x, sqrt(-1/3*(x - 500))), (x, 0, 500), color='red') + plot(points(w))
p.set_aspect_ratio('automatic')
p

m = matrix(GF(5), [[1,0,4],[3,4,1],[3,6,-2]])
m

m.right_kernel()

m.eigenvalues()

m.jordan_form()

m.rref()

v = vector(GF(5), [1,4,2])
v

m * v

4 * v

v * M

A = matrix(Zmod(8), [[1,2,0,4],[4,2,3,1],[0,3,4,1]])
A

B = matrix(Zmod(8), [[4,2,1],[3,8,2],[0,2,0],[2,2,4]])
B

A * B

B * A

((A*B).det(), (B*A).det())

Working with symbolic variables (i.e., by using var()) has a lot of flexibility, but because of that generality, they don't have a lot of structure and methods. The computations can also become very slow for more complicated expressions. Instead, if you want to work in a polynomial ring, then we should take advantage of that extra structure.

R.<x> = QQ[]
R
x

R == PolynomialRing(QQ, 'x')
R == QQ['x']

p = (x^3 + x^2) * (x^2 - x - 2)
p

p.factor()

p(2)

p(x + 2)

p.roots()

p / (x - 1)

Notice above that we also have rational functions and have created one just like we did with the rational numbers. We can also do quotients with remainders.

p / (x - 1) in R

p / (x - 2) in R

p / (x - 1) in R.fraction_field()

parent(p / (x - 2))

R(p / (x-2))

p // (x - 2)

parent(p // (x-2))

p // (x - 1)

p.quo_rem(x-1)

r = p / (x-1)
r.partial_fraction_decomposition()

r = p / ((x-1) * (x^2-3))
r.partial_fraction_decomposition()

To create a multivariate polynomial ring is no different than for univariate, but if we want an array of variables, here is an easier construction.

S = PolynomialRing(R, 'a', 10)
S

S.inject_variables()
(x * a0 + a1 - 1/2*a9) * (a2^3 * a6 - a7 + (x - 2*x^3) * a8)

Sage can also perform more advanced computations, such as Gröbner bases and quotient rings.

R.<x,y,z> = QQ[]
I = R.ideal([x^2 - y, y^2 - z, x*y^4*z + 2*x*y])

I.groebner_basis()

Q = R.quotient(I)
xb,yb,zb = Q.gens()
xb^2

zb^5

There are a number of other basic algebraic objects in Sage. We give a few of them here:

$\QQbar$ = QQbar
$p$ -adics Zp(p)
(universal) cyclotomic fields
power series (either with finite precision or lazily computed)
Laurent polynomials
finite fields

You can always start typing something you want and hit <tab> to see if it is in Sage. Note that capitalization matters!.

As previously mentioned, Sage is built upon Python, which means we can write our own functions and loops to test things and construct examples. Below, we compute the number of polynomials with increasing coefficients in degree 1 to 20 over $\GF{3}$ that are divisible by $x - 1$ .

We can then take this data and plug it into OEIS and find that it is http://oeis.org/A007997. If we do $\GF{5}$ , then we get http://oeis.org/A008646. For $\GF{7}$ , we obtain http://oeis.org/A032192. This then suggests we should try it for $\ZZ / 4\ZZ$ , and we get http://oeis.org/A008610.

R.<x> = Zmod(4)[]
data = []
import itertools
def test(p):
    return (x-1).divides(p)

for d in range(1,21):
    total = 0
    for t in itertools.combinations_with_replacement(R.base_ring(), d):
        if test(R(t)):
            total += 1
    data.append(total)
data

Sage also has a lot of builtin features for many different fields in mathematics, including:

(algebraic) combinatroics (including matroid theory, graph theory)
number theory
(algebraic) geometry
group theory

P = Permutations(3)
view(P.bruhat_poset())

G = groups.permutation.DiCyclic(3)
G.character_table()

Next, we turn to representation theory of the Lie algebra

\mathfrak{sl}_n

. When Sage contains a number of Lie algebras and related constructions by default. However, we implement our own

\mathfrak{sl}_2 = \{ M \in \mathrm{Mat}_{2 \times 2} \mid \operatorname{tr}(M) = 0 \}

. This has a natural basis of

E = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}, \qquad F = \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix}, \qquad H = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}.

We also implement its symmetric powers representations

V_m

, which are also the only finite-dimensional irreducible representations. Specifically, the irreducible representations

V_m

for

m \in \ZZ_{>0}

have a basis

\{v_{-m}, v_{-m+2}, \dotsc, v_{m-2}, v_{m}\}

with actions

\begin{align*} E v_k & = \left(\frac{k+m}{2} + 1\right) v_{k+2}, \\ F v_k & = \left(\frac{m-k}{2} + 1\right) v_{k-2}, \\ H v_k & = k v_k. \end{align*}

from sage.structure.element_wrapper import ElementWrapper
from sage.structure.parent import Parent
from sage.structure.unique_representation import UniqueRepresentation
from sage.categories.modules import Modules

class sl2(Parent, UniqueRepresentation):
    def __init__(self, R):
        Parent.__init__(self, R, category=Modules(R).WithBasis().FiniteDimensional())

    def _repr_(self):
        return "sl2 over {}".format(self.base_ring())

    @cached_method
    def basis(self):
        MS = MatrixSpace(self.base_ring(), 2)
        return Family({'E': self.element_class(self, MS([[0,1],[0,0]])),
                       'F': self.element_class(self, MS([[0,0],[1,0]])),
                       'H': self.element_class(self, MS([[1,0],[0,-1]]))})

    class Element(ElementWrapper):
        def _add_(self, rhs):
            return type(self)(self.parent(), self.value + rhs.value)

        def _sub_(self, rhs):
            return type(self)(self.parent(), self.value - rhs.value)

        def _neg_(self):
            return type(self)(self.parent(), -self.value)

        def _acted_upon_(self, scalar, self_on_left=True):
            if scalar in self.base_ring():
                return type(self)(self.parent(), scalar * self.value)
            return None

        def bracket(self, rhs):
            return type(self)(self.parent(), self.value * rhs.value - rhs.value * self.value)

# Note that to define the representation, we do not actually need to define the algebra.
# We only need to define its action!
class sl2IrreducibleRepresentation(CombinatorialFreeModule):
    def __init__(self, R, m):
        self._m = m
        CombinatorialFreeModule.__init__(self, R, range(-m, m+1, 2), prefix='v')

    def _repr_(self):
        return "Irreducible representation of highest weight {} of sl2 over {}".format(self._m, self.base_ring())

    def highest_weight_vector(self):
        return self.basis()[self._m]

    class Element(CombinatorialFreeModule.Element):
        def E(self):
            P = self.parent()
            d = self.monomial_coefficients(copy=False)
            return P._from_dict({k+2: c * (k+P._m+2) // 2 for k,c in d.items() if k < P._m})

        def F(self):
            P = self.parent()
            d = self.monomial_coefficients(copy=False)
            return P._from_dict({k-2: c * (P._m-k+2) // 2 for k,c in d.items() if k > -P._m})

        def H(self):
            d = self.monomial_coefficients(copy=False)
            return self.parent()._from_dict({k: c*k for k,c in d.items()})

        def _acted_upon_(self, scalar, self_on_left=True):
            if isinstance(scalar, sl2.Element):
                ec = scalar.value[0,1]
                fc = scalar.value[1,0]
                hc = scalar.value[0,0]
                return ec*self.E() + fc*self.F() + hc*self.H()
            return super(sl2IrreducibleRepresentation.Element, self)._acted_upon_(scalar, self_on_left=True)

L = sl2(QQ)
B = L.basis()
E,F,H = B['E'], B['F'], B['H']

E.bracket(F) == H

V = sl2IrreducibleRepresentation(QQ, 1)
v = V.highest_weight_vector()
v

vm = v.F()
vm

V = sl2IrreducibleRepresentation(QQ, 4)
v = V.highest_weight_vector()
v

v.F().F().F().F().E()

(E + 7*F - H) * v.F()

all(b.E().F() - b.F().E() == -b.H() for b in V.basis())

Next, if we want to generalize this to $\mathfrak{sl}_n$ (or any simple Lie algebra), we need to discuss root systems: a finite set of vectors closed under reflections in $\RR^n$ with no scalar multiples. The root systems for $\RR^2$ are $A_1 \times A_1$ , $A_2$ , $B_2$ , and $G_2$ .

You can find an extended tutorial for the visualization of root systems here: http://doc.sagemath.org/html/en/reference/combinat/sage/combinat/root_system/plot.html

RootSystem('A2').ambient_space().plot(roots='all')

RootSystem('B2').ambient_space().plot(roots='all')

 RootSystem('G2').ambient_space().plot(roots='all')

To every (crystallographic) root system, there is an associated Lie algebra. However, the $\mathfrak{sl}_2$ example we did before was a bit complicated. So we can simplify things by using crystals, which are certain edge-colored directed graphs. For $\mathfrak{sl}_n$ , they have an interpretation using partitions and semistandard tableaux with entries at most $n$ .

We first use Sage to list all such semistandard tableaux.

ascii_art(list(SemistandardTableaux([2,1], max_entry=3)))

C = crystals.Tableaux(['A',2], shape=[2,1])
latex.eval(latex(C))

This is built up by using the tensor product rule from the vector representation of $\mathfrak{sl}_n$ .

B = crystals.Letters(['A',2])
T = tensor([B,B,B])
latex.eval(latex(T))

T = crystalssor([C,C])
ascii_art(T.highest_weight_vectors())

In fact, as the partition gets larger (by containment), we have an embedding of crystals. The direct limit is called $B(\infty)$ and the crystal of every finite-dimensional representation is contained within $B(\infty)$ .

B = crystals.infinity.Tableaux(['A',2])
latex.eval(latex(B.subcrystal(max_depth=3)))

If we consider the (stable) characters of $\mathfrak{sl}_n$ , these are specific symmetric functions called Schur functions. More generally, Sage provides a very expansive implementation of symmetric functions, including all of the classical bases, the character bases of $S_n$ of recent work of Orellana-Zabrocki (i.e., structure coefficients are stable Kronecker coefficients), Hall-Littlewood polynomials, Jack polynomials, and Macdonald polynomials. It can also computer other products, such as the Kronecker product, plythesyms, Hopf algebra computations, and other related objects such as $k$ -Schur functions and Stanley symmetric functions.

For a more detailed tutorial, see this SageMath tutorial. We do some basic computations here to demonstrate how to work with symmetric functions.

R = FractionField(ZZ['t'])
t = R.gen()
Sym = SymmetricFunctions(R)
e = Sym.e()
e

s = Sym.s()
s

s[2,1] * s[1]

s[1]^3

s[2,1]^2

    e[2,1] + t^2*s[2]

e(s[3,2,1,1] + t*s[[]])   # The s[[]] needs the extra [] because Sage uses Python

s(t*e[4,1] + t^2*e[2,2,1] + (t - t^4)/(t-2*t^2) * e[5])

P = Sym.hall_littlewood(t=t).P()
P

s(P[2,1])

P(s[3,3,1])

Qp = Sym.hall_littlewood(t=t).Qp()
s(Qp[3,3,1])

P(Qp[2,1].omega())

s[3,1].plethysm(t*s[2] + t^2*s[1,1])

P[2,1].expand(4)

Sage also has generalizations of the ring of symmetric functions, including quasisymmetric functions (QSym), non-commutative symmetric functions (NCSF, which is Hopf dual to QSym), and symmetric functions in non-commutative variables (NCSym) (and its Hopf dual). As they have multiple bases, like symmetric functions, we need to specify the bases we want to work with.

For more detailed information, see

Product

Resources

Company