Coverage for /usr/lib/python3/dist-packages/sympy/solvers/ode/lie_group.py: 5%

1r""" 

2This module contains the implementation of the internal helper functions for the lie_group hint for 

3dsolve. These helper functions apply different heuristics on the given equation 

4and return the solution. These functions are used by :py:meth:`sympy.solvers.ode.single.LieGroup` 

5 

6References 

7========= 

8 

9- `abaco1_simple`, `function_sum` and `chi` are referenced from E.S Cheb-Terrab, L.G.S Duarte 

10and L.A,C.P da Mota, Computer Algebra Solving of First Order ODEs Using 

11Symmetry Methods, pp. 7 - pp. 8 

12 

13- `abaco1_product`, `abaco2_similar`, `abaco2_unique_unknown`, `linear` and `abaco2_unique_general` 

14are referenced from E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order 

15ODE Patterns, pp. 7 - pp. 12 

16 

17- `bivariate` from Lie Groups and Differential Equations pp. 327 - pp. 329 

18 

19""" 

20from itertools import islice 

21 

22from sympy.core import Add, S, Mul, Pow 

23from sympy.core.exprtools import factor_terms 

24from sympy.core.function import Function, AppliedUndef, expand 

25from sympy.core.relational import Equality, Eq 

26from sympy.core.symbol import Symbol, Wild, Dummy, symbols 

27from sympy.functions import exp, log 

28from sympy.integrals.integrals import integrate 

29from sympy.polys import Poly 

30from sympy.polys.polytools import cancel, div 

31from sympy.simplify import (collect, powsimp, # type: ignore 

32 separatevars, simplify) 

33from sympy.solvers import solve 

34from sympy.solvers.pde import pdsolve 

35 

36from sympy.utilities import numbered_symbols 

37from sympy.solvers.deutils import _preprocess, ode_order 

38from .ode import checkinfsol 

39 

40 

41lie_heuristics = ( 

42 "abaco1_simple", 

43 "abaco1_product", 

44 "abaco2_similar", 

45 "abaco2_unique_unknown", 

46 "abaco2_unique_general", 

47 "linear", 

48 "function_sum", 

49 "bivariate", 

50 "chi" 

51 ) 

52 

53 

54def _ode_lie_group_try_heuristic(eq, heuristic, func, match, inf): 

55 

56 xi = Function("xi") 

57 eta = Function("eta") 

58 f = func.func 

59 x = func.args[0] 

60 y = match['y'] 

61 h = match['h'] 

62 tempsol = [] 

63 if not inf: 

64 try: 

65 inf = infinitesimals(eq, hint=heuristic, func=func, order=1, match=match) 

66 except ValueError: 

67 return None 

68 for infsim in inf: 

69 xiinf = (infsim[xi(x, func)]).subs(func, y) 

70 etainf = (infsim[eta(x, func)]).subs(func, y) 

71 # This condition creates recursion while using pdsolve. 

72 # Since the first step while solving a PDE of form 

73 # a*(f(x, y).diff(x)) + b*(f(x, y).diff(y)) + c = 0 

74 # is to solve the ODE dy/dx = b/a 

75 if simplify(etainf/xiinf) == h: 

76 continue 

77 rpde = f(x, y).diff(x)*xiinf + f(x, y).diff(y)*etainf 

78 r = pdsolve(rpde, func=f(x, y)).rhs 

79 s = pdsolve(rpde - 1, func=f(x, y)).rhs 

80 newcoord = [_lie_group_remove(coord) for coord in [r, s]] 

81 r = Dummy("r") 

82 s = Dummy("s") 

83 C1 = Symbol("C1") 

84 rcoord = newcoord[0] 

85 scoord = newcoord[-1] 

86 try: 

87 sol = solve([r - rcoord, s - scoord], x, y, dict=True) 

88 if sol == []: 

89 continue 

90 except NotImplementedError: 

91 continue 

92 else: 

93 sol = sol[0] 

94 xsub = sol[x] 

95 ysub = sol[y] 

96 num = simplify(scoord.diff(x) + scoord.diff(y)*h) 

97 denom = simplify(rcoord.diff(x) + rcoord.diff(y)*h) 

