Ask Your Question

Daniel L's profile - activity

2020-03-09 21:34:19 -0500 commented question Unimodular matrices with additional restrictions

Good point. This certainly works if the dimension is 8 e.g. the $E_8$ root lattice gram matrix: E8 lattice wiki. I edited the post. Thank you

2020-03-09 12:56:08 -0500 asked a question Unimodular matrices with additional restrictions

I'd like to generate some unimodular matrices over ZZ with some restrictions:

  1. the entries are between -2 and 2
  2. the matrix is symmetric: $A = A^t$

Item (1) above can be addressed with the option upper_bound. Simply searching for matrices satisfying (1) and (2),

for j in [1..1000]:
A = random_matrix(ZZ,3,3, algorithm = 'unimodular', upper_bound = 3, max_tries = 1000)
    if A == A.transpose():
        A

becomes extremely inefficient, especially as the dimension grows. Is there a way to enforce items (1),(2) within the unimodular algorithm? Or another workaround?

Edit: I originally also wanted the diagonal entries fixed at 2, but this is a bit too restrictive.

2020-03-09 12:34:07 -0500 received badge  Notable Question (source)
2020-02-26 11:30:01 -0500 commented answer continuous Fourier transform

Thank you!

2020-02-26 11:29:50 -0500 commented answer continuous Fourier transform

Thank you!

2020-02-21 10:07:18 -0500 received badge  Nice Question (source)
2020-02-19 13:38:32 -0500 asked a question continuous Fourier transform

I have some basic code to plot a "nice" real-valued function and its continuous Fourier transform. For example,

F(x) = e^(-x^2)*(1-x^2)
plot2 = F.plot(xmin = -5, xmax = 5, color = 'red')
### continuous Fourier transform on RR
var('xi')
fhat(xi) = integral((F(x)*exp(-2*pi*I*x*xi)),(x,-oo,oo))
plot1 = fhat(xi).real().plot(xmin = -5, xmax = 5)   #put fhat and f in the same plot
show(plot1 + plot2, xmin = -5, xmax = 5, ymin = -1, ymax = 2)

f and fourier transform

Like many things in Sage, I appreciate it when the code looks like the formal math. However, this code is not robust. For instance, if we change the function to

F(x) = e^(-x^2)*sin(x)/x

Then I get an error message "Variable 'xi' not found". I'd like to improve this code to work with a wider class of functions.

Thanks,

2019-09-12 07:01:42 -0500 received badge  Popular Question (source)
2019-06-06 14:29:14 -0500 received badge  Popular Question (source)
2019-04-29 20:42:53 -0500 received badge  Popular Question (source)
2019-03-01 05:40:51 -0500 received badge  Nice Question (source)
2018-07-01 16:32:37 -0500 received badge  Enthusiast
2018-06-30 17:13:48 -0500 commented answer Solving a power series equation

The workaround on the above thread works well. Thanks for the help!

2018-06-30 15:58:02 -0500 commented answer Solving a power series equation

When I run this code I get the following traceback error: ValueError: rtol too small (4.5e-16 < 8.88178e-16) I'm not sure what this means or why it's happening.

2018-06-29 12:52:53 -0500 asked a question Solving a power series equation

I'd like to solve the equation f(t) == 1 where f is a power series:

var('n,t')
a1 = float(4 + sqrt(8)) 
b1 = float(3 + sqrt(8)) 
c1 = b1
an = 4*n*(n-1)+a1*n
an1 = 4*n*(n+1) + a1*(n+1)
anr = b1 + 2*n*(4*n+a1)
bn = 1+2*n*(n-1)+n*(b1-1)
bn1 = 1+2*n*(n+1)+(n+1)*(b1-1)
bnr = c1 + 2*n*(2*n+b1-1)

def f(t):
return (1+b1)^-t + 3/2*a1^-t + 1/2*(2*b1)^-t + b1^-t  + 2*sum(an^-t +an1^-t + anr^-t + 2*bn^-t +2*bn1^-t + bnr^-t ,n,1,1000).n()

I've tried simply plugging in values for t and I estimate that if f(t) = 1 then t is approximately 1.51. I'd like to possibly do something like

solve(f(t) == 1,t)

but this appears to lock sage up and I have to interrupt the process. The true f function is actually an infinite series, but I am truncating it to the first 1000 terms. I know there are methods of solving series equations such as the method of regula falsi. Does Sage have anything like this? Thanks!

2018-04-22 13:45:13 -0500 asked a question problem loading large file

I have a file which was generated and saved from a Sage notebook worksheet:

save(G,'G4D200')

where G is a list of vectors of length roughly 20 million, and the file 'G4D200.sobj` is about 370Mb.

