3619.sagews

n = 3
R = PolynomialRing(RR,[['r','t'][cmp(int(i/n),i%n)]+'_'+str(1+int(i/n))+str(1+i%n) for i in range(n^2)])
S = matrix(R,3,R.gens())
#print S*S
from operator import add
#print [S[int(k/n)][k%n] for k in range(n^2)]
#print flatten([[[(i,int(k/n),k%n,j,1*(int(k/n)!=k%n))for k in range(n^2)]for i in range(n)]for j in range(n)])
def Sh(i,j):
    return reduce(add,[S[i][k//n]*S[i][k%n]*S[j][k%n]*S[j][k//n]*((k//n)!=k%n) for k in range(n^2)])
Shot = matrix(R,n,n, Sh)
show(Shot)

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rrr} 2.00000000000000 r_{11}^{2} t_{12}^{2} + 2.00000000000000 r_{11}^{2} t_{13}^{2} + 2.00000000000000 t_{12}^{2} t_{13}^{2} & 2.00000000000000 r_{11} t_{12} t_{21} r_{22} + 2.00000000000000 r_{11} t_{13} t_{21} t_{23} + 2.00000000000000 t_{12} t_{13} r_{22} t_{23} & 2.00000000000000 r_{11} t_{12} t_{31} t_{32} + 2.00000000000000 r_{11} t_{13} t_{31} r_{33} + 2.00000000000000 t_{12} t_{13} t_{32} r_{33} \\ 2.00000000000000 r_{11} t_{12} t_{21} r_{22} + 2.00000000000000 r_{11} t_{13} t_{21} t_{23} + 2.00000000000000 t_{12} t_{13} r_{22} t_{23} & 2.00000000000000 t_{21}^{2} r_{22}^{2} + 2.00000000000000 t_{21}^{2} t_{23}^{2} + 2.00000000000000 r_{22}^{2} t_{23}^{2} & 2.00000000000000 t_{21} r_{22} t_{31} t_{32} + 2.00000000000000 t_{21} t_{23} t_{31} r_{33} + 2.00000000000000 r_{22} t_{23} t_{32} r_{33} \\ 2.00000000000000 r_{11} t_{12} t_{31} t_{32} + 2.00000000000000 r_{11} t_{13} t_{31} r_{33} + 2.00000000000000 t_{12} t_{13} t_{32} r_{33} & 2.00000000000000 t_{21} r_{22} t_{31} t_{32} + 2.00000000000000 t_{21} t_{23} t_{31} r_{33} + 2.00000000000000 r_{22} t_{23} t_{32} r_{33} & 2.00000000000000 t_{31}^{2} t_{32}^{2} + 2.00000000000000 t_{31}^{2} r_{33}^{2} + 2.00000000000000 t_{32}^{2} r_{33}^{2} \end{array}$$\right)

def kron(i,j):
    if i==j: return 1
    else: return 0
aa = [var('alpha_%d' % (i+1)) for i in range(n)]
show(aa)
#F=matrix(R,n,n,)

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\alpha_{1}, \alpha_{2}, \alpha_{3}\right]

help(sage.matrix.matrix_generic_dense.Matrix_generic_dense)

Click to open help window

xx = [var('alpha_%d'% (i+1)) for i in range(n)]
n = 3
R = PolynomialRing(RR,[['r','t'][cmp(int(i/n),i%n)]+'_'+str(1+int(i/n))+str(1+i%n) for i in range(n^2)])
S = matrix(R,3,R.gens())
def kron(i,j):
    if i==j: return 1
    else: return 0
def func(S,xx):
    return [[(-1/2)*xx[i]*S[i][j]*kron(i,j) for i in range(n)]for j in range(n)]
f = matrix(func(S,xx))
conjugate(f)*f.simplify()

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rrr} \frac{1}{4} \, \alpha_{1} r_{11} \overline{\alpha_{1}} \overline{r_{11}} & 0 & 0 \\ 0 & \frac{1}{4} \, \alpha_{2} r_{22} \overline{\alpha_{2}} \overline{r_{22}} & 0 \\ 0 & 0 & \frac{1}{4} \, \alpha_{3} r_{33} \overline{\alpha_{3}} \overline{r_{33}} \end{array}$$\right)

mat = matrix([[0,1],[1,0]])
u = exp(i*mat)
show(((u.conjugate_transpose())*u).n())

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rr} 1.00000000000000 & 0 \\ 0 & 1.00000000000000 \end{array}$$\right)

#fixed points
n=3
#print latex(matrix([[1*(i==j)for i in range(n)]for j in range(n)]))
def f(i,j,k):
    if (i!=k)and(i!=j):
        return 1*(i!=j)
    elif (k==i)and (i==j):
        return 1
    else :
        return 0
#print latex(matrix([[f(i,j,2) for i in range(n)] for j in range(n)]))
matrix([[1*(i==((j+1)%n))for i in range(n)]for j in range(n)])

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rrr} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{array}$$\right)