98 if num and denom: 

99 diffeq = simplify((num/denom).subs([(x, xsub), (y, ysub)])) 

100 sep = separatevars(diffeq, symbols=[r, s], dict=True) 

101 if sep: 

102 # Trying to separate, r and s coordinates 

103 deq = integrate((1/sep[s]), s) + C1 - integrate(sep['coeff']*sep[r], r) 

104 # Substituting and reverting back to original coordinates 

105 deq = deq.subs([(r, rcoord), (s, scoord)]) 

106 try: 

107 sdeq = solve(deq, y) 

108 except NotImplementedError: 

109 tempsol.append(deq) 

110 else: 

111 return [Eq(f(x), sol) for sol in sdeq] 

112 

113 

114 elif denom: # (ds/dr) is zero which means s is constant 

115 return [Eq(f(x), solve(scoord - C1, y)[0])] 

116 

117 elif num: # (dr/ds) is zero which means r is constant 

118 return [Eq(f(x), solve(rcoord - C1, y)[0])] 

119 

120 # If nothing works, return solution as it is, without solving for y 

121 if tempsol: 

122 return [Eq(sol.subs(y, f(x)), 0) for sol in tempsol] 

123 return None 

124 

125 

126def _ode_lie_group( s, func, order, match): 

127 

128 heuristics = lie_heuristics 

129 inf = {} 

130 f = func.func 

131 x = func.args[0] 

132 df = func.diff(x) 

133 xi = Function("xi") 

134 eta = Function("eta") 

135 xis = match['xi'] 

136 etas = match['eta'] 

137 y = match.pop('y', None) 

138 if y: 

139 h = -simplify(match[match['d']]/match[match['e']]) 

140 y = y 

141 else: 

142 y = Dummy("y") 

143 h = s.subs(func, y) 

144 

145 if xis is not None and etas is not None: 

146 inf = [{xi(x, f(x)): S(xis), eta(x, f(x)): S(etas)}] 

147 

148 if checkinfsol(Eq(df, s), inf, func=f(x), order=1)[0][0]: 

149 heuristics = ["user_defined"] + list(heuristics) 

150 

151 match = {'h': h, 'y': y} 

152 

153 # This is done so that if any heuristic raises a ValueError 

154 # another heuristic can be used. 

155 sol = None 

156 for heuristic in heuristics: 

157 sol = _ode_lie_group_try_heuristic(Eq(df, s), heuristic, func, match, inf) 

158 if sol: 

159 return sol 

160 return sol 

161 

162 

163def infinitesimals(eq, func=None, order=None, hint='default', match=None): 