When I try to load this file later to do some analysis

new_G = load('/path_of_file.../G4D200.sobj')

my computer goes deep into swap memory (with python eating it all up) and basically freezes. I have done the exact same thing for smaller file sizes (<10Mb) with no issues. I've tried waiting over 30 minutes and nothing seems to happen. Is this possibly a memory issue? Something with Sage? Possible workaround? I'm running the Sage notebook in fedora linux with an older laptop with 8Gb of RAM.

Thanks!

2018-04-01 22:13:06 -0500 commented answer 3d plot, slow

Great, thanks for the help!

2018-03-22 15:52:29 -0500 commented question 3d plot, slow

Also, is there a way to crank up the resolution so the spheres don't look pixelated?

2018-03-22 13:47:01 -0500 asked a question 3d plot, slow

I have some graphical objects (spheres and planes) and I'm trying to create a nice picture to export as a .png (other options?). The problem I'm having is that the interactive mode (right clicking on the graphic) is so slow that the menu options freeze up. I'd like to preview the image (rotate it around, zoom, etc.) before saving so that I can make sure I've got the right picture. What is a more efficient/ faster way to create and export a 3d picture? Thanks!

var('x,y,z')

E1 = implicit_plot3d(x^2+y^2+z^2==1/4, (x,-1/2,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'gray', opacity = 0.6)

E2 = implicit_plot3d((x-1)^2+y^2+z^2==1/4, (x,-1/2,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'gray')

E3 = implicit_plot3d((x-1/2)^2+(y-sqrt(3)/2)^2+z^2==1/4, (x,-1/2,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'gray')

E4 = implicit_plot3d(z == -1/2, (x,-1,4), (y,-1,4), (z,-0.6,0.5), color = 'gray', opacity = 0.2)

E6 = implicit_plot3d((x+1/2)^2+(y-sqrt(3)/2)^2+z^2==1/4, (x,-1,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'gray')

S1 = implicit_plot3d((x-1)^2+y^2+(z-1/2)^2==1, (x,-1,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'blue', opacity = 0.2)

S2 = implicit_plot3d((x+1/2)^2+(y-sqrt(3)/2)^2+(z+1/2)^2==1, (x,-1,3/2), (y,-1/2,sqrt(3)/2+1/2), (z,-0.5,0.5), color = 'blue', opacity = 0.2)

show(E1+E2+E3+E6+M1+S1+S2, aspect_ratio = 1, frame = false)

Also, if I try

t = E1+E2+E3+E6+M1+S1+S2

t.save('spheres.png', aspect_ratio = 1, frame = false) then the picture is quite distorted. Is aspect_ratio not what I want for 3d plots?

2018-02-01 09:25:51 -0500 commented answer random matrix with determinant +- 1

This works. Nice simple solution. Thanks Dan.

2018-01-31 23:27:00 -0500 commented answer random matrix with determinant +- 1

Right, that should work!...I'll let you know...

2018-01-31 22:15:26 -0500 asked a question random matrix with determinant +- 1

I want to generate a random 4x4 matrix with integer entries and determinant either 1 or -1. I know that you can use

random_matrix(ZZ,4,4, algorithm = 'unimodular')

to generate matrices with determinant 1 (so in the special linear group). However, I'm actually more interested in the matrices with determinant -1.

Is there a 'Sage' way to do this? Or are there other functions/routines out there I should look at?

Thanks!

2018-01-24 09:01:36 -0500 received badge  Popular Question (source)
2017-11-09 14:31:28 -0500 commented answer tree of vectors

On the last line tikz.pdf(), I don't get any output. I have dot2tex, graphviz, and slabbe installed. Any idea why this may be so?

2017-11-08 02:14:23 -0500 received badge  Nice Question (source)
2017-11-07 16:18:12 -0500 received badge  Commentator
2017-11-07 16:18:12 -0500 commented answer tree of vectors

Ok got it. In some cases, my thin_group may have loops and in other cases not. The visualization is very nice!

2017-11-07 15:49:01 -0500 commented answer tree of vectors

Very nice code, thank you. Can I ask, in line 4, why did you use random matrices?

2017-11-07 15:47:21 -0500 commented answer tree of vectors

Dan: very nice code. I'm still working on checking it with what I have here, but it appears to agree with my hand calculation for vectors of ht < 20. Thank you

2017-11-07 12:01:32 -0500 asked a question tree of vectors

This is a more complicated version of apply functions iteratively

I have a starting "seed" vector, say e1 = matrix(4,1,(1,0,0,0)) Next, I have a collection of 4x4 matrices

thin_group = [T1,T2,T3,T4,T5]

where

$T1 = \left(\begin{array}{rrrr} -1 & 0 & 0 & 0 \\ 2 & 1 & 0 & 0 \\ 4 & 0 & 1 & 0 \\ 4 & 0 & 0 & 1 \end{array}\right), T2 = \left(\begin{array}{rrrr} 1 & 2 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 4 & 1 & 0 \\ 0 & 4 & 0 & 1 \end{array}\right), T3 = \left(\begin{array}{rrrr} 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \end{array}\right), T4 = \left(\begin{array}{rrrr} 1 & 6 & 0 & 6 \\ 0 & 3 & 0 & 2 \\ 0 & 4 & 1 & 4 \\ 0 & -4 & 0 & -3 \end{array}\right), T5 = \left(\begin{array}{rrrr} 3 & 0 & 0 & 2 \\ 6 & 1 & 0 & 6 \\ 4 & 0 & 1 & 4 \\ -4 & 0 & 0 & -3 \end{array}\right)$

These matrices either fix e1 or send to a different vector. I have a "height" function, which just measures how large the vector entries are getting.

def h(X): return max(abs(X[0,0],abs(X[1,0]), abs(X[2,0]),abs(X[3,0]))

I want to develop the following procedure:

hmax = 50 Begin applying T1,...,T5 to e1. As soon as I get a different vector (T1*e1 = (-1,2,4,4)), measure the height of that vector. If the height is below hmax (h(-1,2,4,4) = 4 < hmax = 50), store that vector somewhere. Now apply another matrix from thin_group to (-1,2,4,4), (but not T1 again). If the new vector has a height below hmax, store it, and apply another matrix from thin_group, just not the previous one. Continue this process until we reach a vector with a height larger than hmax. Once this happens, go back to the previous vector and apply a different matrix from thin_group i.e. go up a different branch of the tree. Go up this branch until hmax is exceeded, then go up a different branch, etc.. Eventually, this process will generate all thin_group images of e1 with hmax < 50. (this is probably about 50,000 vectors)

This procedure should "remember" the vectors that were counted/stored, so I'm thinking that I may be able to index each vector in the tree with a sequence (a1,a2,a3,...,an) where aj=1,2,3,4,5 to denote which matrix T1,...,T5 we applied previously.

Any/all ideas or partial solutions greatly appreciated. Thanks!

2017-06-13 19:45:46 -0500 commented answer Help with unstable code

This is a nice solution to the problem. Thanks.

2017-06-13 19:45:04 -0500 commented answer Help with unstable code

Very nice code. I like how you made use of projective space. In doing so, you remind the user/reader that a FLT is just matrix multiplication. Thanks

2017-06-10 14:16:34 -0500 asked a question Help with unstable code

This basically continues from a previous post: generation of certain matrices

I'm trying to draw circles in the plane by using fractional linear transformations.

Here is my code:

K = NumberField(x^2 + 2, 's')
OK = K.ring_of_integers()
def FLT(M,z):
"""takes the fractional linear transformation/Mobius transformation of a z in CC"""
    if z == 'infinity' and M[1,0] != 0:
        return M[0,0]/M[1,0]
    elif z == 'infinity' and M[1,0] == 0:
        return 'infinity'
    elif M[1,0] != 0 and z == -M[1,1]/M[1,0]:
        return 'infinity'
    else:
        return (M[0,0]*z+M[0,1])/(M[1,0]*z + M[1,1])

Now generate circles based on 3 image points of a matrix M with entries in the ring OK via FLT:

var('a','b','r','x','y')
func_cir=(x-a)**2+(y-b)**2==r**2
circles = set()
j = 0
while j < 10:
    M = random_matrix(OK,2,algorithm = 'unimodular')
    if FLT(M,-1) != 'infinity' and FLT(M,0) != 'infinity' and FLT(M,1) != 'infinity':
    pta=[FLT(M,-1)[0],FLT(M,-1)[1]]
    ptb=[FLT(M,0)[0],FLT(M,0)[1]]
    ptc=[FLT(M,1)[0],FLT(M,1)[1]]    
    eq1 = func_cir.subs(x==pta[0]).subs(y==pta[1])
    eq2 = func_cir.subs(x==ptb[0]).subs(y==ptb[1])
    eq3 = func_cir.subs(x==ptc[0]).subs(y==ptc[1])
    sol = solve([eq1,eq2,eq3],a,b,r)
    C = circle((sol[1][0].rhs(),sol[1][1].rhs()),sol[1][2].rhs())    
    circles.add(C)
    j += 1

This second chunk of code throws the error "list index out of range" about 50% of the time. I thought that it may have been due to division by zero, but now I'm not so sure. I've also tried using a for loop with the exact same result. Thank you very much for the help!