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

1r""" 

2This module is intended for solving recurrences or, in other words, 

3difference equations. Currently supported are linear, inhomogeneous 

4equations with polynomial or rational coefficients. 

5 

6The solutions are obtained among polynomials, rational functions, 

7hypergeometric terms, or combinations of hypergeometric term which 

8are pairwise dissimilar. 

9 

10``rsolve_X`` functions were meant as a low level interface 

11for ``rsolve`` which would use Mathematica's syntax. 

12 

13Given a recurrence relation: 

14 

15 .. math:: a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) + 

16 ... + a_{0}(n) y(n) = f(n) 

17 

18where `k > 0` and `a_{i}(n)` are polynomials in `n`. To use 

19``rsolve_X`` we need to put all coefficients in to a list ``L`` of 

20`k+1` elements the following way: 

21 

22 ``L = [a_{0}(n), ..., a_{k-1}(n), a_{k}(n)]`` 

23 

24where ``L[i]``, for `i=0, \ldots, k`, maps to 

25`a_{i}(n) y(n+i)` (`y(n+i)` is implicit). 

26 

27For example if we would like to compute `m`-th Bernoulli polynomial 

28up to a constant (example was taken from rsolve_poly docstring), 

29then we would use `b(n+1) - b(n) = m n^{m-1}` recurrence, which 

30has solution `b(n) = B_m + C`. 

31 

32Then ``L = [-1, 1]`` and `f(n) = m n^(m-1)` and finally for `m=4`: 

33 

34>>> from sympy import Symbol, bernoulli, rsolve_poly 

35>>> n = Symbol('n', integer=True) 

36 

37>>> rsolve_poly([-1, 1], 4*n**3, n) 

38C0 + n**4 - 2*n**3 + n**2 

39 

40>>> bernoulli(4, n) 

41n**4 - 2*n**3 + n**2 - 1/30 

42 

43For the sake of completeness, `f(n)` can be: 

44 

45 [1] a polynomial -> rsolve_poly 

46 [2] a rational function -> rsolve_ratio 

47 [3] a hypergeometric function -> rsolve_hyper 

48""" 

49from collections import defaultdict 

50 

51from sympy.concrete import product 

52from sympy.core.singleton import S 

53from sympy.core.numbers import Rational, I 

54from sympy.core.symbol import Symbol, Wild, Dummy 

55from sympy.core.relational import Equality 

56from sympy.core.add import Add 

57from sympy.core.mul import Mul 

58from sympy.core.sorting import default_sort_key 

59from sympy.core.sympify import sympify 

60 

61from sympy.simplify import simplify, hypersimp, hypersimilar # type: ignore 

62from sympy.solvers import solve, solve_undetermined_coeffs 

63from sympy.polys import Poly, quo, gcd, lcm, roots, resultant 

64from sympy.functions import binomial, factorial, FallingFactorial, RisingFactorial 

65from sympy.matrices import Matrix, casoratian 

66from sympy.utilities.iterables import numbered_symbols 

67 

68 

69def rsolve_poly(coeffs, f, n, shift=0, **hints): 

