Revision history [back]
 | 1 | initial version |
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...)