164 r""" 

165 The infinitesimal functions of an ordinary differential equation, `\xi(x,y)` 

166 and `\eta(x,y)`, are the infinitesimals of the Lie group of point transformations 

167 for which the differential equation is invariant. So, the ODE `y'=f(x,y)` 

168 would admit a Lie group `x^*=X(x,y;\varepsilon)=x+\varepsilon\xi(x,y)`, 

169 `y^*=Y(x,y;\varepsilon)=y+\varepsilon\eta(x,y)` such that `(y^*)'=f(x^*, y^*)`. 

170 A change of coordinates, to `r(x,y)` and `s(x,y)`, can be performed so this Lie group 

171 becomes the translation group, `r^*=r` and `s^*=s+\varepsilon`. 

172 They are tangents to the coordinate curves of the new system. 

173 

174 Consider the transformation `(x, y) \to (X, Y)` such that the 

175 differential equation remains invariant. `\xi` and `\eta` are the tangents to 

176 the transformed coordinates `X` and `Y`, at `\varepsilon=0`. 

177 

178 .. math:: \left(\frac{\partial X(x,y;\varepsilon)}{\partial\varepsilon 

179 }\right)|_{\varepsilon=0} = \xi, 

180 \left(\frac{\partial Y(x,y;\varepsilon)}{\partial\varepsilon 

181 }\right)|_{\varepsilon=0} = \eta, 

182 

183 The infinitesimals can be found by solving the following PDE: 

184 

185 >>> from sympy import Function, Eq, pprint 

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

187 >>> xi, eta, h = map(Function, ['xi', 'eta', 'h']) 

188 >>> h = h(x, y) # dy/dx = h 

189 >>> eta = eta(x, y) 

190 >>> xi = xi(x, y) 

191 >>> genform = Eq(eta.diff(x) + (eta.diff(y) - xi.diff(x))*h 

192 ... - (xi.diff(y))*h**2 - xi*(h.diff(x)) - eta*(h.diff(y)), 0) 

193 >>> pprint(genform) 

194 /d d \ d 2 d 

195 |--(eta(x, y)) - --(xi(x, y))|*h(x, y) - eta(x, y)*--(h(x, y)) - h (x, y)*--(x 

196 \dy dx / dy dy 

197 <BLANKLINE> 

198 d d 

199 i(x, y)) - xi(x, y)*--(h(x, y)) + --(eta(x, y)) = 0 

200 dx dx 

201 

202 Solving the above mentioned PDE is not trivial, and can be solved only by 

203 making intelligent assumptions for `\xi` and `\eta` (heuristics). Once an 

204 infinitesimal is found, the attempt to find more heuristics stops. This is done to 

205 optimise the speed of solving the differential equation. If a list of all the 

206 infinitesimals is needed, ``hint`` should be flagged as ``all``, which gives 

207 the complete list of infinitesimals. If the infinitesimals for a particular 

208 heuristic needs to be found, it can be passed as a flag to ``hint``. 

209 

210 Examples 

211 ======== 

212 

213 >>> from sympy import Function 

214 >>> from sympy.solvers.ode.lie_group import infinitesimals 

215 >>> from sympy.abc import x 

216 >>> f = Function('f') 

217 >>> eq = f(x).diff(x) - x**2*f(x) 

218 >>> infinitesimals(eq) 

219 [{eta(x, f(x)): exp(x**3/3), xi(x, f(x)): 0}] 

220 

221 References 

222 ========== 

223 

224 - Solving differential equations by Symmetry Groups, 

225 John Starrett, pp. 1 - pp. 14 

226 

227 """ 

228 

229 if isinstance(eq, Equality): 

230 eq = eq.lhs - eq.rhs 

231 if not func: 

232 eq, func = _preprocess(eq) 

233 variables = func.args 

234 if len(variables) != 1: 

235 raise ValueError("ODE's have only one independent variable") 

236 else: 

237 x = variables[0] 

238 if not order: 

239 order = ode_order(eq, func) 

240 if order != 1: 

241 raise NotImplementedError("Infinitesimals for only " 

242 "first order ODE's have been implemented") 

243 else: 

244 df = func.diff(x) 

245 # Matching differential equation of the form a*df + b 

246 a = Wild('a', exclude = [df]) 

247 b = Wild('b', exclude = [df]) 

248 if match: # Used by lie_group hint 

249 h = match['h'] 

250 y = match['y'] 

251 else: 

252 match = collect(expand(eq), df).match(a*df + b) 

253 if match: 

254 h = -simplify(match[b]/match[a]) 

255 else: 

256 try: 

257 sol = solve(eq, df) 

258 except NotImplementedError: 

259 raise NotImplementedError("Infinitesimals for the " 

260 "first order ODE could not be found") 

261 else: 

262 h = sol[0] # Find infinitesimals for one solution 

263 y = Dummy("y") 

264 h = h.subs(func, y) 

265 

266 u = Dummy("u") 

267 hx = h.diff(x) 

268 hy = h.diff(y) 

269 hinv = ((1/h).subs([(x, u), (y, x)])).subs(u, y) # Inverse ODE 

270 match = {'h': h, 'func': func, 'hx': hx, 'hy': hy, 'y': y, 'hinv': hinv} 

271 if hint == 'all': 

272 xieta = [] 

273 for heuristic in lie_heuristics: 

274 function = globals()['lie_heuristic_' + heuristic] 

275 inflist = function(match, comp=True) 

276 if inflist: 

277 xieta.extend([inf for inf in inflist if inf not in xieta]) 

278 if xieta: 

279 return xieta 

280 else: 