70 r""" 

71 Given linear recurrence operator `\operatorname{L}` of order 

72 `k` with polynomial coefficients and inhomogeneous equation 

73 `\operatorname{L} y = f`, where `f` is a polynomial, we seek for 

74 all polynomial solutions over field `K` of characteristic zero. 

75 

76 The algorithm performs two basic steps: 

77 

78 (1) Compute degree `N` of the general polynomial solution. 

79 (2) Find all polynomials of degree `N` or less 

80 of `\operatorname{L} y = f`. 

81 

82 There are two methods for computing the polynomial solutions. 

83 If the degree bound is relatively small, i.e. it's smaller than 

84 or equal to the order of the recurrence, then naive method of 

85 undetermined coefficients is being used. This gives a system 

86 of algebraic equations with `N+1` unknowns. 

87 

88 In the other case, the algorithm performs transformation of the 

89 initial equation to an equivalent one for which the system of 

90 algebraic equations has only `r` indeterminates. This method is 

91 quite sophisticated (in comparison with the naive one) and was 

92 invented together by Abramov, Bronstein and Petkovsek. 

93 

94 It is possible to generalize the algorithm implemented here to 

95 the case of linear q-difference and differential equations. 

96 

97 Lets say that we would like to compute `m`-th Bernoulli polynomial 

98 up to a constant. For this we can use `b(n+1) - b(n) = m n^{m-1}` 

99 recurrence, which has solution `b(n) = B_m + C`. For example: 

100 

101 >>> from sympy import Symbol, rsolve_poly 

102 >>> n = Symbol('n', integer=True) 

103 

104 >>> rsolve_poly([-1, 1], 4*n**3, n) 

105 C0 + n**4 - 2*n**3 + n**2 

106 

107 References 

108 ========== 

109 

110 .. [1] S. A. Abramov, M. Bronstein and M. Petkovsek, On polynomial 

111 solutions of linear operator equations, in: T. Levelt, ed., 

112 Proc. ISSAC '95, ACM Press, New York, 1995, 290-296. 

113 

114 .. [2] M. Petkovsek, Hypergeometric solutions of linear recurrences 

115 with polynomial coefficients, J. Symbolic Computation, 

116 14 (1992), 243-264. 

117 

118 .. [3] M. Petkovsek, H. S. Wilf, D. Zeilberger, A = B, 1996. 

119 

120 """ 

121 f = sympify(f) 

122 

123 if not f.is_polynomial(n): 

124 return None 

125 

126 homogeneous = f.is_zero 

127 

128 r = len(coeffs) - 1 

129 

130 coeffs = [Poly(coeff, n) for coeff in coeffs] 

131 

132 polys = [Poly(0, n)]*(r + 1) 

133 terms = [(S.Zero, S.NegativeInfinity)]*(r + 1) 

134 

135 for i in range(r + 1): 

136 for j in range(i, r + 1): 

137 polys[i] += coeffs[j]*(binomial(j, i).as_poly(n)) 

138 

139 if not polys[i].is_zero: 

140 (exp,), coeff = polys[i].LT() 

141 terms[i] = (coeff, exp) 

142 

143 d = b = terms[0][1] 

144 

145 for i in range(1, r + 1): 

146 if terms[i][1] > d: 

147 d = terms[i][1] 

148 

149 if terms[i][1] - i > b: 

150 b = terms[i][1] - i 

151 

152 d, b = int(d), int(b) 

153 

154 x = Dummy('x') 

155 

156 degree_poly = S.Zero 

157 

158 for i in range(r + 1): 

159 if terms[i][1] - i == b: 

160 degree_poly += terms[i][0]*FallingFactorial(x, i) 

161 

162 nni_roots = list(roots(degree_poly, x, filter='Z', 

163 predicate=lambda r: r >= 0).keys()) 

164 

165 if nni_roots: 

166 N = [max(nni_roots)] 

167 else: 

168 N = [] 

169 

170 if homogeneous: 

171 N += [-b - 1] 

172 else: 

173 N += [f.as_poly(n).degree() - b, -b - 1] 

174 

175 N = int(max(N)) 

176 

177 if N < 0: 

178 if homogeneous: 

179 if hints.get('symbols', False): 

180 return (S.Zero, []) 

181 else: 

182 return S.Zero 

183 else: 

184 return None 

185 

186 if N <= r: 

187 C = [] 

188 y = E = S.Zero 

189 

190 for i in range(N + 1): 

191 C.append(Symbol('C' + str(i + shift))) 

192 y += C[i] * n**i 

193 

194 for i in range(r + 1): 

195 E += coeffs[i].as_expr()*y.subs(n, n + i) 

196 

197 solutions = solve_undetermined_coeffs(E - f, C, n) 

198 

199 if solutions is not None: 

200 _C = C 

201 C = [c for c in C if (c not in solutions)] 

202 result = y.subs(solutions) 

203 else: 

204 return None # TBD 

205 else: 

206 A = r 

207 U = N + A + b + 1 

208 

209 nni_roots = list(roots(polys[r], filter='Z', 

210 predicate=lambda r: r >= 0).keys()) 

211 

212 if nni_roots != []: 

213 a = max(nni_roots) + 1 

214 else: 

215 a = S.Zero 

216 

217 def _zero_vector(k): 

218 return [S.Zero] * k 

219 

220 def _one_vector(k): 

221 return [S.One] * k 

222 

