Ask Your Question
1

Solving a system of 5 polynomial equations

asked 2017-12-29 03:25:53 +0100

david_c gravatar image

updated 2018-01-01 16:21:33 +0100

I have a system of 5 polynomial equations. I solved the equations in Maple but I'm trying to solve them also in sagemath since it's open source. I've been trying for a long time (weeks) and still couldn't solve the equations in sagemath (Maple finds all the solutions in 10 seconds). I don't need the complex solutions (if exist). only the real solutions.

I have posted before a similar question but with 18 equations and it took Maple 40 minutes to solve them. But since then I was able to reduce the number of equations from 18 to 5. This time Maple solved the equations much faster (10 seconds instead of 40 minutes). But unfortunately, still can't solve the equations with sagemath.

The goal is to find all the possible (real) solutions of the equations. Maple did it and found 20 solutions (the code is attached bellow). But I wasn't able to solve the equations in sagemath at all (I attached the code for sagemath bellow as well).

I tried to use groebner basis in sagemath but there are so many different implementations with different order of the polynomials ('lex', 'degrevlex',etc). I tried many of them and none of them could solve the equations.

I will appreciate if anyone can suggest other ways to solve the equations. Maybe not groebner? any method that works and returns all the 20 solutions.

This is the code in Maple that defines the 5 equations and finds all the solutions:

restart;
eq1:=-18889.48706*qd+467.29186*qb+2982.844413*qd*qc*qb+30136.14351*qc-1115.07186*qc^2*qd-45629.75749*qa*qd*qb-23697.88597*qd^2*qc+4316.66628*qa*qc*qb-15135.14056*qb^2*qc-21851.63993*qa*qd*qc+7902.042467*qb^2*qd-5483.40637*qd^2*qb-11776.45729*qc^2*qb+40621.75090*qa+11348.58548*qd^3-18144.72833*qb^3-19544.23125*qc^3+2*b*qa+32521.65211*qa*qc^2+25623.87525*qa*qd^2+26225.65173*qa*qb^2;
eq2:=-29579.52853*qd+73144.02842*qb+8748.09605*qd*qc*qb+2240.035087*qc+26343.71460*qc^2*qd+2344.316752*qa*qd*qb-25331.91138*qd^2*qc+1444.808993*qa*qc*qb+6290.87549*qb^2*qc+2982.844337*qa*qd*qc+49189.98407*qb^2*qd-29432.25353*qd^2*qb-17559.89085*qc^2*qb-7338.23254*qa+2*b*qb+21944.51817*qd^3-11930.76271*qb^3+2394.88100*qc^3-3970.93290*qa*qc^2+2322.11803*qa*qd^2-23212.08740*qa*qb^2;
eq3:=-3375.795773*qd+2240.03510*qb+7057.67167*qd*qc*qb+56197.14569*qc+52577.03189*qc^2*qd+2982.84431*qa*qd*qb-17977.24342*qd^2*qc-136.341409*qa*qc*qb-23855.89112*qb^2*qc-15689.91200*qa*qd*qc+15299.86799*qb^2*qd-25331.91138*qd^2*qb+11501.30926*qc^2*qb+19564.44679*qa-963.74805*qd^3+658.069730*qb^3-24654.25363*qc^3-16345.90692*qa*qc^2-13126.18932*qa*qd^2-4563.44385*qa*qb^2+2*b*qc;
eq4:=37246.67872*qd-29579.52855*qb-46347.15644*qd*qc*qb-3375.795773*qc-11079.46666*qc^2*qd+12449.76046*qa*qd*qb-24742.88407*qd^2*qc+2982.844222*qa*qc*qb+15299.86800*qb^2*qc-15680.68187*qa*qd*qc-28830.47715*qb^2*qd+20203.79697*qd^2*qb+26343.71462*qc^2*qb-14402.89763*qa-14953.58829*qd^3+31606.58056*qb^3+24809.55728*qc^3-5601.66129*qa*qc^2+16099.39874*qa*qd^2+3415.45305*qa*qb^2+2*b*qd;
eq5:=qa^2+qb^2+qc^2+qd^2-1;
vars:=[b, qa, qb, qc, qd];
polysys:={eq1,eq2,eq3,eq4,eq5};
sols:=CodeTools:-Usage(RootFinding:-Isolate(polysys, vars, output = numeric, method = RS));

These are all the solutions that Maple found (20 solutions):