281 raise NotImplementedError("Infinitesimals could not be found for " 

282 "the given ODE") 

283 

284 elif hint == 'default': 

285 for heuristic in lie_heuristics: 

286 function = globals()['lie_heuristic_' + heuristic] 

287 xieta = function(match, comp=False) 

288 if xieta: 

289 return xieta 

290 

291 raise NotImplementedError("Infinitesimals could not be found for" 

292 " the given ODE") 

293 

294 elif hint not in lie_heuristics: 

295 raise ValueError("Heuristic not recognized: " + hint) 

296 

297 else: 

298 function = globals()['lie_heuristic_' + hint] 

299 xieta = function(match, comp=True) 

300 if xieta: 

301 return xieta 

302 else: 

303 raise ValueError("Infinitesimals could not be found using the" 

304 " given heuristic") 

305 

306 

307def lie_heuristic_abaco1_simple(match, comp=False): 

308 r""" 

309 The first heuristic uses the following four sets of 

310 assumptions on `\xi` and `\eta` 

311 

312 .. math:: \xi = 0, \eta = f(x) 

313 

314 .. math:: \xi = 0, \eta = f(y) 

315 

316 .. math:: \xi = f(x), \eta = 0 

317 

318 .. math:: \xi = f(y), \eta = 0 

319 

320 The success of this heuristic is determined by algebraic factorisation. 

321 For the first assumption `\xi = 0` and `\eta` to be a function of `x`, the PDE 

322 

323 .. math:: \frac{\partial \eta}{\partial x} + (\frac{\partial \eta}{\partial y} 

324 - \frac{\partial \xi}{\partial x})*h 

325 - \frac{\partial \xi}{\partial y}*h^{2} 

326 - \xi*\frac{\partial h}{\partial x} - \eta*\frac{\partial h}{\partial y} = 0 

327 

328 reduces to `f'(x) - f\frac{\partial h}{\partial y} = 0` 

329 If `\frac{\partial h}{\partial y}` is a function of `x`, then this can usually 

330 be integrated easily. A similar idea is applied to the other 3 assumptions as well. 

331 

332 

333 References 

334 ========== 

335 

336 - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra 

337 Solving of First Order ODEs Using Symmetry Methods, pp. 8 

338 

339 

340 """ 

341 

342 xieta = [] 

343 y = match['y'] 

344 h = match['h'] 

345 func = match['func'] 

346 x = func.args[0] 

347 hx = match['hx'] 

348 hy = match['hy'] 

349 xi = Function('xi')(x, func) 

350 eta = Function('eta')(x, func) 

351 

352 hysym = hy.free_symbols 

353 if y not in hysym: 

354 try: 

355 fx = exp(integrate(hy, x)) 

356 except NotImplementedError: 

357 pass 

358 else: 

359 inf = {xi: S.Zero, eta: fx} 

360 if not comp: 

361 return [inf] 

362 if comp and inf not in xieta: 

363 xieta.append(inf) 

364 

365 factor = hy/h 

366 facsym = factor.free_symbols 

367 if x not in facsym: 

368 try: 

369 fy = exp(integrate(factor, y)) 

370 except NotImplementedError: 

371 pass 

372 else: 

373 inf = {xi: S.Zero, eta: fy.subs(y, func)} 

374 if not comp: 

375 return [inf] 

376 if comp and inf not in xieta: 

377 xieta.append(inf) 

378 

379 factor = -hx/h 

380 facsym = factor.free_symbols 

381 if y not in facsym: 

382 try: 

383 fx = exp(integrate(factor, x)) 

384 except NotImplementedError: 

385 pass 

386 else: 

387 inf = {xi: fx, eta: S.Zero} 

388 if not comp: 

389 return [inf] 

390 if comp and inf not in xieta: 

391 xieta.append(inf) 

392 

393 factor = -hx/(h**2) 

394 facsym = factor.free_symbols 

395 if x not in facsym: 

396 try: 

397 fy = exp(integrate(factor, y)) 

398 except NotImplementedError: 

399 pass 

400 else: 

401 inf = {xi: fy.subs(y, func), eta: S.Zero} 

402 if not comp: 

403 return [inf] 