223 def _delta(p, k): 

224 B = S.One 

225 D = p.subs(n, a + k) 

226 

227 for i in range(1, k + 1): 

228 B *= Rational(i - k - 1, i) 

229 D += B * p.subs(n, a + k - i) 

230 

231 return D 

232 

233 alpha = {} 

234 

235 for i in range(-A, d + 1): 

236 I = _one_vector(d + 1) 

237 

238 for k in range(1, d + 1): 

239 I[k] = I[k - 1] * (x + i - k + 1)/k 

240 

241 alpha[i] = S.Zero 

242 

243 for j in range(A + 1): 

244 for k in range(d + 1): 

245 B = binomial(k, i + j) 

246 D = _delta(polys[j].as_expr(), k) 

247 

248 alpha[i] += I[k]*B*D 

249 

250 V = Matrix(U, A, lambda i, j: int(i == j)) 

251 

252 if homogeneous: 

253 for i in range(A, U): 

254 v = _zero_vector(A) 

255 

256 for k in range(1, A + b + 1): 

257 if i - k < 0: 

258 break 

259 

260 B = alpha[k - A].subs(x, i - k) 

261 

262 for j in range(A): 

263 v[j] += B * V[i - k, j] 

264 

265 denom = alpha[-A].subs(x, i) 

266 

267 for j in range(A): 

268 V[i, j] = -v[j] / denom 

269 else: 

270 G = _zero_vector(U) 

271 

272 for i in range(A, U): 

273 v = _zero_vector(A) 

274 g = S.Zero 

275 

276 for k in range(1, A + b + 1): 

277 if i - k < 0: 

278 break 

279 

280 B = alpha[k - A].subs(x, i - k) 

281 

282 for j in range(A): 

283 v[j] += B * V[i - k, j] 

284 

285 g += B * G[i - k] 

286 

287 denom = alpha[-A].subs(x, i) 

288 

289 for j in range(A): 

290 V[i, j] = -v[j] / denom 

291 

292 G[i] = (_delta(f, i - A) - g) / denom 

293 

294 P, Q = _one_vector(U), _zero_vector(A) 

295 

296 for i in range(1, U): 

297 P[i] = (P[i - 1] * (n - a - i + 1)/i).expand() 

298 

299 for i in range(A): 

300 Q[i] = Add(*[(v*p).expand() for v, p in zip(V[:, i], P)]) 

301 

302 if not homogeneous: 

303 h = Add(*[(g*p).expand() for g, p in zip(G, P)]) 

304 

305 C = [Symbol('C' + str(i + shift)) for i in range(A)] 

306 

307 g = lambda i: Add(*[c*_delta(q, i) for c, q in zip(C, Q)]) 

308 

309 if homogeneous: 

310 E = [g(i) for i in range(N + 1, U)] 

311 else: 

312 E = [g(i) + _delta(h, i) for i in range(N + 1, U)] 

313 

314 if E != []: 

315 solutions = solve(E, *C) 

316 

317 if not solutions: 

318 if homogeneous: 

319 if hints.get('symbols', False): 

320 return (S.Zero, []) 

321 else: 

322 return S.Zero 

323 else: 

324 return None 

325 else: 

326 solutions = {} 

327 

328 if homogeneous: 

329 result = S.Zero 

330 else: 

331 result = h 

332 

333 _C = C[:] 

334 for c, q in list(zip(C, Q)): 

335 if c in solutions: 

336 s = solutions[c]*q 

337 C.remove(c) 

338 else: 

339 s = c*q 

340 

341 result += s.expand() 

342 

343 if C != _C: 

344 # renumber so they are contiguous 

345 result = result.xreplace(dict(zip(C, _C))) 

346 C = _C[:len(C)] 

347 

348 if hints.get('symbols', False): 

349 return (result, C) 

350 else: 

351 return result 

352 

353 

354def rsolve_ratio(coeffs, f, n, **hints): 