sols := 
[[b = -6104.435348, qa = -.2144080500, qb = -.2880353207, qc = -.1836576407, qd = -.9150599506], 
[b = -12366.14419, qa = -.5402828747, qb = -.3634451559, qc = .1841090309, qd = -.7362784111], 
[b = -24212.37926, qa = .2546299302, qb = .6836963191, qc = -.4499174642, qd = -.5150701091], 
[b = -23961.92853, qa = .2585879982, qb = .7021280140, qc = -.4278033463, qd = -.5070826324], 
[b = -6007.960508, qa = -.3010769760, qb = -.1196072448, qc = .8107192013, qd = -.4876280735], 
[b = -22102.25884, qa = .7018528580, qb = -.1525453855, qc = -.5362892523, qd = -.4433128793], 
[b = -26791.96626, qa = .3384082623, qb = -.1163579699, qc = -.8287561155, qd = -.4302371112], 
[b = -26003.78615, qa = 0.7746836766e-1, qb = -.1265934010, qc = -.9030243373, qd = -.4031374569], 
[b = -27290.86178, qa = .3563297191, qb = .8341751821, qc = .4195778650, qd = -0.3369439196e-1], 
[b = -39895.97736, qa = -.5758752754, qb = .7983140794, qc = -.1758240948, qd = -0.1217314478e-1], 
[b = -39895.97736, qa = .5758752754, qb = -.7983140794, qc = .1758240948, qd = 0.1217314478e-1], 
[b = -27290.86178, qa = -.3563297191, qb = -.8341751821, qc = -.4195778650, qd = 0.3369439196e-1], 
[b = -26003.78615, qa = -0.7746836766e-1, qb = .1265934010, qc = .9030243373, qd = .4031374569], 
[b = -26791.96626, qa = -.3384082623, qb = .1163579699, qc = .8287561155, qd = .4302371112], 
[b = -22102.25884, qa = -.7018528580, qb = .1525453855, qc = .5362892523, qd = .4433128793], 
[b = -6007.960508, qa = .3010769760, qb = .1196072448, qc = -.8107192013, qd = .4876280735], 
[b = -23961.92853, qa = -.2585879982, qb = -.7021280140, qc = .4278033463, qd = .5070826324], 
[b = -24212.37926, qa = -.2546299302, qb = -.6836963191, qc = .4499174642, qd = .5150701091], 
[b = -12366.14419, qa = .5402828747, qb = .3634451559, qc = -.1841090309, qd = .7362784111], 
[b = -6104.435348, qa = .2144080500, qb = .2880353207, qc = .1836576407, qd = .9150599506]]

And this is the code in sagemath that I couldn't solve. Same equations as before, only converted them to rational coefficients using eq_i=P(eq_i). It is stuck in the last line of code (I.variety(RR)).

P.<b, qa, qb, qc, qd>=PolynomialRing(QQ,order='degrevlex')
eq1=-18889.48706*qd+467.29186*qb+2982.844413*qd*qc*qb+30136.14351*qc-1115.07186*qc^2*qd-45629.75749*qa*qd*qb-23697.88597*qd^2*qc+4316.66628*qa*qc*qb-15135.14056*qb^2*qc-21851.63993*qa*qd*qc+7902.042467*qb^2*qd-5483.40637*qd^2*qb-11776.45729*qc^2*qb+40621.75090*qa+11348.58548*qd^3-18144.72833*qb^3-19544.23125*qc^3+2*b*qa+32521.65211*qa*qc^2+25623.87525*qa*qd^2+26225.65173*qa*qb^2
eq2=-29579.52853*qd+73144.02842*qb+8748.09605*qd*qc*qb+2240.035087*qc+26343.71460*qc^2*qd+2344.316752*qa*qd*qb-25331.91138*qd^2*qc+1444.808993*qa*qc*qb+6290.87549*qb^2*qc+2982.844337*qa*qd*qc+49189.98407*qb^2*qd-29432.25353*qd^2*qb-17559.89085*qc^2*qb-7338.23254*qa+2*b*qb+21944.51817*qd^3-11930.76271*qb^3+2394.88100*qc^3-3970.93290*qa*qc^2+2322.11803*qa*qd^2-23212.08740*qa*qb^2
eq3=-3375.795773*qd+2240.03510*qb+7057.67167*qd*qc*qb+56197.14569*qc+52577.03189*qc^2*qd+2982.84431*qa*qd*qb-17977.24342*qd^2*qc-136.341409*qa*qc*qb-23855.89112*qb^2*qc-15689.91200*qa*qd*qc+15299.86799*qb^2*qd-25331.91138*qd^2*qb+11501.30926*qc^2*qb+19564.44679*qa-963.74805*qd^3+658.069730*qb^3-24654.25363*qc^3-16345.90692*qa*qc^2-13126.18932*qa*qd^2-4563.44385*qa*qb^2+2*b*qc
eq4=37246.67872*qd-29579.52855*qb-46347.15644*qd*qc*qb-3375.795773*qc-11079.46666*qc^2*qd+12449.76046*qa*qd*qb-24742.88407*qd^2*qc+2982.844222*qa*qc*qb+15299.86800*qb^2*qc-15680.68187*qa*qd*qc-28830.47715*qb^2*qd+20203.79697*qd^2*qb+26343.71462*qc^2*qb-14402.89763*qa-14953.58829*qd^3+31606.58056*qb^3+24809.55728*qc^3-5601.66129*qa*qc^2+16099.39874*qa*qd^2+3415.45305*qa*qb^2+2*b*qd
eq5=qa^2+qb^2+qc^2+qd^2-1
eq1=P(eq1)
eq2=P(eq2)
eq3=P(eq3)
eq4=P(eq4)
eq5=P(eq5)
I=Ideal(eq1, eq2, eq3, eq4, eq5)
I.groebner_basis('libsingular:std')
I.variety(RR)