404 if comp and inf not in xieta: 

405 xieta.append(inf) 

406 

407 if xieta: 

408 return xieta 

409 

410def lie_heuristic_abaco1_product(match, comp=False): 

411 r""" 

412 The second heuristic uses the following two assumptions on `\xi` and `\eta` 

413 

414 .. math:: \eta = 0, \xi = f(x)*g(y) 

415 

416 .. math:: \eta = f(x)*g(y), \xi = 0 

417 

418 The first assumption of this heuristic holds good if 

419 `\frac{1}{h^{2}}\frac{\partial^2}{\partial x \partial y}\log(h)` is 

420 separable in `x` and `y`, then the separated factors containing `x` 

421 is `f(x)`, and `g(y)` is obtained by 

422 

423 .. math:: e^{\int f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)\,dy} 

424 

425 provided `f\frac{\partial}{\partial x}\left(\frac{1}{f*h}\right)` is a function 

426 of `y` only. 

427 

428 The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as 

429 `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption 

430 satisfies. After obtaining `f(x)` and `g(y)`, the coordinates are again 

431 interchanged, to get `\eta` as `f(x)*g(y)` 

432 

433 

434 References 

435 ========== 

436 - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order 

437 ODE Patterns, pp. 7 - pp. 8 

438 

439 """ 

440 

441 xieta = [] 

442 y = match['y'] 

443 h = match['h'] 

444 hinv = match['hinv'] 

445 func = match['func'] 

446 x = func.args[0] 

447 xi = Function('xi')(x, func) 

448 eta = Function('eta')(x, func) 

449 

450 

451 inf = separatevars(((log(h).diff(y)).diff(x))/h**2, dict=True, symbols=[x, y]) 

452 if inf and inf['coeff']: 

453 fx = inf[x] 

454 gy = simplify(fx*((1/(fx*h)).diff(x))) 

455 gysyms = gy.free_symbols 

456 if x not in gysyms: 

457 gy = exp(integrate(gy, y)) 

458 inf = {eta: S.Zero, xi: (fx*gy).subs(y, func)} 

459 if not comp: 

460 return [inf] 

461 if comp and inf not in xieta: 

462 xieta.append(inf) 

463 

464 u1 = Dummy("u1") 

465 inf = separatevars(((log(hinv).diff(y)).diff(x))/hinv**2, dict=True, symbols=[x, y]) 

466 if inf and inf['coeff']: 

467 fx = inf[x] 

468 gy = simplify(fx*((1/(fx*hinv)).diff(x))) 

469 gysyms = gy.free_symbols 

470 if x not in gysyms: 

471 gy = exp(integrate(gy, y)) 

472 etaval = fx*gy 

473 etaval = (etaval.subs([(x, u1), (y, x)])).subs(u1, y) 

474 inf = {eta: etaval.subs(y, func), xi: S.Zero} 

475 if not comp: 

476 return [inf] 

477 if comp and inf not in xieta: 

478 xieta.append(inf) 

479 

480 if xieta: 

481 return xieta 

482 

483def lie_heuristic_bivariate(match, comp=False): 

484 r""" 

485 The third heuristic assumes the infinitesimals `\xi` and `\eta` 

486 to be bi-variate polynomials in `x` and `y`. The assumption made here 

487 for the logic below is that `h` is a rational function in `x` and `y` 

488 though that may not be necessary for the infinitesimals to be 

489 bivariate polynomials. The coefficients of the infinitesimals 

490 are found out by substituting them in the PDE and grouping similar terms 

491 that are polynomials and since they form a linear system, solve and check 

492 for non trivial solutions. The degree of the assumed bivariates 

493 are increased till a certain maximum value. 

494 

495 References 

496 ========== 

497 - Lie Groups and Differential Equations 

498 pp. 327 - pp. 329 

499 

500 """ 

501 

502 h = match['h'] 

503 hx = match['hx'] 

504 hy = match['hy'] 

505 func = match['func'] 

506 x = func.args[0] 

507 y = match['y'] 

508 xi = Function('xi')(x, func) 

509 eta = Function('eta')(x, func) 

510 

511 if h.is_rational_function(): 

