Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

"""Solvers of systems of polynomial equations. """ 

 

from __future__ import print_function, division 

 

from sympy.core import S 

from sympy.polys import Poly, groebner, roots 

from sympy.polys.polytools import parallel_poly_from_expr 

from sympy.polys.polyerrors import (ComputationFailed, 

PolificationFailed, CoercionFailed, PolynomialError) 

from sympy.simplify import rcollect 

from sympy.utilities import default_sort_key, postfixes 

 

 

class SolveFailed(Exception): 

"""Raised when solver's conditions weren't met. """ 

 

 

def solve_poly_system(seq, *gens, **args): 

""" 

Solve a system of polynomial equations. 

 

Examples 

======== 

 

>>> from sympy import solve_poly_system 

>>> from sympy.abc import x, y 

 

>>> solve_poly_system([x*y - 2*y, 2*y**2 - x**2], x, y) 

[(0, 0), (2, -sqrt(2)), (2, sqrt(2))] 

 

""" 

try: 

polys, opt = parallel_poly_from_expr(seq, *gens, **args) 

except PolificationFailed as exc: 

raise ComputationFailed('solve_poly_system', len(seq), exc) 

 

if len(polys) == len(opt.gens) == 2: 

f, g = polys 

 

if all(i <= 2 for i in f.degree_list() + g.degree_list()): 

try: 

return solve_biquadratic(f, g, opt) 

except SolveFailed: 

pass 

 

return solve_generic(polys, opt) 

 

 

def solve_biquadratic(f, g, opt): 

"""Solve a system of two bivariate quadratic polynomial equations. 

 

Examples 

======== 

 

>>> from sympy.polys import Options, Poly 

>>> from sympy.abc import x, y 

>>> from sympy.solvers.polysys import solve_biquadratic 

>>> NewOption = Options((x, y), {'domain': 'ZZ'}) 

 

>>> a = Poly(y**2 - 4 + x, y, x, domain='ZZ') 

>>> b = Poly(y*2 + 3*x - 7, y, x, domain='ZZ') 

>>> solve_biquadratic(a, b, NewOption) 

[(1/3, 3), (41/27, 11/9)] 

 

>>> a = Poly(y + x**2 - 3, y, x, domain='ZZ') 

>>> b = Poly(-y + x - 4, y, x, domain='ZZ') 

>>> solve_biquadratic(a, b, NewOption) 

[(-sqrt(29)/2 + 7/2, -sqrt(29)/2 - 1/2), (sqrt(29)/2 + 7/2, -1/2 + \ 

sqrt(29)/2)] 

""" 

G = groebner([f, g]) 

 

if len(G) == 1 and G[0].is_ground: 

return None 

 

if len(G) != 2: 

raise SolveFailed 

 

x, y = opt.gens 

p, q = G 

if not p.gcd(q).is_ground: 

# not 0-dimensional 

raise SolveFailed 

 

p = Poly(p, x, expand=False) 

p_roots = [ rcollect(expr, y) for expr in roots(p).keys() ] 

 

q = q.ltrim(-1) 

q_roots = list(roots(q).keys()) 

 

solutions = [] 

 

for q_root in q_roots: 

for p_root in p_roots: 

solution = (p_root.subs(y, q_root), q_root) 

solutions.append(solution) 

 

return sorted(solutions, key=default_sort_key) 

 

 

def solve_generic(polys, opt): 

""" 

Solve a generic system of polynomial equations. 

 

Returns all possible solutions over C[x_1, x_2, ..., x_m] of a 

set F = { f_1, f_2, ..., f_n } of polynomial equations, using 

Groebner basis approach. For now only zero-dimensional systems 

are supported, which means F can have at most a finite number 

of solutions. 

 

The algorithm works by the fact that, supposing G is the basis 

of F with respect to an elimination order (here lexicographic 

order is used), G and F generate the same ideal, they have the 

same set of solutions. By the elimination property, if G is a 

reduced, zero-dimensional Groebner basis, then there exists an 

univariate polynomial in G (in its last variable). This can be 

solved by computing its roots. Substituting all computed roots 

for the last (eliminated) variable in other elements of G, new 

polynomial system is generated. Applying the above procedure 

recursively, a finite number of solutions can be found. 

 

The ability of finding all solutions by this procedure depends 

on the root finding algorithms. If no solutions were found, it 

means only that roots() failed, but the system is solvable. To 

overcome this difficulty use numerical algorithms instead. 

 

References 

========== 

 

.. [Buchberger01] B. Buchberger, Groebner Bases: A Short 

Introduction for Systems Theorists, In: R. Moreno-Diaz, 

B. Buchberger, J.L. Freire, Proceedings of EUROCAST'01, 

February, 2001 

 

.. [Cox97] D. Cox, J. Little, D. O'Shea, Ideals, Varieties 

and Algorithms, Springer, Second Edition, 1997, pp. 112 

 

Examples 

======== 

 

>>> from sympy.polys import Poly, Options 

>>> from sympy.solvers.polysys import solve_generic 

>>> from sympy.abc import x, y 

>>> NewOption = Options((x, y), {'domain': 'ZZ'}) 

 

>>> a = Poly(x - y + 5, x, y, domain='ZZ') 

>>> b = Poly(x + y - 3, x, y, domain='ZZ') 

>>> solve_generic([a, b], NewOption) 

[(-1, 4)] 

 

>>> a = Poly(x - 2*y + 5, x, y, domain='ZZ') 

>>> b = Poly(2*x - y - 3, x, y, domain='ZZ') 

>>> solve_generic([a, b], NewOption) 

[(11/3, 13/3)] 

 

>>> a = Poly(x**2 + y, x, y, domain='ZZ') 

>>> b = Poly(x + y*4, x, y, domain='ZZ') 

>>> solve_generic([a, b], NewOption) 

[(0, 0), (1/4, -1/16)] 

""" 