Anyone knows how to make it work?

Thanks

edit retag flag offensive close merge delete

Comments

For reference, the previous question by @david_c about a system of 18 polynomial equations is https://ask.sagemath.org/question/39998/solving-a-system-of-18-polynomial-equations-in-sagemath/ .

slelievre gravatar imageslelievre ( 2017-12-29 20:21:48 +0100 )edit

2 Answers

Sort by » oldest newest most voted
0

answered 2018-01-17 12:42:09 +0100

dan_fulea gravatar image

The following answer is rather showing why it is so hard to get the solutions by using algebra, then getting the numerical solutions. The questions insists to get an an algebraic solution, but the related algebra seems to become more and more complicated, so that it approaches the limits of the computable horizon. I will try to show why this is so. The code is maybe also useful for a numerical further refinement, although i had to truncate randomly to some point.

Algebraic approach

First of all, let us define some class, and methods that may be useful. My post uses the previous data, where the variables t11, t12, t13 were also present.

Some globals could also have been implemented as attributes of the class, but i tried to get a simple, quick way to proceed.

# ring with coefficients using a given approximation
FIELD = RealField( 200 )
R.<b,qa,qb,qc,qd,t11,t12,t13> = PolynomialRing( FIELD )

# ring with exact coefficients
S.<b,qa,qb,qc,qd,t11,t12,t13> = PolynomialRing( QQ )

# ==============================

f1 = ( -41.62655310*qb^2-17.84524397*qd*qc-54.90839708*qa*qd+80.88591073*qc*qb-36.57989776*qa*qc
       -38.37186364*qd*qb+103.9400553*qa*qb-43.34367897*qd^2+3.30995945*qc^2
       -3.030271028*t13-4.168000558*t12+13.08022556*t11-397.1837038 )

f2 = ( -2.8905538*qb^2-69.7800994*qd*qc+27.87345418*qa*qd+116.1536456*qc*qb-40.97241364*qa*qc
       -80.81288558*qd*qb+113.6742097*qa*qb+26.45817378*qd^2-55.32624148*qc^2
       -5.174596465*t13+8.630792519*t12-4.168000558*t11+126.0503892 )

f3 = ( -151.1500904*qb^2+13.88726142*qd*qc-62.96764170*qa*qd+23.12716969*qc*qb+68.70359049*qa*qc
       -73.73067590*qd*qb-70.78770770*qa*qb-21.05107218*qd^2-128.3070522*qc^2
       +10.28898193*t13-5.174596465*t12-3.030271028*t11+252.2534640 )

f4 = ( 48643.74394*qa*qd*qb+55369.81474*qa*qd*qc+83823.49379*qa*qb*qc+1068.121378*qc*qb*qd
       +2*b*qa+57179.40910*qa*qd^2-6744.998133*qc*qb^2-32360.07260*qc^2*qb-7128.296550*qb^2*qd
       +7808.793570*qd^2*qb+9053.831581*qd^2*qc+60527.22650*qa*qb^2+5696.900092*qd*qc^2
       +27.87345422*t12*qd+113.6742097*t12*qb-40.97241363*t12*qc+87381.51742*qa*qc^2
       -54.90839700*t11*qd+103.9400554*t11*qb-36.57989777*t11*qc-62.96764170*t13*qd-70.78770775*t13*qb
       +68.70359050*t13*qc-20550.50538*qb^3-22371.46454*qd^3+15873.58643*qc^3+16945.27914*qb
       +22351.84382*qd+15214.00155*qc )