355 r""" 

356 Given linear recurrence operator `\operatorname{L}` of order `k` 

357 with polynomial coefficients and inhomogeneous equation 

358 `\operatorname{L} y = f`, where `f` is a polynomial, we seek 

359 for all rational solutions over field `K` of characteristic zero. 

360 

361 This procedure accepts only polynomials, however if you are 

362 interested in solving recurrence with rational coefficients 

363 then use ``rsolve`` which will pre-process the given equation 

364 and run this procedure with polynomial arguments. 

365 

366 The algorithm performs two basic steps: 

367 

368 (1) Compute polynomial `v(n)` which can be used as universal 

369 denominator of any rational solution of equation 

370 `\operatorname{L} y = f`. 

371 

372 (2) Construct new linear difference equation by substitution 

373 `y(n) = u(n)/v(n)` and solve it for `u(n)` finding all its 

374 polynomial solutions. Return ``None`` if none were found. 

375 

376 The algorithm implemented here is a revised version of the original 

377 Abramov's algorithm, developed in 1989. The new approach is much 

378 simpler to implement and has better overall efficiency. This 

379 method can be easily adapted to the q-difference equations case. 

380 

381 Besides finding rational solutions alone, this functions is 

382 an important part of Hyper algorithm where it is used to find 

383 a particular solution for the inhomogeneous part of a recurrence. 

384 

385 Examples 

386 ======== 

387 

388 >>> from sympy.abc import x 

389 >>> from sympy.solvers.recurr import rsolve_ratio 

390 >>> rsolve_ratio([-2*x**3 + x**2 + 2*x - 1, 2*x**3 + x**2 - 6*x, 

391 ... - 2*x**3 - 11*x**2 - 18*x - 9, 2*x**3 + 13*x**2 + 22*x + 8], 0, x) 

392 C0*(2*x - 3)/(2*(x**2 - 1)) 

393 

394 References 

395 ========== 

396 

397 .. [1] S. A. Abramov, Rational solutions of linear difference 

398 and q-difference equations with polynomial coefficients, 

399 in: T. Levelt, ed., Proc. ISSAC '95, ACM Press, New York, 

400 1995, 285-289 

401 

402 See Also 

403 ======== 

404 

405 rsolve_hyper 

406 """ 

407 f = sympify(f) 

408 

409 if not f.is_polynomial(n): 

410 return None 

411 

412 coeffs = list(map(sympify, coeffs)) 

413 

414 r = len(coeffs) - 1 

415 

416 A, B = coeffs[r], coeffs[0] 

417 A = A.subs(n, n - r).expand() 

418 

419 h = Dummy('h') 

420 

421 res = resultant(A, B.subs(n, n + h), n) 

422 

423 if not res.is_polynomial(h): 

424 p, q = res.as_numer_denom() 

425 res = quo(p, q, h) 

426 

427 nni_roots = list(roots(res, h, filter='Z', 

428 predicate=lambda r: r >= 0).keys()) 

429 

430 if not nni_roots: 

431 return rsolve_poly(coeffs, f, n, **hints) 

432 else: 

433 C, numers = S.One, [S.Zero]*(r + 1) 

434 

435 for i in range(int(max(nni_roots)), -1, -1): 

436 d = gcd(A, B.subs(n, n + i), n) 

437 

438 A = quo(A, d, n) 

439 B = quo(B, d.subs(n, n - i), n) 

440 

441 C *= Mul(*[d.subs(n, n - j) for j in range(i + 1)]) 

442 

443 denoms = [C.subs(n, n + i) for i in range(r + 1)] 

444 

445 for i in range(r + 1): 

446 g = gcd(coeffs[i], denoms[i], n) 

447 

448 numers[i] = quo(coeffs[i], g, n) 

449 denoms[i] = quo(denoms[i], g, n) 

450 

451 for i in range(r + 1): 

452 numers[i] *= Mul(*(denoms[:i] + denoms[i + 1:])) 

453 

454 result = rsolve_poly(numers, f * Mul(*denoms), n, **hints) 

455 

456 if result is not None: 

457 if hints.get('symbols', False): 

458 return (simplify(result[0] / C), result[1]) 

459 else: 

460 return simplify(result / C) 

461 else: 

462 return None 

463 

464 

465def rsolve_hyper(coeffs, f, n, **hints): 