def _is_univariate(f): 

"""Returns True if 'f' is univariate in its last variable. """ 

for monom in f.monoms(): 

if any(m for m in monom[:-1]): 

return False 

 

return True 

 

def _subs_root(f, gen, zero): 

"""Replace generator with a root so that the result is nice. """ 

p = f.as_expr({gen: zero}) 

 

if f.degree(gen) >= 2: 

p = p.expand(deep=False) 

 

return p 

 

def _solve_reduced_system(system, gens, entry=False): 

"""Recursively solves reduced polynomial systems. """ 

if len(system) == len(gens) == 1: 

zeros = list(roots(system[0], gens[-1]).keys()) 

return [ (zero,) for zero in zeros ] 

 

basis = groebner(system, gens, polys=True) 

 

if len(basis) == 1 and basis[0].is_ground: 

if not entry: 

return [] 

else: 

return None 

 

univariate = list(filter(_is_univariate, basis)) 

 

if len(univariate) == 1: 

f = univariate.pop() 

else: 

raise NotImplementedError("only zero-dimensional systems supported (finite number of solutions)") 

 

gens = f.gens 

gen = gens[-1] 

 

zeros = list(roots(f.ltrim(gen)).keys()) 

 

if not zeros: 

return [] 

 

if len(basis) == 1: 

return [ (zero,) for zero in zeros ] 

 

solutions = [] 

 

for zero in zeros: 

new_system = [] 

new_gens = gens[:-1] 

 

for b in basis[:-1]: 

eq = _subs_root(b, gen, zero) 

 

if eq is not S.Zero: 

new_system.append(eq) 

 

for solution in _solve_reduced_system(new_system, new_gens): 

solutions.append(solution + (zero,)) 

 

return solutions 

 

try: 

result = _solve_reduced_system(polys, opt.gens, entry=True) 

except CoercionFailed: 

raise NotImplementedError 

 

if result is not None: 

return sorted(result, key=default_sort_key) 

else: 

return None 

 

 

def solve_triangulated(polys, *gens, **args): 

""" 

Solve a polynomial system using Gianni-Kalkbrenner algorithm. 

 

The algorithm proceeds by computing one Groebner basis in the ground 

domain and then by iteratively computing polynomial factorizations in 

appropriately constructed algebraic extensions of the ground domain. 

 

Examples 

======== 

 

>>> from sympy.solvers.polysys import solve_triangulated 

>>> from sympy.abc import x, y, z 

 

>>> F = [x**2 + y + z - 1, x + y**2 + z - 1, x + y + z**2 - 1] 

 

>>> solve_triangulated(F, x, y, z) 

[(0, 0, 1), (0, 1, 0), (1, 0, 0)] 

 

References 

========== 

 

1. Patrizia Gianni, Teo Mora, Algebraic Solution of System of 

Polynomial Equations using Groebner Bases, AAECC-5 on Applied Algebra, 

Algebraic Algorithms and Error-Correcting Codes, LNCS 356 247--257, 1989 

 

""" 

G = groebner(polys, gens, polys=True) 

G = list(reversed(G)) 

 

domain = args.get('domain') 

 

if domain is not None: 

for i, g in enumerate(G): 

G[i] = g.set_domain(domain) 

 

f, G = G[0].ltrim(-1), G[1:] 

dom = f.get_domain() 

 

zeros = f.ground_roots() 

solutions = set([]) 

 

for zero in zeros: 

solutions.add(((zero,), dom)) 

 

var_seq = reversed(gens[:-1]) 

vars_seq = postfixes(gens[1:]) 

 

for var, vars in zip(var_seq, vars_seq): 

_solutions = set([]) 

 

for values, dom in solutions: 

H, mapping = [], list(zip(vars, values)) 

 

for g in G: 

_vars = (var,) + vars 

 

if g.has_only_gens(*_vars) and g.degree(var) != 0: 

h = g.ltrim(var).eval(dict(mapping)) 

 

if g.degree(var) == h.degree(): 

H.append(h) 

 

p = min(H, key=lambda h: h.degree()) 

zeros = p.ground_roots() 

 

for zero in zeros: 

if not zero.is_Rational: 

dom_zero = dom.algebraic_field(zero) 

else: 

dom_zero = dom 

 

_solutions.add(((zero,) + values, dom_zero)) 

 

solutions = _solutions 

 

solutions = list(solutions) 

 

for i, (solution, _) in enumerate(solutions): 

solutions[i] = solution 

 

return sorted(solutions, key=default_sort_key)