f5 = ( 41911.74689*qa^2*qc+103.9400554*t11*qa+60527.22650*qa^2*qb+113.6742097*t12*qa
       -70.78770775*t13*qa+2*b*qb+24321.87195*qa^2*qd-14256.59311*qa*qd*qb+1068.121383*qa*qd*qc
       -13489.99629*qa*qb*qc+94398.45549*qc*qb*qd+7808.793561*qa*qd^2+18413.07448*qc*qb^2
       +190741.1431*qc^2*qb+27014.36446*qb^2*qd+181674.7133*qd^2*qb+16197.41252*qd^2*qc
       -61651.51619*qa*qb^2+23001.42975*qd*qc^2-80.81288568*t12*qd-5.781107439*t12*qb+116.1536457*t12*qc
       -32360.07262*qa*qc^2-38.37186355*t11*qd-83.25310616*t11*qb+80.88591081*t11*qc
       -73.73067575*t13*qd-302.3001806*t13*qb+23.12716970*t13*qc+196624.8969*qb^3+8254.040572*qd^3
       +11741.93330*qc^3-126998.3098*qb+16945.27914*qa-9234.551898*qd-19188.23950*qc )

f6 = ( 87381.51743*qa^2*qc-36.57989776*t11*qa+41911.74688*qa^2*qb-40.97241364*t12*qa
       +68.70359049*t13*qa+27684.90736*qa^2*qd+1068.121370*qa*qd*qb+11393.80020*qa*qd*qc
       -64720.14523*qa*qb*qc+46002.85945*qc*qb*qd+9053.831568*qa*qd^2+190741.1430*qc*qb^2
       +35225.79992*qc^2*qb+47199.22777*qb^2*qd+16197.41253*qd^2*qb+190164.0459*qd^2*qc
       -6744.998130*qa*qb^2+94531.22698*qd*qc^2-69.7800994*t12*qd+116.1536456*t12*qb
       -110.6524831*t12*qc+47620.75929*qa*qc^2-17.84524397*t11*qd+80.88591073*t11*qb
       +6.61991877*t11*qc+13.88726142*t13*qd+23.12716969*t13*qb-256.6141048*t13*qc
       +6137.691497*qb^3+43885.84331*qd^3+135209.5818*qc^3-19188.23950*qb+15214.00155*qa
       -46604.79008*qd-78396.52112*qc+2*b*qc )

f7 = ( 27684.90736*qa^2*qc-54.90839710*t11*qa+24321.87196*qa^2*qb+27.87345416*t12*qa
       -62.96764170*t13*qa+57179.40910*qa^2*qd+15617.58713*qa*qd*qb+18107.66316*qa*qd*qc
       +1068.121368*qa*qb*qc+32394.82511*qc*qb*qd-67114.39361*qa*qd^2+47199.22776*qc*qb^2
       +23001.42974*qc^2*qb+181674.7134*qb^2*qd+24762.12165*qd^2*qb+131657.5299*qd^2*qc
       -7128.296555*qa*qb^2+190164.0459*qd*qc^2+52.91634761*t12*qd-80.81288562*t12*qb
       -69.78009942*t12*qc+5696.900066*qa*qc^2-86.68735798*t11*qd-38.37186373*t11*qb
       -17.84524400*t11*qc-42.10214425*t13*qd-73.73067582*t13*qb+13.88726138*t13*qc+9004.788156*qb^3
       +133077.8789*qd^3+31510.40898*qc^3-9234.551898*qb+22351.84381*qa-92497.15079*qd
       -46604.79008*qc+2*b*qd )

f8 = ( qa^2+qb^2+qc^2+qd^2-1 )

# ==============================