## thus it begins!
# Define S, F and the Governing Equations
xx = [var('alpha_%d'% (i+1)) for i in range(n)]
n = 3
R = PolynomialRing(RR,[['r','t'][cmp(int(i/n),i%n)]+'_'+str(1+int(i/n))+str(1+i%n) for i in range(n^2)])
S = matrix(R,3,R.gens())
def kron(i,j):
    if i==j: return 1
    else: return 0
def func(S,xx):
    return [[(-1/2)*xx[i]*S[i][j]*kron(i,j) for i in range(n)]for j in range(n)]
f = matrix(func(S,xx))
dSdl = S*f.conjugate_transpose()*S - f

# use subs via a dict as below
##print diff_S({S[0][0]:0})

n = 3
vnames = [('r' if i==j else 't') +'_%d%d' % (i,j) for i in [1..n] for j in [1..n]]
R = PolynomialRing(RR, vnames)
S = matrix(R, 3, R.gens())
f = S

mat = [[1*(i==((j+1)%n)) for i in [0..n-1]]for j in [0..n-1]]
#parent(mat)
#create dict
SubDict = dict([(S[p//n][p%n],mat[p//n][p%n]) for p in [0..((n^2)-1)]])
print f.subs(SubDict) #yay
# generate list of fixed point FPList
#define a function dSdl

[ 0 1.00000000000000 0] [ 0 0 1.00000000000000] [1.00000000000000 0 0]

#type(r_11)
var('a_45')
[a_45==0]
var('r_11')
S[0][0]==r_11

\newcommand{\Bold}[1]{\mathbf{#1}}r_{11} = r_{11}

def StabDirections(S_fp,dSdl,delta,n):
    PertDirList = [matrix([[delta*((i*n+j)==p) for i in range(n)] for j in range(n)]) for p in range(n^2)]
    return [(dSdl(S_fp + Pert) - dSdl(S_fp))/delta for Pert in PertDirList]
def PrettyPrintThroughFPList():
    return none
n=3
FPList = [matrix([[1*(i==j)for i in range(n)]for j in range(n)])]
def f(i,j,k):
    if (i!=k)and(j!=k):
        return 1*(i!=j)
    elif (k==i)and (i==j):
        return 1
    else :
        return 0
[FPList.append(matrix([[f(i,j,k) for i in range(n)] for j in range(n)])) for k in range(n)]
[FPList.append(matrix([[1*(i==((j+r)%n))for i in range(n)]for j in range(n)])) for r in [1,2]]
FPList
xx = [var('alpha_%d'% (i+1)) for i in range(n)]
def kron(i,j):
    if i==j: return 1
    else: return 0
a = 1/sum([1/i for i in xx ])
def listfunc(xx):
    return [[ -a*kron(i,j)/xx[i] + sqrt((1-a/xx[i])*(1-a/xx[j]))*(1-kron(i,j)) for i in range(n)]for j in range(n)]
FPList.append(matrix(listfunc(xx)))
FPList

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left($$\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}$$\right), \left($$\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array}$$\right), \left($$\begin{array}{rrr} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{array}$$\right), \left($$\begin{array}{rrr} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}$$\right), \left($$\begin{array}{rrr} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{array}$$\right), \left($$\begin{array}{rrr} 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}$$\right), \left($$\begin{array}{rrr} -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} \\ \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)}} \\ \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)}} & -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} \end{array}$$\right)\right]

xx = [var('alpha_%d'% (i+1)) for i in range(n)]
def kron(i,j):
    if i==j: return 1
    else: return 0
a = 1/sum([1/i for i in xx ])
def func(xx):
    return [[ -a*kron(i,j)/xx[i] + sqrt((1-a/xx[i])*(1-a/xx[j]))*(1-kron(i,j)) for i in range(n)]for j in range(n)]
matrix(func(xx))

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rrr} -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} \\ \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)}} \\ \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{1}} - 1\right)}} & \sqrt{{\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} - 1\right)} {\left(\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{2}} - 1\right)}} & -\frac{1}{{\left(\frac{1}{\alpha_{1}} + \frac{1}{\alpha_{2}} + \frac{1}{\alpha_{3}}\right)} \alpha_{3}} \end{array}$$\right)

#analytical fixed point calculation
# loaded previous constructs
col = matrix([S[k//n][k%n] for k in range(n^2)])
stab = matrix([dSdl[k//n][k%n] for k in range(n^2)])

stabdiff = [[ k[0].derivative(S[j//n][j%n]) for j in range(n^2)] for k in stab]
#stabdiff[0][0][0]

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{2} \, r_{11}^{2} \overline{\alpha_{1}} D[0]\left({\rm conjugate}\right)\left(r_{11}\right) - r_{11} \overline{\alpha_{1}} \overline{r_{11}} + \frac{1}{2} \, \alpha_{1}

Click to open help window

[S[0][0]==0]

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_8.py", line 10, in <module> exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("W1NbMF1bMF09PTBd"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmp9t6DtK/___code___.py", line 3, in <module> exec compile(u'[S[_sage_const_0 ][_sage_const_0 ]==_sage_const_0 ]' + '\n', '', 'single') File "", line 1, in <module> NameError: name 'S' is not defined