Ask Your Question
0

Help with unstable code

asked 2017-06-10 21:16:34 +0100

Daniel L gravatar image

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!

edit retag flag offensive close merge delete

2 Answers

Sort by » oldest newest most voted
1

answered 2017-06-13 20:36:38 +0100

dan_fulea gravatar image

updated 2017-06-13 20:52:33 +0100

The following alternative code is running without errors. The Möbius transform is slightly changed, and the $\sqrt 2=|s|$ was added. (I think it is also wanted.)

K.<s> = QuadraticField( -2 )
OK    = K.ring_of_integers()
IP    = ProjectiveSpace( 1, K )

points = ( -1, 0, 1 )
def lift_point( x ):
    return vector( [ x, 1 ] )

def FLT( M, z ):
    """takes the fractional linear transformation/Mobius transformation of a z in CC"""
    Mz = IP( ( M * lift_point(z) ).list() )
    if Mz == IP( (1,0) ):
        return 'infinity'
    return Mz[0]

def components( z ):
    if z == 'infinity':
        return None
    return ( z[0], z[1]*sqrt(2) )

var('a','b','r','x','y');

func_cir = (x-a)**2 + (y-b)**2 == r**2

def get_circles( number_of_circles=10, show=False ):
    j = 0
    circles  = []
    while j < number_of_circles:

        M = random_matrix( OK, 2, algorithm = 'unimodular' )
        A, B, C = [ components( FLT( M, z ) ) for z in points ]

        print ( "\n%s\nM is %s\nA,B,C are %s, %s, %s ."
                % ( 30*'=', M.list(), A, B, C ) )
        noneOccured = bool( None in ( A, B, C ) ) 
        colinearityMatrix = matrix( 3,3, [ A[0], A[1], 1 ,
                                           B[0], B[1], 1 ,
                                           C[0], C[1], 1 ] ) if not noneOccured else None

        if noneOccured:
            print "!!! INFINITY... IGNORED !!!\n"

        elif bool( colinearityMatrix.det() == 0 ):
            # print ( "Are A,B,C colinear? %s" % colinearity )
            print "!!! COLINEAR... IGNORED !!!\n"

        else:
            j += 1

            eq1 = func_cir.subs( x==A[0] ).subs( y==A[1] )
            eq2 = func_cir.subs( x==B[0] ).subs( y==B[1] )
            eq3 = func_cir.subs( x==C[0] ).subs( y==C[1] )

            sol = solve( [ eq1, eq2, eq3 ], a,b,r )

            print "%s\n[%s]\nSolutions:\n%s" % ( 30*'-', j, sol )
            c = circle( ( sol[0][0].rhs(),
                          sol[0][1].rhs() ),
                          abs( sol[0][2].rhs() ) )
            circles.append( c )

    if show:
        plot_object = plot( [] )
        for c in circles:
            plot_object += plot( c )
        plot_object.show()
    return circles

Now the following calls did the job for me without errors:

get_circles( 1000 );

There is a verbose, aggressive print show in between, but the message is, there are no errors.

Let us also generate three circles, and show them:

get_circles( 3, show=True  );
edit flag offensive delete link more

Comments

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

Daniel L gravatar imageDaniel L ( 2017-06-14 02:45:04 +0100 )edit
1

answered 2017-06-12 11:38:27 +0100

ndomes gravatar image

The error "list index out of range" occurs if solve returns an empty list. You can capture this exception with a try except statement:

    try:
        C = circle((sol[1][0].rhs(),sol[1][1].rhs()),sol[1][2].rhs())    
    except:
        pass
edit flag offensive delete link more

Comments

This is a nice solution to the problem. Thanks.

Daniel L gravatar imageDaniel L ( 2017-06-14 02:45:46 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2017-06-10 21:16:34 +0100

Seen: 737 times

Last updated: Jun 13 '17