class MyEquations(object):
    """simple rudimentary class to solve a system of equations
    """

    def __init__( self, exactRing, numericalRing, equations, usedEquations ):
        """Doc string by example:
        The exact ring is something like
        R.<x,y,z> = PolynomialRing( QQ )

        The *in*exact, "numerical" ring is something like
        S.<x,y,z> = PolynomialRing( RR )

        The caller of this constructor method has to take care that the 
        names of the variables concide, and have the same order.

        The "equations" are in fact polynomials, something like:
        ( x^2+y^2+z^2 - 1, x^3+y^3+z^3 - 1, x*y*z - 1 )
        and we consider the corresponding system, where each polynomial vanishes.

        The important method implemented in this class is a poor man elimination w.r.t. one variable,
        *without* building groebner bases. After we eliminate w.r.t. z by using x*y*z-1 == 0 for instance,
        we will remember this in usedEquations...

        """

        self.exactRing     = exactRing
        self.numericalRing = numericalRing
        self.equations     = equations
        self.usedEquations = usedEquations

        self.exactRingGens = self.exactRing.gens()

        return

    def eliminateVariable( self, variable, givenEquation, verbose=False ):
        """associate an instance of the same class,
        by using the poor man elimination (resultant) with respect to the given variable.
        """

        if givenEquation not in self.equations:
            print "Given equation is not accepted."
            return self

        S  = self.exactRing
        R  = self.numericalRing

        if S(variable) not in self.exactRingGens:
            print "Variable %s not found among the generators, returning self..." % variable
            return self

        placeOfVariable = S.gens().index( S(variable) )
        names_new = ','.join( [ str(gen) for gen in S.gens() if variable != gen ] )    # sorry for this quick+dirty solution
        if verbose:
            print ( "... eliminating variable %s in place %s from %s"
                    % ( variable, placeOfVariable, str(S.gens()) ) )
            print ( "... new names: %s" % names_new )

        Snew = PolynomialRing(    QQ, names=names_new )
        Rnew = PolynomialRing( FIELD, names=names_new )

        newEquations = [ eq.resultant( givenEquation, variable )
                         for eq in self.equations
                         if  eq != givenEquation ]

        newEquations = [ self.rewrite( eq, Snew, placeOfVariable )
                         for eq in newEquations ]    # this maps the eq as an element of Snew, polynomial ring with one less variable

        newEquations = [ Snew( Rnew( eq ) )
                         for eq in newEquations ]    # this gets a rough version of the coeffs

        newUsedEquations = {}
        for k, val in self.usedEquations.iteritems():
            newUsedEquations[k] = val
        newUsedEquations[ variable ] = S(givenEquation)

        return MyEquations( Snew, Rnew, newEquations, newUsedEquations )

    def rewrite( self, eq, Snew, placeOfVariable ):
        """eq should not contain the variable with position <placeOfVariable>.
        Then rewrite eq as an element of Snew...
        """
        S = self.exactRing
        Snew_gens = Snew.gens()

        eqnew = Snew(0)
        for powersTuple, coeff in eq.dict().items():
            powersTuple_new = [ powersTuple[k]
                                for k in range(len(powersTuple))
                                if  k != placeOfVariable ]
            if powersTuple[placeOfVariable]:
                print ( "*** the polynomial depends on the variable %s at place %s :: ERROR ***"
                        % ( S.gens()[placeOfVariable], placeOfVariable ) )
                return
            eqnew += coeff * prod( [ Snew_gens[k]^powersTuple_new[k]
                                     for k in range(len(powersTuple_new)) ] )
        return eqnew

    def adjustEquations( self ):
        """Sometimes, the coefficients of self.equations get so big,
        that exact computations take too long. Then we will renorm the coefficients,
        so that the minimal coefficient, taken in modulus, becomes one.
        This makes it simpler to proceed...
        """
        S = self.exactRing
        self.equations = [ S(self.renorm(eq)) for eq in self.equations ]

    def renorm( self, eq ):
        minimalAbsCoefficient = min( [ abs(c) for c in eq.coefficients() ] )
        S, R = self.exactRing, self.numericalRing
        if minimalAbsCoefficient:
            return S( R( eq/minimalAbsCoefficient ) )
        return eq

    def truncateEquationsByTotalDegree( self, totalDegree ):
        """replace the polynomials (i.e. equations) by truncated polynomials
        This may be ok for solutions in (-1,+1)
        """
        self.equations = [ self.truncate( eq, totalDegree )
                           for eq in self.equations ]

    def truncate( self, eq, totalDegree ):
        """Here there is some kind of a sage bug when the monomial is the "1",
        it delivers as coefficient the polynomial itself...
        In order to fix by patchwork the bug, 
        i'll use the method constant_coefficient for the mono = 1.

        There is a method eq.truncate, but it uses one single variable. 
        """
        S     = self.exactRing
        eqnew = S(0)
        gens  = S.gens()

        for powersTuple, coeff in eq.dict().items():
            if sum(powersTuple) > totalDegree:
                continue
            eqnew += coeff * prod( [ gens[k]^powersTuple[k] for k in range(len(powersTuple)) ] )

        return eqnew

    def roughInformationOnEquations( self, label='' ):
        count = 0
        print ( "%sInformation on %s\nGenerators: %s\nEquations degrees:\n"
                % ( label + ' :: ' if label else '', self, self.exactRing.gens() ) )
        for eq in self.equations:
            count += 1
            print "eq%s :: degrees = %s" % ( count, eq.degrees() )

rawEquations = [f1, f2, f3, f4, f5, f6, f7, f8]
equations    = [ S(eq) for eq in rawEquations ]

In this situation, we start with the instance:

X8 = MyEquations( S, R, equations, {} )
S8 = X8.exactRing

and perform explicit eliminations. The class was in fact designed to "avoid Gröbner" bases, since their computation has no in between information, and instead use the "poor man elimination" for such a case,which is building the resultant.

After the eliminations:

X7 = X8.eliminateVariable( S8(t11), X8.equations[0] )    # eliminate t11 using f1
S7 = X7.exactRing

X6 = X7.eliminateVariable( S7(t12), X7.equations[0] )    # eliminate t12 using the new f2
S6 = X6.exactRing