466 r""" 

467 Given linear recurrence operator `\operatorname{L}` of order `k` 

468 with polynomial coefficients and inhomogeneous equation 

469 `\operatorname{L} y = f` we seek for all hypergeometric solutions 

470 over field `K` of characteristic zero. 

471 

472 The inhomogeneous part can be either hypergeometric or a sum 

473 of a fixed number of pairwise dissimilar hypergeometric terms. 

474 

475 The algorithm performs three basic steps: 

476 

477 (1) Group together similar hypergeometric terms in the 

478 inhomogeneous part of `\operatorname{L} y = f`, and find 

479 particular solution using Abramov's algorithm. 

480 

481 (2) Compute generating set of `\operatorname{L}` and find basis 

482 in it, so that all solutions are linearly independent. 

483 

484 (3) Form final solution with the number of arbitrary 

485 constants equal to dimension of basis of `\operatorname{L}`. 

486 

487 Term `a(n)` is hypergeometric if it is annihilated by first order 

488 linear difference equations with polynomial coefficients or, in 

489 simpler words, if consecutive term ratio is a rational function. 

490 

491 The output of this procedure is a linear combination of fixed 

492 number of hypergeometric terms. However the underlying method 

493 can generate larger class of solutions - D'Alembertian terms. 

494 

495 Note also that this method not only computes the kernel of the 

496 inhomogeneous equation, but also reduces in to a basis so that 

497 solutions generated by this procedure are linearly independent 

498 

499 Examples 

500 ======== 

501 

502 >>> from sympy.solvers import rsolve_hyper 

503 >>> from sympy.abc import x 

504 

505 >>> rsolve_hyper([-1, -1, 1], 0, x) 

506 C0*(1/2 - sqrt(5)/2)**x + C1*(1/2 + sqrt(5)/2)**x 

507 

508 >>> rsolve_hyper([-1, 1], 1 + x, x) 

509 C0 + x*(x + 1)/2 

510 

511 References 

512 ========== 

513 

514 .. [1] M. Petkovsek, Hypergeometric solutions of linear recurrences 

515 with polynomial coefficients, J. Symbolic Computation, 

516 14 (1992), 243-264. 

517 

518 .. [2] M. Petkovsek, H. S. Wilf, D. Zeilberger, A = B, 1996. 

519 """ 

520 coeffs = list(map(sympify, coeffs)) 

521 

522 f = sympify(f) 

523 

524 r, kernel, symbols = len(coeffs) - 1, [], set() 

525 

526 if not f.is_zero: 

527 if f.is_Add: 

528 similar = {} 

529 

530 for g in f.expand().args: 

531 if not g.is_hypergeometric(n): 

532 return None 

533 

534 for h in similar.keys(): 

535 if hypersimilar(g, h, n): 

536 similar[h] += g 

537 break 

538 else: 

539 similar[g] = S.Zero 

540 

541 inhomogeneous = [g + h for g, h in similar.items()] 

542 elif f.is_hypergeometric(n): 

543 inhomogeneous = [f] 

544 else: 

545 return None 

546 

547 for i, g in enumerate(inhomogeneous): 

548 coeff, polys = S.One, coeffs[:] 

549 denoms = [S.One]*(r + 1) 

550 

551 s = hypersimp(g, n) 

552 

553 for j in range(1, r + 1): 

554 coeff *= s.subs(n, n + j - 1) 

555 

556 p, q = coeff.as_numer_denom() 

557 

558 polys[j] *= p 

559 denoms[j] = q 

560 

561 for j in range(r + 1): 

562 polys[j] *= Mul(*(denoms[:j] + denoms[j + 1:])) 

563 

564 # FIXME: The call to rsolve_ratio below should suffice (rsolve_poly 

565 # call can be removed) but the XFAIL test_rsolve_ratio_missed must 

566 # be fixed first. 

567 R = rsolve_ratio(polys, Mul(*denoms), n, symbols=True) 

568 if R is not None: 

569 R, syms = R 

570 if syms: 

571 R = R.subs(zip(syms, [0]*len(syms))) 

572 else: 

573 R = rsolve_poly(polys, Mul(*denoms), n) 

574 

575 if R: 

576 inhomogeneous[i] *= R 

577 else: 

578 return None 

579 

580 result = Add(*inhomogeneous) 

581 result = simplify(result) 

582 else: 

583 result = S.Zero 

584 

585 Z = Dummy('Z') 

586 

587 p, q = coeffs[0], coeffs[r].subs(n, n - r + 1) 