512 # The maximum degree that the infinitesimals can take is 

513 # calculated by this technique. 

514 etax, etay, etad, xix, xiy, xid = symbols("etax etay etad xix xiy xid") 

515 ipde = etax + (etay - xix)*h - xiy*h**2 - xid*hx - etad*hy 

516 num, denom = cancel(ipde).as_numer_denom() 

517 deg = Poly(num, x, y).total_degree() 

518 deta = Function('deta')(x, y) 

519 dxi = Function('dxi')(x, y) 

520 ipde = (deta.diff(x) + (deta.diff(y) - dxi.diff(x))*h - (dxi.diff(y))*h**2 

521 - dxi*hx - deta*hy) 

522 xieq = Symbol("xi0") 

523 etaeq = Symbol("eta0") 

524 

525 for i in range(deg + 1): 

526 if i: 

527 xieq += Add(*[ 

528 Symbol("xi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) 

529 for power in range(i + 1)]) 

530 etaeq += Add(*[ 

531 Symbol("eta_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) 

532 for power in range(i + 1)]) 

533 pden, denom = (ipde.subs({dxi: xieq, deta: etaeq}).doit()).as_numer_denom() 

534 pden = expand(pden) 

535 

536 # If the individual terms are monomials, the coefficients 

537 # are grouped 

538 if pden.is_polynomial(x, y) and pden.is_Add: 

539 polyy = Poly(pden, x, y).as_dict() 

540 if polyy: 

541 symset = xieq.free_symbols.union(etaeq.free_symbols) - {x, y} 

542 soldict = solve(polyy.values(), *symset) 

543 if isinstance(soldict, list): 

544 soldict = soldict[0] 

545 if any(soldict.values()): 

546 xired = xieq.subs(soldict) 

547 etared = etaeq.subs(soldict) 

548 # Scaling is done by substituting one for the parameters 

549 # This can be any number except zero. 

550 dict_ = {sym: 1 for sym in symset} 

551 inf = {eta: etared.subs(dict_).subs(y, func), 

552 xi: xired.subs(dict_).subs(y, func)} 

553 return [inf] 

554 

555def lie_heuristic_chi(match, comp=False): 

556 r""" 

557 The aim of the fourth heuristic is to find the function `\chi(x, y)` 

558 that satisfies the PDE `\frac{d\chi}{dx} + h\frac{d\chi}{dx} 

559 - \frac{\partial h}{\partial y}\chi = 0`. 

560 

561 This assumes `\chi` to be a bivariate polynomial in `x` and `y`. By intuition, 

562 `h` should be a rational function in `x` and `y`. The method used here is 

563 to substitute a general binomial for `\chi` up to a certain maximum degree 

564 is reached. The coefficients of the polynomials, are calculated by by collecting 

565 terms of the same order in `x` and `y`. 

566 

567 After finding `\chi`, the next step is to use `\eta = \xi*h + \chi`, to 

568 determine `\xi` and `\eta`. This can be done by dividing `\chi` by `h` 

569 which would give `-\xi` as the quotient and `\eta` as the remainder. 

570 

571 

572 References 

573 ========== 

574 - E.S Cheb-Terrab, L.G.S Duarte and L.A,C.P da Mota, Computer Algebra 

575 Solving of First Order ODEs Using Symmetry Methods, pp. 8 

576 

577 """ 

578 

579 h = match['h'] 

580 hy = match['hy'] 

581 func = match['func'] 

582 x = func.args[0] 

583 y = match['y'] 

584 xi = Function('xi')(x, func) 

585 eta = Function('eta')(x, func) 

586 

587 if h.is_rational_function(): 

588 schi, schix, schiy = symbols("schi, schix, schiy") 

589 cpde = schix + h*schiy - hy*schi 

590 num, denom = cancel(cpde).as_numer_denom() 

591 deg = Poly(num, x, y).total_degree() 

592 

593 chi = Function('chi')(x, y) 

594 chix = chi.diff(x) 

595 chiy = chi.diff(y) 

596 cpde = chix + h*chiy - hy*chi 

597 chieq = Symbol("chi") 

598 for i in range(1, deg + 1): 

599 chieq += Add(*[ 

600 Symbol("chi_" + str(power) + "_" + str(i - power))*x**power*y**(i - power) 

601 for power in range(i + 1)]) 

602 cnum, cden = cancel(cpde.subs({chi : chieq}).doit()).as_numer_denom() 

603 cnum = expand(cnum) 

604 if cnum.is_polynomial(x, y) and cnum.is_Add: 

605 cpoly = Poly(cnum, x, y).as_dict() 

606 if cpoly: 

607 solsyms = chieq.free_symbols - {x, y} 

608 soldict = solve(cpoly.values(), *solsyms) 

609 if isinstance(soldict, list): 

610 soldict = soldict[0] 

611 if any(soldict.values()): 

612 chieq = chieq.subs(soldict) 

613 dict_ = {sym: 1 for sym in solsyms} 

614 chieq = chieq.subs(dict_) 

615 # After finding chi, the main aim is to find out 

616 # eta, xi by the equation eta = xi*h + chi 

617 # One method to set xi, would be rearranging it to 

618 # (eta/h) - xi = (chi/h). This would mean dividing 

619 # chi by h would give -xi as the quotient and eta 

620 # as the remainder. Thanks to Sean Vig for suggesting 

621 # this method. 

622 xic, etac = div(chieq, h) 

623 inf = {eta: etac.subs(y, func), xi: -xic.subs(y, func)} 

624 return [inf] 

625 

626def lie_heuristic_function_sum(match, comp=False): 

627 r""" 

628 This heuristic uses the following two assumptions on `\xi` and `\eta` 

629 

630 .. math:: \eta = 0, \xi = f(x) + g(y) 

631 

632 .. math:: \eta = f(x) + g(y), \xi = 0 

633 

634 The first assumption of this heuristic holds good if 

635 

636 .. math:: \frac{\partial}{\partial y}[(h\frac{\partial^{2}}{ 

637 \partial x^{2}}(h^{-1}))^{-1}] 

638 

639 is separable in `x` and `y`, 

640 

641 1. The separated factors containing `y` is `\frac{\partial g}{\partial y}`. 

642 From this `g(y)` can be determined. 

643 2. The separated factors containing `x` is `f''(x)`. 

644 3. `h\frac{\partial^{2}}{\partial x^{2}}(h^{-1})` equals 

645 `\frac{f''(x)}{f(x) + g(y)}`. From this `f(x)` can be determined. 

646 

647 The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as 

648 `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first 

649 assumption satisfies. After obtaining `f(x)` and `g(y)`, the coordinates 

650 are again interchanged, to get `\eta` as `f(x) + g(y)`. 

651 

652 For both assumptions, the constant factors are separated among `g(y)` 

653 and `f''(x)`, such that `f''(x)` obtained from 3] is the same as that 

654 obtained from 2]. If not possible, then this heuristic fails. 

655 

656 

657 References 

658 ========== 

659 - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order 

660 ODE Patterns, pp. 7 - pp. 8 

661 

662 """ 

663 

664 xieta = [] 

665 h = match['h'] 

666 func = match['func'] 

667 hinv = match['hinv'] 

668 x = func.args[0] 

669 y = match['y'] 

670 xi = Function('xi')(x, func) 

671 eta = Function('eta')(x, func) 

672 

673 for odefac in [h, hinv]: 

674 factor = odefac*((1/odefac).diff(x, 2)) 

675 sep = separatevars((1/factor).diff(y), dict=True, symbols=[x, y]) 

676 if sep and sep['coeff'] and sep[x].has(x) and sep[y].has(y): 

677 k = Dummy("k") 

678 try: 

679 gy = k*integrate(sep[y], y) 

680 except NotImplementedError: 

681 pass 

682 else: 

683 fdd = 1/(k*sep[x]*sep['coeff']) 

684 fx = simplify(fdd/factor - gy) 

685 check = simplify(fx.diff(x, 2) - fdd) 

686 if fx: 

687 if not check: 

688 fx = fx.subs(k, 1) 

689 gy = (gy/k) 

690 else: 

691 sol = solve(check, k) 

692 if sol: 

693 sol = sol[0] 

694 fx = fx.subs(k, sol) 

695 gy = (gy/k)*sol 

696 else: 

697 continue 

698 if odefac == hinv: # Inverse ODE 

699 fx = fx.subs(x, y) 

700 gy = gy.subs(y, x) 

701 etaval = factor_terms(fx + gy) 

702 if etaval.is_Mul: 

703 etaval = Mul(*[arg for arg in etaval.args if arg.has(x, y)]) 

704 if odefac == hinv: # Inverse ODE 

705 inf = {eta: etaval.subs(y, func), xi : S.Zero} 

706 else: 

707 inf = {xi: etaval.subs(y, func), eta : S.Zero} 

708 if not comp: 

709 return [inf] 

710 else: 

711 xieta.append(inf) 

712 

713 if xieta: 

714 return xieta 

715 

716def lie_heuristic_abaco2_similar(match, comp=False): 

717 r""" 

718 This heuristic uses the following two assumptions on `\xi` and `\eta` 

719 

720 .. math:: \eta = g(x), \xi = f(x) 

721 

722 .. math:: \eta = f(y), \xi = g(y) 

723 

724 For the first assumption, 

725 

726 1. First `\frac{\frac{\partial h}{\partial y}}{\frac{\partial^{2} h}{ 

727 \partial yy}}` is calculated. Let us say this value is A 

728 

729 2. If this is constant, then `h` is matched to the form `A(x) + B(x)e^{ 

730 \frac{y}{C}}` then, `\frac{e^{\int \frac{A(x)}{C} \,dx}}{B(x)}` gives `f(x)` 

731 and `A(x)*f(x)` gives `g(x)` 

732 

733 3. Otherwise `\frac{\frac{\partial A}{\partial X}}{\frac{\partial A}{ 

734 \partial Y}} = \gamma` is calculated. If 

735 

736 a] `\gamma` is a function of `x` alone 

737 

738 b] `\frac{\gamma\frac{\partial h}{\partial y} - \gamma'(x) - \frac{ 

739 \partial h}{\partial x}}{h + \gamma} = G` is a function of `x` alone. 

740 then, `e^{\int G \,dx}` gives `f(x)` and `-\gamma*f(x)` gives `g(x)` 

741 

742 The second assumption holds good if `\frac{dy}{dx} = h(x, y)` is rewritten as 

743 `\frac{dy}{dx} = \frac{1}{h(y, x)}` and the same properties of the first assumption 

744 satisfies. After obtaining `f(x)` and `g(x)`, the coordinates are again 

745 interchanged, to get `\xi` as `f(x^*)` and `\eta` as `g(y^*)` 

746 

747 References 

748 ========== 

749 - E.S. Cheb-Terrab, A.D. Roche, Symmetries and First Order 

750 ODE Patterns, pp. 10 - pp. 12 

751 

752 """ 

753 

754 h = match['h'] 

755 hx = match['hx'] 

756 hy = match['hy'] 

757 func = match['func'] 

758 hinv = match['hinv'] 

759 x = func.args[0] 

760 y = match['y'] 

761 xi = Function('xi')(x, func) 

762 eta = Function('eta')(x, func) 

763 

764 factor = cancel(h.diff(y)/h.diff(y, 2)) 

765 factorx = factor.diff(x) 

766 factory = factor.diff(y) 

767 if not factor.has(x) and not factor.has(y): 

768 A = Wild('A', exclude=[y]) 

769 B = Wild('B', exclude=[y]) 

770 C = Wild('C', exclude=[x, y]) 

771 match = h.match(A + B*exp(y/C)) 

772 try: 

773 tau = exp(-integrate(match[A]/match[C]), x)/match[B] 

774 except NotImplementedError: 

775 pass 

776 else: 

777 gx = match[A]*tau 

778 return [{xi: tau, eta: gx}] 

779 

780 else: 

781 gamma = cancel(factorx/factory) 

782 if not gamma.has(y): 

783 tauint = cancel((gamma*hy - gamma.diff(x) - hx)/(h + gamma)) 

784 if not tauint.has(y): 

785 try: 

786 tau = exp(integrate(tauint, x))