X5 = X6.eliminateVariable( S6(t13), X6.equations[0] )    # eliminate t13 using the new f3
S5 = X5.exactRing

X4 = X5.eliminateVariable( S5(b)  , X5.equations[0] )    # eliminate   b using the new f4
S4 = X4.exactRing

X3 = X4.eliminateVariable( S4(qa) , X4.equations[-1] )   # eliminate  qa using f8
S3 = X3.exactRing

X2 = X3.eliminateVariable( S3(qb) , X3.equations[0] )    # eliminate  qb using the new f5
S2 = X2.exactRing

We get the following

sage: X2.roughInformationOnEquations()
....: 
Information on <__main__.MyEquations object at 0x7f5abf7fc2d0>
Generators: (qc, qd)
Equations degrees:

eq1 :: degrees = (64, 64)
eq2 :: degrees = (64, 64)

A further elimination (of either variable) would produce a polynomial (with only even degree monomials) of degree $$ 64^2=4096\ ,$$ which is beyond my own time and computer resources. This is also the reason why building the Groebner basis did not work (in time). The variety could also not be computed for the same reason.

The following step distroys by design the hard earned equations, so we can store them for an other use...

X2equations = X2.equations

Since the solutions we are seeking for are in $[-1,1]$ (and maple does also show solutions near zero), i decided to implement a truncation method. This is no longer accurate, but i wanted only to proceed somehow. Please do something better instead. After:

X2 . truncateEquationsByTotalDegree( 23 )                # this is too rough, i know...
X2 . adjustEquations()

X2.roughInformationOnEquations()

We get:

Information on <__main__.MyEquations object at 0x7f5abf7fc2d0>
Generators: (qc, qd)
Equations degrees:

eq1 :: degrees = (22, 22)
eq2 :: degrees = (22, 22)

and there was a chance to eliminate qc. (Looong time...)

(One can also pass to a lower precision "numerical real field" at this point, to have a computation with smaller coefficients.)

X1 = X2.eliminateVariable( S2(qc) , X2.equations[0] )    # eliminate  qc using the new f6
S1 = X1.exactRing

Now we have only one equation, X1.equations[0], in only one remained variable, qd, and let us give the letter g for this one equation.

g = X1.equations[0]

It is now a simple matter to ask for the numerical real solutions of the one equation. After this, we have a cascade of equations, by substituting each qd-root in the polynomial joining qc and qd, then substituting each such obtained pair ( qc, qd ) to get all corresponding values for qb, than we use all triples ( qb, qc, qd ) ... and so on.

sage: g.roots( ring=AA )
[(-0.0270?, 1),
 (-0.03?, 1),
 (-0.0106105558612?, 1),
 (-0.0076740752917721?, 1),
 (0.0076740752917721?, 1),
 (0.0106105558612?, 1),
 (0.0252?, 1),
 (0.03?, 1)]

This is not satisfactory, but this is all i have "algebraically", and have to resign algebraically at this point.

Numerical approach