588 

589 p_factors = list(roots(p, n).keys()) 

590 q_factors = list(roots(q, n).keys()) 

591 

592 factors = [(S.One, S.One)] 

593 

594 for p in p_factors: 

595 for q in q_factors: 

596 if p.is_integer and q.is_integer and p <= q: 

597 continue 

598 else: 

599 factors += [(n - p, n - q)] 

600 

601 p = [(n - p, S.One) for p in p_factors] 

602 q = [(S.One, n - q) for q in q_factors] 

603 

604 factors = p + factors + q 

605 

606 for A, B in factors: 

607 polys, degrees = [], [] 

608 D = A*B.subs(n, n + r - 1) 

609 

610 for i in range(r + 1): 

611 a = Mul(*[A.subs(n, n + j) for j in range(i)]) 

612 b = Mul(*[B.subs(n, n + j) for j in range(i, r)]) 

613 

614 poly = quo(coeffs[i]*a*b, D, n) 

615 polys.append(poly.as_poly(n)) 

616 

617 if not poly.is_zero: 

618 degrees.append(polys[i].degree()) 

619 

620 if degrees: 

621 d, poly = max(degrees), S.Zero 

622 else: 

623 return None 

624 

625 for i in range(r + 1): 

626 coeff = polys[i].nth(d) 

627 

628 if coeff is not S.Zero: 

629 poly += coeff * Z**i 

630 

631 for z in roots(poly, Z).keys(): 

632 if z.is_zero: 

633 continue 

634 

635 recurr_coeffs = [polys[i].as_expr()*z**i for i in range(r + 1)] 

636 if d == 0 and 0 != Add(*[recurr_coeffs[j]*j for j in range(1, r + 1)]): 

637 # faster inline check (than calling rsolve_poly) for a 

638 # constant solution to a constant coefficient recurrence. 

639 sol = [Symbol("C" + str(len(symbols)))] 

640 else: 

641 sol, syms = rsolve_poly(recurr_coeffs, 0, n, len(symbols), symbols=True) 

642 sol = sol.collect(syms) 

643 sol = [sol.coeff(s) for s in syms] 

644 

645 for C in sol: 

646 ratio = z * A * C.subs(n, n + 1) / B / C 

647 ratio = simplify(ratio) 

648 # If there is a nonnegative root in the denominator of the ratio, 

649 # this indicates that the term y(n_root) is zero, and one should 

650 # start the product with the term y(n_root + 1). 

651 n0 = 0 

652 for n_root in roots(ratio.as_numer_denom()[1], n).keys(): 

653 if n_root.has(I): 

654 return None 

655 elif (n0 < (n_root + 1)) == True: 

656 n0 = n_root + 1 

657 K = product(ratio, (n, n0, n - 1)) 

658 if K.has(factorial, FallingFactorial, RisingFactorial): 

659 K = simplify(K) 

660 

661 if casoratian(kernel + [K], n, zero=False) != 0: 

662 kernel.append(K) 

663 

664 kernel.sort(key=default_sort_key) 

665 sk = list(zip(numbered_symbols('C'), kernel)) 

666 

667 for C, ker in sk: 

668 result += C * ker 

669 

670 if hints.get('symbols', False): 

671 # XXX: This returns the symbols in a non-deterministic order 

672 symbols |= {s for s, k in sk} 

673 return (result, list(symbols)) 

674 else: 

675 return result 

676 

677 

678def rsolve(f, y, init=None): 

679 r""" 

680 Solve univariate recurrence with rational coefficients. 

681 

682 Given `k`-th order linear recurrence `\operatorname{L} y = f`, 

683 or equivalently: 

684 

685 .. math:: a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) + 

686 \cdots + a_{0}(n) y(n) = f(n) 

687 

688 where `a_{i}(n)`, for `i=0, \ldots, k`, are polynomials or rational 

689 functions in `n`, and `f` is a hypergeometric function or a sum 

690 of a fixed number of pairwise dissimilar hypergeometric terms in 

691 `n`, finds all solutions or returns ``None``, if none were found. 

692 

693 Initial conditions can be given as a dictionary in two forms: 

694 

695 (1) ``{ n_0 : v_0, n_1 : v_1, ..., n_m : v_m}`` 

696 (2) ``{y(n_0) : v_0, y(n_1) : v_1, ..., y(n_m) : v_m}`` 

697 

698 or as a list ``L`` of values: 

699 

700 ``L = [v_0, v_1, ..., v_m]`` 

701 

702 where ``L[i] = v_i``, for `i=0, \ldots, m`, maps to `y(n_i)`. 

703 

704 Examples 

705 ======== 

706 

707 Lets consider the following recurrence: 

708 

709 .. math:: (n - 1) y(n + 2) - (n^2 + 3 n - 2) y(n + 1) + 

710 2 n (n + 1) y(n) = 0 

711 

712 >>> from sympy import Function, rsolve 

713 >>> from sympy.abc import n 

714 >>> y = Function('y') 

715 

716 >>> f = (n - 1)*y(n + 2) - (n**2 + 3*n - 2)*y(n + 1) + 2*n*(n + 1)*y(n) 

717 

718 >>> rsolve(f, y(n)) 

719 2**n*C0 + C1*factorial(n) 

720 

721 >>> rsolve(f, y(n), {y(0):0, y(1):3}) 

722 3*2**n - 3*factorial(n) 

723 

724 See Also 

725 ======== 

726 

727 rsolve_poly, rsolve_ratio, rsolve_hyper 

728 

729 """ 

730 if isinstance(f, Equality): 

731 f = f.lhs - f.rhs 

732 

733 n = y.args[0] 

734 k = Wild('k', exclude=(n,)) 

735 

736 # Preprocess user input to allow things like 

737 # y(n) + a*(y(n + 1) + y(n - 1))/2 

738 f = f.expand().collect(y.func(Wild('m', integer=True))) 

739 

740 h_part = defaultdict(list) 

741 i_part = [] 

742 for g in Add.make_args(f): 

743 coeff, dep = g.as_coeff_mul(y.func) 

744 if not dep: 

745 i_part.append(coeff) 

746 continue 

747 for h in dep: 

748 if h.is_Function and h.func == y.func: 

749 result = h.args[0].match(n + k) 

750 if result is not None: 

751 h_part[int(result[k])].append(coeff) 

752 continue 

753 raise ValueError( 

754 "'%s(%s + k)' expected, got '%s'" % (y.func, n, h)) 

755 for k in h_part: 

756 h_part[k] = Add(*h_part[k]) 

757 h_part.default_factory = lambda: 0 

758 i_part = Add(*i_part) 

759 

760 for k, coeff in h_part.items(): 

761 h_part[k] = simplify(coeff) 

762 

763 common = S.One 

764 

765 if not i_part.is_zero and not i_part.is_hypergeometric(n) and \ 

766 not (i_part.is_Add and all((x.is_hypergeometric(n) for x in i_part.expand().args))): 

767 raise ValueError("The independent term should be a sum of hypergeometric functions, got '%s'" % i_part) 

768 

769 for coeff in h_part.values(): 

770 if coeff.is_rational_function(n): 

771 if not coeff.is_polynomial(n): 

772 common = lcm(common, coeff.as_numer_denom()[1], n) 

773 else: 

774 raise ValueError( 

775 "Polynomial or rational function expected, got '%s'" % coeff) 

776 

777 i_numer, i_denom = i_part.as_numer_denom() 

778 

779 if i_denom.is_polynomial(n): 

780 common = lcm(common, i_denom, n) 

781 

782 if common is not S.One: 

783 for k, coeff in h_part.items(): 

784 numer, denom = coeff.as_numer_denom() 

785 h_part[k] = numer*quo(common, denom, n) 

786 

787 i_part = i_numer*quo(common, i_denom, n) 

788 

789 K_min = min(h_part.keys()) 

790 

791 if K_min < 0: 

792 K = abs(K_min) 

793 

794 H_part = defaultdict(lambda: S.Zero) 

795 i_part = i_part.subs(n, n + K).expand() 

796 common = common.subs(n, n + K).expand() 

797 

798 for k, coeff in h_part.items(): 

799 H_part[k + K] = coeff.subs(n, n + K).expand() 

800 else: 

801 H_part = h_part 

802 

803 K_max = max(H_part.keys()) 

804 coeffs = [H_part[i] for i in range(K_max + 1)]