Let us start for simplicity from the last two equations obtained after repeated elimination of the variables. This is only to have a simpler situation to work with. We have saved the relatively accurate two equations in the global variable `X2equations. Let us use it now.

def X2eq( gc, gd ):
    maxCoeff = [ max( [ abs(c)
                        for c in eq.coefficients() ] )
                 for eq in X2equations ]
    g0, g1 = ( X2.equations[0] / maxCoeff[0] ,
               X2.equations[1] / maxCoeff[1] , )
    return [ g0(gc, gd) ,
             g1(gc, gd) , ]

Then we can use mpmath to get a common root for the equations g0 == 0 , g1 == 0, e.g. by starting from $(0,0)$. The following code finds such a tuple:

import mpmath
foundroot = mpmath.findroot( X2eq, (0,0) )
sol = ( RR( foundroot[0,0] ) ,
        RR( foundroot[1,0] ) , )
print sol

And the solution is:

(0.0114346341009851, -0.00153039811718345)

We can check it, by evaluation of X2eq in this tuple:

sage: X2eq( *sol )
[2.69913711193955e-610, -2.91585386517266e-610]

This may encourage us to search numerically also for the solution of the initial problem. Let us find a root for the given four equation, starting from the point ( 1000, 0.1, 0.1, 0.1, 0.1 ):

def f( b, qa, qb, qc, qd ):
    """define a function returning the values of the given equations.
    """
    eq1 = ( -18889.48706*qd + 467.29186*qb + 2982.844413*qd*qc*qb + 30136.14351*qc - 1115.07186*qc^2*qd
            - 45629.75749*qa*qd*qb - 23697.88597*qd^2*qc + 4316.66628*qa*qc*qb - 15135.14056*qb^2*qc
            - 21851.63993*qa*qd*qc + 7902.042467*qb^2*qd - 5483.40637*qd^2*qb - 11776.45729*qc^2*qb
            + 40621.75090*qa + 11348.58548*qd^3 - 18144.72833*qb^3 - 19544.23125*qc^3 + 2*b*qa
            + 32521.65211*qa*qc^2 + 25623.87525*qa*qd^2 + 26225.65173*qa*qb^2 )
    eq2 = ( - 29579.52853*qd + 73144.02842*qb + 8748.09605*qd*qc*qb + 2240.035087*qc + 26343.71460*qc^2*qd
            + 2344.316752*qa*qd*qb - 25331.91138*qd^2*qc + 1444.808993*qa*qc*qb + 6290.87549*qb^2*qc
            + 2982.844337*qa*qd*qc + 49189.98407*qb^2*qd - 29432.25353*qd^2*qb - 17559.89085*qc^2*qb
            - 7338.23254*qa + 2*b*qb + 21944.51817*qd^3 - 11930.76271*qb^3 + 2394.88100*qc^3
            - 3970.93290*qa*qc^2 + 2322.11803*qa*qd^2 - 23212.08740*qa*qb^2 )
    eq3 = ( - 3375.795773*qd + 2240.03510*qb + 7057.67167*qd*qc*qb + 56197.14569*qc + 52577.03189*qc^2*qd
            + 2982.84431*qa*qd*qb - 17977.24342*qd^2*qc - 136.341409*qa*qc*qb - 23855.89112*qb^2*qc
            - 15689.91200*qa*qd*qc + 15299.86799*qb^2*qd - 25331.91138*qd^2*qb + 11501.30926*qc^2*qb
            + 19564.44679*qa - 963.74805*qd^3 + 658.069730*qb^3 - 24654.25363*qc^3 - 16345.90692*qa*qc^2
            - 13126.18932*qa*qd^2 - 4563.44385*qa*qb^2 + 2*b*qc )
    eq4 = ( 37246.67872*qd - 29579.52855*qb - 46347.15644*qd*qc*qb - 3375.795773*qc - 11079.46666*qc^2*qd
            + 12449.76046*qa*qd*qb - 24742.88407*qd^2*qc + 2982.844222*qa*qc*qb + 15299.86800*qb^2*qc
            - 15680.68187*qa*qd*qc - 28830.47715*qb^2*qd + 20203.79697*qd^2*qb + 26343.71462*qc^2*qb
            - 14402.89763*qa - 14953.58829*qd^3 + 31606.58056*qb^3 + 24809.55728*qc^3 - 5601.66129*qa*qc^2
            + 16099.39874*qa*qd^2 + 3415.45305*qa*qb^2 + 2*b*qd )
    eq5 = ( qa^2 + qb^2 + qc^2 + qd^2 - 1.0 )

    M = 10.0 ** 4
    return( eq1/M, eq2/M, eq3/M, eq4/M, eq5/M )

import mpmath
foundroot    = mpmath.findroot( f, (1000,0.1,0.1,0.1,0.1) )
foundroot_RR = [ RR( foundroot[k,0] ) for k in range(5) ]

print "found root:\n  x  = %s\nf(x) = %s" % ( foundroot_RR, f(*foundroot_RR) )

We get the solution:

found root:
  x  = [-12366.1441874536, 0.540282874742273, 0.363445155880032, -0.184109030925670, 0.736278411104913]
f(x) = (-2.20552465179935e-15, -2.70574673777446e-15, 2.09183781407773e-15, -4.72937244921923e-15, 0.000000000000000)

Of course, we can do the same by picking random numbers from $[-10^5,10^5]\times [-1,1]^4$ and ask for the solution obtained by mpmath after starting from this point. I tried:

import random

def getRandomRoot():
    root = mpmath.findroot( f,
                            [ random.uniform( -10**5, 10**5 ) , ] +
                            [ random.uniform( -1    , 1     ) for k in range(4) ] )
    return [ RR(root[k,0]) for k in range(5) ]

randomRoots = [] 
for _ in range(100):
    root = None
    try:
        root = getRandomRoot()
    except Exception:
        print '*',
    if root and root not in randomRoots:
        randomRoots.append( root )
print
for root in randomRoots:
    print root

This gave me some roots, some of them repeat themselves up to some numerical tolerance, so they should be considered only once. Maybe not all roots are "attractors" for the iteration used by mpmath.

randomRoots.sort( lambda r, s: cmp( tuple( [ abs( r[k] ) for k in [1..4] ] ) ,
                                    tuple( [ abs( s[k] ) for k in [1..4] ] ) ) )

for root in randomRoots:
    print root

gives

[-26003.7861514443, 0.0774683676048174, -0.126593401011236, -0.903024337294491, -0.403137456845570]
[-26003.7861515540, -0.0774683676568282, 0.126593401010424, 0.903024337278419, 0.403137456852120]
[-26003.7861515662, -0.0774683676608112, 0.126593401010295, 0.903024337277287, 0.403137456852538]
[-26003.7861515662, -0.0774683676608113, 0.126593401010295, 0.903024337277287, 0.403137456852538]
[-26003.7861515662, 0.0774683676608113, -0.126593401010295, -0.903024337277287, -0.403137456852538]
[-26003.7861515662, 0.0774683676608117, -0.126593401010295, -0.903024337277288, -0.403137456852538]
[-6104.43534771621, -0.214408050049576, -0.288035320711193, -0.183657640670246, -0.915059950560715]
[-6104.43534771621, 0.214408050049576, 0.288035320711193, 0.183657640670246, 0.915059950560715]
[-24212.3792558492, -0.254629930206882, -0.683696319082174, 0.449917464205924, 0.515070109129656]
[-23961.9285295158, -0.258587998194318, -0.702128013983053, 0.427803346276447, 0.507082632402988]
[-23961.9285286599, -0.258587998207774, -0.702128014045988, 0.427803346200981, 0.507082632375747]
[-6007.96050751224, -0.301076976045761, -0.119607244780659, 0.810719201312935, -0.487628073549540]
[-6007.96050751224, 0.301076976045761, 0.119607244780659, -0.810719201312935, 0.487628073549540]
[-6007.96050751223, 0.301076976045761, 0.119607244780659, -0.810719201312936, 0.487628073549540]
[-26791.9662641179, 0.338408262304109, -0.116357969919450, -0.828756115483396, -0.430237111241477]
[-26791.9662641179, -0.338408262304109, 0.116357969919450, 0.828756115483396, 0.430237111241477]
[-27290.8617785685, 0.356329719114358, 0.834175182081041, 0.419577865033598, -0.0336943919635282]
[-27290.8617785685, -0.356329719114358, -0.834175182081041, -0.419577865033598, 0.0336943919635282]
[-27290.8617785685, -0.356329719114358, -0.834175182081042, -0.419577865033596, 0.0336943919635290]
[-27290.8617785685, 0.356329719114358, 0.834175182081042, 0.419577865033598, -0.0336943919635283]
[-12366.1441874536, 0.540282874742273, 0.363445155880032, -0.184109030925670, 0.736278411104913]
[-12366.1441874536, -0.540282874742273, -0.363445155880032, 0.184109030925670, -0.736278411104913]
[-12366.1441876125, -0.540282874753202, -0.363445155884730, 0.184109030921865, -0.736278411105592]
[-39895.9773624027, -0.575875275404261, 0.798314079436882, -0.175824094756561, -0.0121731447845349]
[-39895.9773624027, -0.575875275404262, 0.798314079436882, -0.175824094756561, -0.0121731447845347]
[-39895.9773624027, 0.575875275404262, -0.798314079436882, 0.175824094756561, 0.0121731447845347]
[-22102.2588367256, 0.701852857997968, -0.152545385483075, -0.536289252273274, -0.443312879334761]
[-22102.2588367256, -0.701852857997968, 0.152545385483075, 0.536289252273274, 0.443312879334761]

(Sorting was done w.r.t. the absolute values of the numbers in positions 1,2,3,4,)

I have to stop here with the conclusion that only a numerical insight (possibly combined with some algebraic initial work) can deliver all solutions. (Cheating by starting with the one or the other root delivered by maple should be the last thought. But the idea of starting with one maple root to see if the iteration stabilizes or not at the same point may give further interesting stuff to mention. Note that my definition of $f$ also involves some constant M, i needed it to some point to insure convergence while code bug-fixes were done, a posteriori this is possibly no longer needed, or the weighting may be adapted...)

edit flag offensive delete link more
0

answered 2018-01-22 21:36:20 +0100

tmonteil gravatar image

Looking at this example and your previous one, it seems to me that Groebner basis algorithms provided by Sage (singular, giac) are not powerful enough compared to the current state-of-art. It is possible that open-source software that are not (not yet, or not anymore) shipped with Sage can do the job. Could you please test the following two software and report if they succeed ?

If one of them work, i might be tempted to package them for Sage and add an interface.

edit flag offensive delete link more

Comments

Note that an interface to Macaulay2 exists in Sage. It will work if you have Macaulay2 installed and Sage can find it. See the wiki page on Sage interfaces.

slelievre gravatar imageslelievre ( 2018-01-29 18:32:21 +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-12-29 03:25:53 +0100

Seen: 1,341 times

Last updated: Jan 22 '18