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================================================== 

3Legendre Series (:mod:`numpy.polynomial.legendre`) 

4================================================== 

5 

6This module provides a number of objects (mostly functions) useful for 

7dealing with Legendre series, including a `Legendre` class that 

8encapsulates the usual arithmetic operations. (General information 

9on how this module represents and works with such polynomials is in the 

10docstring for its "parent" sub-package, `numpy.polynomial`). 

11 

12Classes 

13------- 

14.. autosummary:: 

15 :toctree: generated/ 

16 

17 Legendre 

18 

19Constants 

20--------- 

21 

22.. autosummary:: 

23 :toctree: generated/ 

24 

25 legdomain 

26 legzero 

27 legone 

28 legx 

29 

30Arithmetic 

31---------- 

32 

33.. autosummary:: 

34 :toctree: generated/ 

35 

36 legadd 

37 legsub 

38 legmulx 

39 legmul 

40 legdiv 

41 legpow 

42 legval 

43 legval2d 

44 legval3d 

45 leggrid2d 

46 leggrid3d 

47 

48Calculus 

49-------- 

50 

51.. autosummary:: 

52 :toctree: generated/ 

53 

54 legder 

55 legint 

56 

57Misc Functions 

58-------------- 

59 

60.. autosummary:: 

61 :toctree: generated/ 

62 

63 legfromroots 

64 legroots 

65 legvander 

66 legvander2d 

67 legvander3d 

68 leggauss 

69 legweight 

70 legcompanion 

71 legfit 

72 legtrim 

73 legline 

74 leg2poly 

75 poly2leg 

76 

77See also 

78-------- 

79numpy.polynomial 

80 

81""" 

82import numpy as np 

83import numpy.linalg as la 

84from numpy.core.multiarray import normalize_axis_index 

85 

86from . import polyutils as pu 

87from ._polybase import ABCPolyBase 

88 

89__all__ = [ 

90 'legzero', 'legone', 'legx', 'legdomain', 'legline', 'legadd', 

91 'legsub', 'legmulx', 'legmul', 'legdiv', 'legpow', 'legval', 'legder', 

92 'legint', 'leg2poly', 'poly2leg', 'legfromroots', 'legvander', 

93 'legfit', 'legtrim', 'legroots', 'Legendre', 'legval2d', 'legval3d', 

94 'leggrid2d', 'leggrid3d', 'legvander2d', 'legvander3d', 'legcompanion', 

95 'leggauss', 'legweight'] 

96 

97legtrim = pu.trimcoef 

98 

99 

100def poly2leg(pol): 

101 """ 

102 Convert a polynomial to a Legendre series. 

103 

104 Convert an array representing the coefficients of a polynomial (relative 

105 to the "standard" basis) ordered from lowest degree to highest, to an 

106 array of the coefficients of the equivalent Legendre series, ordered 

107 from lowest to highest degree. 

108 

109 Parameters 

110 ---------- 

111 pol : array_like 

112 1-D array containing the polynomial coefficients 

113 

114 Returns 

115 ------- 

116 c : ndarray 

117 1-D array containing the coefficients of the equivalent Legendre 

118 series. 

119 

120 See Also 

121 -------- 

122 leg2poly 

123 

124 Notes 

125 ----- 

126 The easy way to do conversions between polynomial basis sets 

127 is to use the convert method of a class instance. 

128 

129 Examples 

130 -------- 

131 >>> from numpy import polynomial as P 

132 >>> p = P.Polynomial(np.arange(4)) 

133 >>> p 

134 Polynomial([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1]) 

135 >>> c = P.Legendre(P.legendre.poly2leg(p.coef)) 

136 >>> c 

137 Legendre([ 1. , 3.25, 1. , 0.75], domain=[-1, 1], window=[-1, 1]) # may vary 

138 

139 """ 

140 [pol] = pu.as_series([pol]) 

141 deg = len(pol) - 1 

142 res = 0 

143 for i in range(deg, -1, -1): 

144 res = legadd(legmulx(res), pol[i]) 

145 return res 

146 

147 

148def leg2poly(c): 

149 """ 

150 Convert a Legendre series to a polynomial. 

151 

152 Convert an array representing the coefficients of a Legendre series, 

153 ordered from lowest degree to highest, to an array of the coefficients 

154 of the equivalent polynomial (relative to the "standard" basis) ordered 

155 from lowest to highest degree. 

156 

157 Parameters 

158 ---------- 

159 c : array_like 

160 1-D array containing the Legendre series coefficients, ordered 

161 from lowest order term to highest. 

162 

163 Returns 

164 ------- 

165 pol : ndarray 

166 1-D array containing the coefficients of the equivalent polynomial 

167 (relative to the "standard" basis) ordered from lowest order term 

168 to highest. 

169 

170 See Also 

171 -------- 

172 poly2leg 

173 

174 Notes 

175 ----- 

176 The easy way to do conversions between polynomial basis sets 

177 is to use the convert method of a class instance. 

178 

179 Examples 

180 -------- 

181 >>> from numpy import polynomial as P 

182 >>> c = P.Legendre(range(4)) 

183 >>> c 

184 Legendre([0., 1., 2., 3.], domain=[-1, 1], window=[-1, 1]) 

185 >>> p = c.convert(kind=P.Polynomial) 

186 >>> p 

187 Polynomial([-1. , -3.5, 3. , 7.5], domain=[-1., 1.], window=[-1., 1.]) 

188 >>> P.leg2poly(range(4)) 

189 array([-1. , -3.5, 3. , 7.5]) 

190 

191 

192 """ 

193 from .polynomial import polyadd, polysub, polymulx 

194 

195 [c] = pu.as_series([c]) 

196 n = len(c) 

197 if n < 3: 

198 return c 

199 else: 

200 c0 = c[-2] 

201 c1 = c[-1] 

202 # i is the current degree of c1 

203 for i in range(n - 1, 1, -1): 

204 tmp = c0 

205 c0 = polysub(c[i - 2], (c1*(i - 1))/i) 

206 c1 = polyadd(tmp, (polymulx(c1)*(2*i - 1))/i) 

207 return polyadd(c0, polymulx(c1)) 

208 

209# 

210# These are constant arrays are of integer type so as to be compatible 

211# with the widest range of other types, such as Decimal. 

212# 

213 

214# Legendre 

215legdomain = np.array([-1, 1]) 

216 

217# Legendre coefficients representing zero. 

218legzero = np.array([0]) 

219 

220# Legendre coefficients representing one. 

221legone = np.array([1]) 

222 

223# Legendre coefficients representing the identity x. 

224legx = np.array([0, 1]) 

225 

226 

227def legline(off, scl): 

228 """ 

229 Legendre series whose graph is a straight line. 

230 

231 

232 

233 Parameters 

234 ---------- 

235 off, scl : scalars 

236 The specified line is given by ``off + scl*x``. 

237 

238 Returns 

239 ------- 

240 y : ndarray 

241 This module's representation of the Legendre series for 

242 ``off + scl*x``. 

243 

244 See Also 

245 -------- 

246 polyline, chebline 

247 

248 Examples 

249 -------- 

250 >>> import numpy.polynomial.legendre as L 

251 >>> L.legline(3,2) 

252 array([3, 2]) 

253 >>> L.legval(-3, L.legline(3,2)) # should be -3 

254 -3.0 

255 

256 """ 

257 if scl != 0: 

258 return np.array([off, scl]) 

259 else: 

260 return np.array([off]) 

261 

262 

263def legfromroots(roots): 

264 """ 

265 Generate a Legendre series with given roots. 

266 

267 The function returns the coefficients of the polynomial 

268 

269 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n), 

270 

271 in Legendre form, where the `r_n` are the roots specified in `roots`. 

272 If a zero has multiplicity n, then it must appear in `roots` n times. 

273 For instance, if 2 is a root of multiplicity three and 3 is a root of 

274 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The 

275 roots can appear in any order. 

276 

277 If the returned coefficients are `c`, then 

278 

279 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x) 

280 

281 The coefficient of the last term is not generally 1 for monic 

282 polynomials in Legendre form. 

283 

284 Parameters 

285 ---------- 

286 roots : array_like 

287 Sequence containing the roots. 

288 

289 Returns 

290 ------- 

291 out : ndarray 

292 1-D array of coefficients. If all roots are real then `out` is a 

293 real array, if some of the roots are complex, then `out` is complex 

294 even if all the coefficients in the result are real (see Examples 

295 below). 

296 

297 See Also 

298 -------- 

299 polyfromroots, chebfromroots, lagfromroots, hermfromroots, hermefromroots 

300 

301 Examples 

302 -------- 

303 >>> import numpy.polynomial.legendre as L 

304 >>> L.legfromroots((-1,0,1)) # x^3 - x relative to the standard basis 

305 array([ 0. , -0.4, 0. , 0.4]) 

306 >>> j = complex(0,1) 

307 >>> L.legfromroots((-j,j)) # x^2 + 1 relative to the standard basis 

308 array([ 1.33333333+0.j, 0.00000000+0.j, 0.66666667+0.j]) # may vary 

309 

310 """ 

311 return pu._fromroots(legline, legmul, roots) 

312 

313 

314def legadd(c1, c2): 

315 """ 

316 Add one Legendre series to another. 

317 

318 Returns the sum of two Legendre series `c1` + `c2`. The arguments 

319 are sequences of coefficients ordered from lowest order term to 

320 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. 

321 

322 Parameters 

323 ---------- 

324 c1, c2 : array_like 

325 1-D arrays of Legendre series coefficients ordered from low to 

326 high. 

327 

328 Returns 

329 ------- 

330 out : ndarray 

331 Array representing the Legendre series of their sum. 

332 

333 See Also 

334 -------- 

335 legsub, legmulx, legmul, legdiv, legpow 

336 

337 Notes 

338 ----- 

339 Unlike multiplication, division, etc., the sum of two Legendre series 

340 is a Legendre series (without having to "reproject" the result onto 

341 the basis set) so addition, just like that of "standard" polynomials, 

342 is simply "component-wise." 

343 

344 Examples 

345 -------- 

346 >>> from numpy.polynomial import legendre as L 

347 >>> c1 = (1,2,3) 

348 >>> c2 = (3,2,1) 

349 >>> L.legadd(c1,c2) 

350 array([4., 4., 4.]) 

351 

352 """ 

353 return pu._add(c1, c2) 

354 

355 

356def legsub(c1, c2): 

357 """ 

358 Subtract one Legendre series from another. 

359 

360 Returns the difference of two Legendre series `c1` - `c2`. The 

361 sequences of coefficients are from lowest order term to highest, i.e., 

362 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. 

363 

364 Parameters 

365 ---------- 

366 c1, c2 : array_like 

367 1-D arrays of Legendre series coefficients ordered from low to 

368 high. 

369 

370 Returns 

371 ------- 

372 out : ndarray 

373 Of Legendre series coefficients representing their difference. 

374 

375 See Also 

376 -------- 

377 legadd, legmulx, legmul, legdiv, legpow 

378 

379 Notes 

380 ----- 

381 Unlike multiplication, division, etc., the difference of two Legendre 

382 series is a Legendre series (without having to "reproject" the result 

383 onto the basis set) so subtraction, just like that of "standard" 

384 polynomials, is simply "component-wise." 

385 

386 Examples 

387 -------- 

388 >>> from numpy.polynomial import legendre as L 

389 >>> c1 = (1,2,3) 

390 >>> c2 = (3,2,1) 

391 >>> L.legsub(c1,c2) 

392 array([-2., 0., 2.]) 

393 >>> L.legsub(c2,c1) # -C.legsub(c1,c2) 

394 array([ 2., 0., -2.]) 

395 

396 """ 

397 return pu._sub(c1, c2) 

398 

399 

400def legmulx(c): 

401 """Multiply a Legendre series by x. 

402 

403 Multiply the Legendre series `c` by x, where x is the independent 

404 variable. 

405 

406 

407 Parameters 

408 ---------- 

409 c : array_like 

410 1-D array of Legendre series coefficients ordered from low to 

411 high. 

412 

413 Returns 

414 ------- 

415 out : ndarray 

416 Array representing the result of the multiplication. 

417 

418 See Also 

419 -------- 

420 legadd, legmul, legmul, legdiv, legpow 

421 

422 Notes 

423 ----- 

424 The multiplication uses the recursion relationship for Legendre 

425 polynomials in the form 

426 

427 .. math:: 

428 

429 xP_i(x) = ((i + 1)*P_{i + 1}(x) + i*P_{i - 1}(x))/(2i + 1) 

430 

431 Examples 

432 -------- 

433 >>> from numpy.polynomial import legendre as L 

434 >>> L.legmulx([1,2,3]) 

435 array([ 0.66666667, 2.2, 1.33333333, 1.8]) # may vary 

436 

437 """ 

438 # c is a trimmed copy 

439 [c] = pu.as_series([c]) 

440 # The zero series needs special treatment 

441 if len(c) == 1 and c[0] == 0: 

442 return c 

443 

444 prd = np.empty(len(c) + 1, dtype=c.dtype) 

445 prd[0] = c[0]*0 

446 prd[1] = c[0] 

447 for i in range(1, len(c)): 

448 j = i + 1 

449 k = i - 1 

450 s = i + j 

451 prd[j] = (c[i]*j)/s 

452 prd[k] += (c[i]*i)/s 

453 return prd 

454 

455 

456def legmul(c1, c2): 

457 """ 

458 Multiply one Legendre series by another. 

459 

460 Returns the product of two Legendre series `c1` * `c2`. The arguments 

461 are sequences of coefficients, from lowest order "term" to highest, 

462 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``. 

463 

464 Parameters 

465 ---------- 

466 c1, c2 : array_like 

467 1-D arrays of Legendre series coefficients ordered from low to 

468 high. 

469 

470 Returns 

471 ------- 

472 out : ndarray 

473 Of Legendre series coefficients representing their product. 

474 

475 See Also 

476 -------- 

477 legadd, legsub, legmulx, legdiv, legpow 

478 

479 Notes 

480 ----- 

481 In general, the (polynomial) product of two C-series results in terms 

482 that are not in the Legendre polynomial basis set. Thus, to express 

483 the product as a Legendre series, it is necessary to "reproject" the 

484 product onto said basis set, which may produce "unintuitive" (but 

485 correct) results; see Examples section below. 

486 

487 Examples 

488 -------- 

489 >>> from numpy.polynomial import legendre as L 

490 >>> c1 = (1,2,3) 

491 >>> c2 = (3,2) 

492 >>> L.legmul(c1,c2) # multiplication requires "reprojection" 

493 array([ 4.33333333, 10.4 , 11.66666667, 3.6 ]) # may vary 

494 

495 """ 

496 # s1, s2 are trimmed copies 

497 [c1, c2] = pu.as_series([c1, c2]) 

498 

499 if len(c1) > len(c2): 

500 c = c2 

501 xs = c1 

502 else: 

503 c = c1 

504 xs = c2 

505 

506 if len(c) == 1: 

507 c0 = c[0]*xs 

508 c1 = 0 

509 elif len(c) == 2: 

510 c0 = c[0]*xs 

511 c1 = c[1]*xs 

512 else: 

513 nd = len(c) 

514 c0 = c[-2]*xs 

515 c1 = c[-1]*xs 

516 for i in range(3, len(c) + 1): 

517 tmp = c0 

518 nd = nd - 1 

519 c0 = legsub(c[-i]*xs, (c1*(nd - 1))/nd) 

520 c1 = legadd(tmp, (legmulx(c1)*(2*nd - 1))/nd) 

521 return legadd(c0, legmulx(c1)) 

522 

523 

524def legdiv(c1, c2): 

525 """ 

526 Divide one Legendre series by another. 

527 

528 Returns the quotient-with-remainder of two Legendre series 

529 `c1` / `c2`. The arguments are sequences of coefficients from lowest 

530 order "term" to highest, e.g., [1,2,3] represents the series 

531 ``P_0 + 2*P_1 + 3*P_2``. 

532 

533 Parameters 

534 ---------- 

535 c1, c2 : array_like 

536 1-D arrays of Legendre series coefficients ordered from low to 

537 high. 

538 

539 Returns 

540 ------- 

541 quo, rem : ndarrays 

542 Of Legendre series coefficients representing the quotient and 

543 remainder. 

544 

545 See Also 

546 -------- 

547 legadd, legsub, legmulx, legmul, legpow 

548 

549 Notes 

550 ----- 

551 In general, the (polynomial) division of one Legendre series by another 

552 results in quotient and remainder terms that are not in the Legendre 

553 polynomial basis set. Thus, to express these results as a Legendre 

554 series, it is necessary to "reproject" the results onto the Legendre 

555 basis set, which may produce "unintuitive" (but correct) results; see 

556 Examples section below. 

557 

558 Examples 

559 -------- 

560 >>> from numpy.polynomial import legendre as L 

561 >>> c1 = (1,2,3) 

562 >>> c2 = (3,2,1) 

563 >>> L.legdiv(c1,c2) # quotient "intuitive," remainder not 

564 (array([3.]), array([-8., -4.])) 

565 >>> c2 = (0,1,2,3) 

566 >>> L.legdiv(c2,c1) # neither "intuitive" 

567 (array([-0.07407407, 1.66666667]), array([-1.03703704, -2.51851852])) # may vary 

568 

569 """ 

570 return pu._div(legmul, c1, c2) 

571 

572 

573def legpow(c, pow, maxpower=16): 

574 """Raise a Legendre series to a power. 

575 

576 Returns the Legendre series `c` raised to the power `pow`. The 

577 argument `c` is a sequence of coefficients ordered from low to high. 

578 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.`` 

579 

580 Parameters 

581 ---------- 

582 c : array_like 

583 1-D array of Legendre series coefficients ordered from low to 

584 high. 

585 pow : integer 

586 Power to which the series will be raised 

587 maxpower : integer, optional 

588 Maximum power allowed. This is mainly to limit growth of the series 

589 to unmanageable size. Default is 16 

590 

591 Returns 

592 ------- 

593 coef : ndarray 

594 Legendre series of power. 

595 

596 See Also 

597 -------- 

598 legadd, legsub, legmulx, legmul, legdiv 

599 

600 Examples 

601 -------- 

602 

603 """ 

604 return pu._pow(legmul, c, pow, maxpower) 

605 

606 

607def legder(c, m=1, scl=1, axis=0): 

608 """ 

609 Differentiate a Legendre series. 

610 

611 Returns the Legendre series coefficients `c` differentiated `m` times 

612 along `axis`. At each iteration the result is multiplied by `scl` (the 

613 scaling factor is for use in a linear change of variable). The argument 

614 `c` is an array of coefficients from low to high degree along each 

615 axis, e.g., [1,2,3] represents the series ``1*L_0 + 2*L_1 + 3*L_2`` 

616 while [[1,2],[1,2]] represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 

617 2*L_0(x)*L_1(y) + 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is 

618 ``y``. 

619 

620 Parameters 

621 ---------- 

622 c : array_like 

623 Array of Legendre series coefficients. If c is multidimensional the 

624 different axis correspond to different variables with the degree in 

625 each axis given by the corresponding index. 

626 m : int, optional 

627 Number of derivatives taken, must be non-negative. (Default: 1) 

628 scl : scalar, optional 

629 Each differentiation is multiplied by `scl`. The end result is 

630 multiplication by ``scl**m``. This is for use in a linear change of 

631 variable. (Default: 1) 

632 axis : int, optional 

633 Axis over which the derivative is taken. (Default: 0). 

634 

635 .. versionadded:: 1.7.0 

636 

637 Returns 

638 ------- 

639 der : ndarray 

640 Legendre series of the derivative. 

641 

642 See Also 

643 -------- 

644 legint 

645 

646 Notes 

647 ----- 

648 In general, the result of differentiating a Legendre series does not 

649 resemble the same operation on a power series. Thus the result of this 

650 function may be "unintuitive," albeit correct; see Examples section 

651 below. 

652 

653 Examples 

654 -------- 

655 >>> from numpy.polynomial import legendre as L 

656 >>> c = (1,2,3,4) 

657 >>> L.legder(c) 

658 array([ 6., 9., 20.]) 

659 >>> L.legder(c, 3) 

660 array([60.]) 

661 >>> L.legder(c, scl=-1) 

662 array([ -6., -9., -20.]) 

663 >>> L.legder(c, 2,-1) 

664 array([ 9., 60.]) 

665 

666 """ 

667 c = np.array(c, ndmin=1, copy=True) 

668 if c.dtype.char in '?bBhHiIlLqQpP': 

669 c = c.astype(np.double) 

670 cnt = pu._deprecate_as_int(m, "the order of derivation") 

671 iaxis = pu._deprecate_as_int(axis, "the axis") 

672 if cnt < 0: 

673 raise ValueError("The order of derivation must be non-negative") 

674 iaxis = normalize_axis_index(iaxis, c.ndim) 

675 

676 if cnt == 0: 

677 return c 

678 

679 c = np.moveaxis(c, iaxis, 0) 

680 n = len(c) 

681 if cnt >= n: 

682 c = c[:1]*0 

683 else: 

684 for i in range(cnt): 

685 n = n - 1 

686 c *= scl 

687 der = np.empty((n,) + c.shape[1:], dtype=c.dtype) 

688 for j in range(n, 2, -1): 

689 der[j - 1] = (2*j - 1)*c[j] 

690 c[j - 2] += c[j] 

691 if n > 1: 

692 der[1] = 3*c[2] 

693 der[0] = c[1] 

694 c = der 

695 c = np.moveaxis(c, 0, iaxis) 

696 return c 

697 

698 

699def legint(c, m=1, k=[], lbnd=0, scl=1, axis=0): 

700 """ 

701 Integrate a Legendre series. 

702 

703 Returns the Legendre series coefficients `c` integrated `m` times from 

704 `lbnd` along `axis`. At each iteration the resulting series is 

705 **multiplied** by `scl` and an integration constant, `k`, is added. 

706 The scaling factor is for use in a linear change of variable. ("Buyer 

707 beware": note that, depending on what one is doing, one may want `scl` 

708 to be the reciprocal of what one might expect; for more information, 

709 see the Notes section below.) The argument `c` is an array of 

710 coefficients from low to high degree along each axis, e.g., [1,2,3] 

711 represents the series ``L_0 + 2*L_1 + 3*L_2`` while [[1,2],[1,2]] 

712 represents ``1*L_0(x)*L_0(y) + 1*L_1(x)*L_0(y) + 2*L_0(x)*L_1(y) + 

713 2*L_1(x)*L_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``. 

714 

715 Parameters 

716 ---------- 

717 c : array_like 

718 Array of Legendre series coefficients. If c is multidimensional the 

719 different axis correspond to different variables with the degree in 

720 each axis given by the corresponding index. 

721 m : int, optional 

722 Order of integration, must be positive. (Default: 1) 

723 k : {[], list, scalar}, optional 

724 Integration constant(s). The value of the first integral at 

725 ``lbnd`` is the first value in the list, the value of the second 

726 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the 

727 default), all constants are set to zero. If ``m == 1``, a single 

728 scalar can be given instead of a list. 

729 lbnd : scalar, optional 

730 The lower bound of the integral. (Default: 0) 

731 scl : scalar, optional 

732 Following each integration the result is *multiplied* by `scl` 

733 before the integration constant is added. (Default: 1) 

734 axis : int, optional 

735 Axis over which the integral is taken. (Default: 0). 

736 

737 .. versionadded:: 1.7.0 

738 

739 Returns 

740 ------- 

741 S : ndarray 

742 Legendre series coefficient array of the integral. 

743 

744 Raises 

745 ------ 

746 ValueError 

747 If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or 

748 ``np.ndim(scl) != 0``. 

749 

750 See Also 

751 -------- 

752 legder 

753 

754 Notes 

755 ----- 

756 Note that the result of each integration is *multiplied* by `scl`. 

757 Why is this important to note? Say one is making a linear change of 

758 variable :math:`u = ax + b` in an integral relative to `x`. Then 

759 :math:`dx = du/a`, so one will need to set `scl` equal to 

760 :math:`1/a` - perhaps not what one would have first thought. 

761 

762 Also note that, in general, the result of integrating a C-series needs 

763 to be "reprojected" onto the C-series basis set. Thus, typically, 

764 the result of this function is "unintuitive," albeit correct; see 

765 Examples section below. 

766 

767 Examples 

768 -------- 

769 >>> from numpy.polynomial import legendre as L 

770 >>> c = (1,2,3) 

771 >>> L.legint(c) 

772 array([ 0.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary 

773 >>> L.legint(c, 3) 

774 array([ 1.66666667e-02, -1.78571429e-02, 4.76190476e-02, # may vary 

775 -1.73472348e-18, 1.90476190e-02, 9.52380952e-03]) 

776 >>> L.legint(c, k=3) 

777 array([ 3.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary 

778 >>> L.legint(c, lbnd=-2) 

779 array([ 7.33333333, 0.4 , 0.66666667, 0.6 ]) # may vary 

780 >>> L.legint(c, scl=2) 

781 array([ 0.66666667, 0.8 , 1.33333333, 1.2 ]) # may vary 

782 

783 """ 

784 c = np.array(c, ndmin=1, copy=True) 

785 if c.dtype.char in '?bBhHiIlLqQpP': 

786 c = c.astype(np.double) 

787 if not np.iterable(k): 

788 k = [k] 

789 cnt = pu._deprecate_as_int(m, "the order of integration") 

790 iaxis = pu._deprecate_as_int(axis, "the axis") 

791 if cnt < 0: 

792 raise ValueError("The order of integration must be non-negative") 

793 if len(k) > cnt: 

794 raise ValueError("Too many integration constants") 

795 if np.ndim(lbnd) != 0: 

796 raise ValueError("lbnd must be a scalar.") 

797 if np.ndim(scl) != 0: 

798 raise ValueError("scl must be a scalar.") 

799 iaxis = normalize_axis_index(iaxis, c.ndim) 

800 

801 if cnt == 0: 

802 return c 

803 

804 c = np.moveaxis(c, iaxis, 0) 

805 k = list(k) + [0]*(cnt - len(k)) 

806 for i in range(cnt): 

807 n = len(c) 

808 c *= scl 

809 if n == 1 and np.all(c[0] == 0): 

810 c[0] += k[i] 

811 else: 

812 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype) 

813 tmp[0] = c[0]*0 

814 tmp[1] = c[0] 

815 if n > 1: 

816 tmp[2] = c[1]/3 

817 for j in range(2, n): 

818 t = c[j]/(2*j + 1) 

819 tmp[j + 1] = t 

820 tmp[j - 1] -= t 

821 tmp[0] += k[i] - legval(lbnd, tmp) 

822 c = tmp 

823 c = np.moveaxis(c, 0, iaxis) 

824 return c 

825 

826 

827def legval(x, c, tensor=True): 

828 """ 

829 Evaluate a Legendre series at points x. 

830 

831 If `c` is of length `n + 1`, this function returns the value: 

832 

833 .. math:: p(x) = c_0 * L_0(x) + c_1 * L_1(x) + ... + c_n * L_n(x) 

834 

835 The parameter `x` is converted to an array only if it is a tuple or a 

836 list, otherwise it is treated as a scalar. In either case, either `x` 

837 or its elements must support multiplication and addition both with 

838 themselves and with the elements of `c`. 

839 

840 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If 

841 `c` is multidimensional, then the shape of the result depends on the 

842 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] + 

843 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that 

844 scalars have shape (,). 

845 

846 Trailing zeros in the coefficients will be used in the evaluation, so 

847 they should be avoided if efficiency is a concern. 

848 

849 Parameters 

850 ---------- 

851 x : array_like, compatible object 

852 If `x` is a list or tuple, it is converted to an ndarray, otherwise 

853 it is left unchanged and treated as a scalar. In either case, `x` 

854 or its elements must support addition and multiplication with 

855 with themselves and with the elements of `c`. 

856 c : array_like 

857 Array of coefficients ordered so that the coefficients for terms of 

858 degree n are contained in c[n]. If `c` is multidimensional the 

859 remaining indices enumerate multiple polynomials. In the two 

860 dimensional case the coefficients may be thought of as stored in 

861 the columns of `c`. 

862 tensor : boolean, optional 

863 If True, the shape of the coefficient array is extended with ones 

864 on the right, one for each dimension of `x`. Scalars have dimension 0 

865 for this action. The result is that every column of coefficients in 

866 `c` is evaluated for every element of `x`. If False, `x` is broadcast 

867 over the columns of `c` for the evaluation. This keyword is useful 

868 when `c` is multidimensional. The default value is True. 

869 

870 .. versionadded:: 1.7.0 

871 

872 Returns 

873 ------- 

874 values : ndarray, algebra_like 

875 The shape of the return value is described above. 

876 

877 See Also 

878 -------- 

879 legval2d, leggrid2d, legval3d, leggrid3d 

880 

881 Notes 

882 ----- 

883 The evaluation uses Clenshaw recursion, aka synthetic division. 

884 

885 Examples 

886 -------- 

887 

888 """ 

889 c = np.array(c, ndmin=1, copy=False) 

890 if c.dtype.char in '?bBhHiIlLqQpP': 

891 c = c.astype(np.double) 

892 if isinstance(x, (tuple, list)): 

893 x = np.asarray(x) 

894 if isinstance(x, np.ndarray) and tensor: 

895 c = c.reshape(c.shape + (1,)*x.ndim) 

896 

897 if len(c) == 1: 

898 c0 = c[0] 

899 c1 = 0 

900 elif len(c) == 2: 

901 c0 = c[0] 

902 c1 = c[1] 

903 else: 

904 nd = len(c) 

905 c0 = c[-2] 

906 c1 = c[-1] 

907 for i in range(3, len(c) + 1): 

908 tmp = c0 

909 nd = nd - 1 

910 c0 = c[-i] - (c1*(nd - 1))/nd 

911 c1 = tmp + (c1*x*(2*nd - 1))/nd 

912 return c0 + c1*x 

913 

914 

915def legval2d(x, y, c): 

916 """ 

917 Evaluate a 2-D Legendre series at points (x, y). 

918 

919 This function returns the values: 

920 

921 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * L_i(x) * L_j(y) 

922 

923 The parameters `x` and `y` are converted to arrays only if they are 

924 tuples or a lists, otherwise they are treated as a scalars and they 

925 must have the same shape after conversion. In either case, either `x` 

926 and `y` or their elements must support multiplication and addition both 

927 with themselves and with the elements of `c`. 

928 

929 If `c` is a 1-D array a one is implicitly appended to its shape to make 

930 it 2-D. The shape of the result will be c.shape[2:] + x.shape. 

931 

932 Parameters 

933 ---------- 

934 x, y : array_like, compatible objects 

935 The two dimensional series is evaluated at the points `(x, y)`, 

936 where `x` and `y` must have the same shape. If `x` or `y` is a list 

937 or tuple, it is first converted to an ndarray, otherwise it is left 

938 unchanged and if it isn't an ndarray it is treated as a scalar. 

939 c : array_like 

940 Array of coefficients ordered so that the coefficient of the term 

941 of multi-degree i,j is contained in ``c[i,j]``. If `c` has 

942 dimension greater than two the remaining indices enumerate multiple 

943 sets of coefficients. 

944 

945 Returns 

946 ------- 

947 values : ndarray, compatible object 

948 The values of the two dimensional Legendre series at points formed 

949 from pairs of corresponding values from `x` and `y`. 

950 

951 See Also 

952 -------- 

953 legval, leggrid2d, legval3d, leggrid3d 

954 

955 Notes 

956 ----- 

957 

958 .. versionadded:: 1.7.0 

959 

960 """ 

961 return pu._valnd(legval, c, x, y) 

962 

963 

964def leggrid2d(x, y, c): 

965 """ 

966 Evaluate a 2-D Legendre series on the Cartesian product of x and y. 

967 

968 This function returns the values: 

969 

970 .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * L_i(a) * L_j(b) 

971 

972 where the points `(a, b)` consist of all pairs formed by taking 

973 `a` from `x` and `b` from `y`. The resulting points form a grid with 

974 `x` in the first dimension and `y` in the second. 

975 

976 The parameters `x` and `y` are converted to arrays only if they are 

977 tuples or a lists, otherwise they are treated as a scalars. In either 

978 case, either `x` and `y` or their elements must support multiplication 

979 and addition both with themselves and with the elements of `c`. 

980 

981 If `c` has fewer than two dimensions, ones are implicitly appended to 

982 its shape to make it 2-D. The shape of the result will be c.shape[2:] + 

983 x.shape + y.shape. 

984 

985 Parameters 

986 ---------- 

987 x, y : array_like, compatible objects 

988 The two dimensional series is evaluated at the points in the 

989 Cartesian product of `x` and `y`. If `x` or `y` is a list or 

990 tuple, it is first converted to an ndarray, otherwise it is left 

991 unchanged and, if it isn't an ndarray, it is treated as a scalar. 

992 c : array_like 

993 Array of coefficients ordered so that the coefficient of the term of 

994 multi-degree i,j is contained in `c[i,j]`. If `c` has dimension 

995 greater than two the remaining indices enumerate multiple sets of 

996 coefficients. 

997 

998 Returns 

999 ------- 

1000 values : ndarray, compatible object 

1001 The values of the two dimensional Chebyshev series at points in the 

1002 Cartesian product of `x` and `y`. 

1003 

1004 See Also 

1005 -------- 

1006 legval, legval2d, legval3d, leggrid3d 

1007 

1008 Notes 

1009 ----- 

1010 

1011 .. versionadded:: 1.7.0 

1012 

1013 """ 

1014 return pu._gridnd(legval, c, x, y) 

1015 

1016 

1017def legval3d(x, y, z, c): 

1018 """ 

1019 Evaluate a 3-D Legendre series at points (x, y, z). 

1020 

1021 This function returns the values: 

1022 

1023 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * L_i(x) * L_j(y) * L_k(z) 

1024 

1025 The parameters `x`, `y`, and `z` are converted to arrays only if 

1026 they are tuples or a lists, otherwise they are treated as a scalars and 

1027 they must have the same shape after conversion. In either case, either 

1028 `x`, `y`, and `z` or their elements must support multiplication and 

1029 addition both with themselves and with the elements of `c`. 

1030 

1031 If `c` has fewer than 3 dimensions, ones are implicitly appended to its 

1032 shape to make it 3-D. The shape of the result will be c.shape[3:] + 

1033 x.shape. 

1034 

1035 Parameters 

1036 ---------- 

1037 x, y, z : array_like, compatible object 

1038 The three dimensional series is evaluated at the points 

1039 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If 

1040 any of `x`, `y`, or `z` is a list or tuple, it is first converted 

1041 to an ndarray, otherwise it is left unchanged and if it isn't an 

1042 ndarray it is treated as a scalar. 

1043 c : array_like 

1044 Array of coefficients ordered so that the coefficient of the term of 

1045 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension 

1046 greater than 3 the remaining indices enumerate multiple sets of 

1047 coefficients. 

1048 

1049 Returns 

1050 ------- 

1051 values : ndarray, compatible object 

1052 The values of the multidimensional polynomial on points formed with 

1053 triples of corresponding values from `x`, `y`, and `z`. 

1054 

1055 See Also 

1056 -------- 

1057 legval, legval2d, leggrid2d, leggrid3d 

1058 

1059 Notes 

1060 ----- 

1061 

1062 .. versionadded:: 1.7.0 

1063 

1064 """ 

1065 return pu._valnd(legval, c, x, y, z) 

1066 

1067 

1068def leggrid3d(x, y, z, c): 

1069 """ 

1070 Evaluate a 3-D Legendre series on the Cartesian product of x, y, and z. 

1071 

1072 This function returns the values: 

1073 

1074 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * L_i(a) * L_j(b) * L_k(c) 

1075 

1076 where the points `(a, b, c)` consist of all triples formed by taking 

1077 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form 

1078 a grid with `x` in the first dimension, `y` in the second, and `z` in 

1079 the third. 

1080 

1081 The parameters `x`, `y`, and `z` are converted to arrays only if they 

1082 are tuples or a lists, otherwise they are treated as a scalars. In 

1083 either case, either `x`, `y`, and `z` or their elements must support 

1084 multiplication and addition both with themselves and with the elements 

1085 of `c`. 

1086 

1087 If `c` has fewer than three dimensions, ones are implicitly appended to 

1088 its shape to make it 3-D. The shape of the result will be c.shape[3:] + 

1089 x.shape + y.shape + z.shape. 

1090 

1091 Parameters 

1092 ---------- 

1093 x, y, z : array_like, compatible objects 

1094 The three dimensional series is evaluated at the points in the 

1095 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a 

1096 list or tuple, it is first converted to an ndarray, otherwise it is 

1097 left unchanged and, if it isn't an ndarray, it is treated as a 

1098 scalar. 

1099 c : array_like 

1100 Array of coefficients ordered so that the coefficients for terms of 

1101 degree i,j are contained in ``c[i,j]``. If `c` has dimension 

1102 greater than two the remaining indices enumerate multiple sets of 

1103 coefficients. 

1104 

1105 Returns 

1106 ------- 

1107 values : ndarray, compatible object 

1108 The values of the two dimensional polynomial at points in the Cartesian 

1109 product of `x` and `y`. 

1110 

1111 See Also 

1112 -------- 

1113 legval, legval2d, leggrid2d, legval3d 

1114 

1115 Notes 

1116 ----- 

1117 

1118 .. versionadded:: 1.7.0 

1119 

1120 """ 

1121 return pu._gridnd(legval, c, x, y, z) 

1122 

1123 

1124def legvander(x, deg): 

1125 """Pseudo-Vandermonde matrix of given degree. 

1126 

1127 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points 

1128 `x`. The pseudo-Vandermonde matrix is defined by 

1129 

1130 .. math:: V[..., i] = L_i(x) 

1131 

1132 where `0 <= i <= deg`. The leading indices of `V` index the elements of 

1133 `x` and the last index is the degree of the Legendre polynomial. 

1134 

1135 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the 

1136 array ``V = legvander(x, n)``, then ``np.dot(V, c)`` and 

1137 ``legval(x, c)`` are the same up to roundoff. This equivalence is 

1138 useful both for least squares fitting and for the evaluation of a large 

1139 number of Legendre series of the same degree and sample points. 

1140 

1141 Parameters 

1142 ---------- 

1143 x : array_like 

1144 Array of points. The dtype is converted to float64 or complex128 

1145 depending on whether any of the elements are complex. If `x` is 

1146 scalar it is converted to a 1-D array. 

1147 deg : int 

1148 Degree of the resulting matrix. 

1149 

1150 Returns 

1151 ------- 

1152 vander : ndarray 

1153 The pseudo-Vandermonde matrix. The shape of the returned matrix is 

1154 ``x.shape + (deg + 1,)``, where The last index is the degree of the 

1155 corresponding Legendre polynomial. The dtype will be the same as 

1156 the converted `x`. 

1157 

1158 """ 

1159 ideg = pu._deprecate_as_int(deg, "deg") 

1160 if ideg < 0: 

1161 raise ValueError("deg must be non-negative") 

1162 

1163 x = np.array(x, copy=False, ndmin=1) + 0.0 

1164 dims = (ideg + 1,) + x.shape 

1165 dtyp = x.dtype 

1166 v = np.empty(dims, dtype=dtyp) 

1167 # Use forward recursion to generate the entries. This is not as accurate 

1168 # as reverse recursion in this application but it is more efficient. 

1169 v[0] = x*0 + 1 

1170 if ideg > 0: 

1171 v[1] = x 

1172 for i in range(2, ideg + 1): 

1173 v[i] = (v[i-1]*x*(2*i - 1) - v[i-2]*(i - 1))/i 

1174 return np.moveaxis(v, 0, -1) 

1175 

1176 

1177def legvander2d(x, y, deg): 

1178 """Pseudo-Vandermonde matrix of given degrees. 

1179 

1180 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample 

1181 points `(x, y)`. The pseudo-Vandermonde matrix is defined by 

1182 

1183 .. math:: V[..., (deg[1] + 1)*i + j] = L_i(x) * L_j(y), 

1184 

1185 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of 

1186 `V` index the points `(x, y)` and the last index encodes the degrees of 

1187 the Legendre polynomials. 

1188 

1189 If ``V = legvander2d(x, y, [xdeg, ydeg])``, then the columns of `V` 

1190 correspond to the elements of a 2-D coefficient array `c` of shape 

1191 (xdeg + 1, ydeg + 1) in the order 

1192 

1193 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ... 

1194 

1195 and ``np.dot(V, c.flat)`` and ``legval2d(x, y, c)`` will be the same 

1196 up to roundoff. This equivalence is useful both for least squares 

1197 fitting and for the evaluation of a large number of 2-D Legendre 

1198 series of the same degrees and sample points. 

1199 

1200 Parameters 

1201 ---------- 

1202 x, y : array_like 

1203 Arrays of point coordinates, all of the same shape. The dtypes 

1204 will be converted to either float64 or complex128 depending on 

1205 whether any of the elements are complex. Scalars are converted to 

1206 1-D arrays. 

1207 deg : list of ints 

1208 List of maximum degrees of the form [x_deg, y_deg]. 

1209 

1210 Returns 

1211 ------- 

1212 vander2d : ndarray 

1213 The shape of the returned matrix is ``x.shape + (order,)``, where 

1214 :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same 

1215 as the converted `x` and `y`. 

1216 

1217 See Also 

1218 -------- 

1219 legvander, legvander3d, legval2d, legval3d 

1220 

1221 Notes 

1222 ----- 

1223 

1224 .. versionadded:: 1.7.0 

1225 

1226 """ 

1227 return pu._vander_nd_flat((legvander, legvander), (x, y), deg) 

1228 

1229 

1230def legvander3d(x, y, z, deg): 

1231 """Pseudo-Vandermonde matrix of given degrees. 

1232 

1233 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample 

1234 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`, 

1235 then The pseudo-Vandermonde matrix is defined by 

1236 

1237 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = L_i(x)*L_j(y)*L_k(z), 

1238 

1239 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading 

1240 indices of `V` index the points `(x, y, z)` and the last index encodes 

1241 the degrees of the Legendre polynomials. 

1242 

1243 If ``V = legvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns 

1244 of `V` correspond to the elements of a 3-D coefficient array `c` of 

1245 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order 

1246 

1247 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},... 

1248 

1249 and ``np.dot(V, c.flat)`` and ``legval3d(x, y, z, c)`` will be the 

1250 same up to roundoff. This equivalence is useful both for least squares 

1251 fitting and for the evaluation of a large number of 3-D Legendre 

1252 series of the same degrees and sample points. 

1253 

1254 Parameters 

1255 ---------- 

1256 x, y, z : array_like 

1257 Arrays of point coordinates, all of the same shape. The dtypes will 

1258 be converted to either float64 or complex128 depending on whether 

1259 any of the elements are complex. Scalars are converted to 1-D 

1260 arrays. 

1261 deg : list of ints 

1262 List of maximum degrees of the form [x_deg, y_deg, z_deg]. 

1263 

1264 Returns 

1265 ------- 

1266 vander3d : ndarray 

1267 The shape of the returned matrix is ``x.shape + (order,)``, where 

1268 :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will 

1269 be the same as the converted `x`, `y`, and `z`. 

1270 

1271 See Also 

1272 -------- 

1273 legvander, legvander3d, legval2d, legval3d 

1274 

1275 Notes 

1276 ----- 

1277 

1278 .. versionadded:: 1.7.0 

1279 

1280 """ 

1281 return pu._vander_nd_flat((legvander, legvander, legvander), (x, y, z), deg) 

1282 

1283 

1284def legfit(x, y, deg, rcond=None, full=False, w=None): 

1285 """ 

1286 Least squares fit of Legendre series to data. 

1287 

1288 Return the coefficients of a Legendre series of degree `deg` that is the 

1289 least squares fit to the data values `y` given at points `x`. If `y` is 

1290 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple 

1291 fits are done, one for each column of `y`, and the resulting 

1292 coefficients are stored in the corresponding columns of a 2-D return. 

1293 The fitted polynomial(s) are in the form 

1294 

1295 .. math:: p(x) = c_0 + c_1 * L_1(x) + ... + c_n * L_n(x), 

1296 

1297 where `n` is `deg`. 

1298 

1299 Parameters 

1300 ---------- 

1301 x : array_like, shape (M,) 

1302 x-coordinates of the M sample points ``(x[i], y[i])``. 

1303 y : array_like, shape (M,) or (M, K) 

1304 y-coordinates of the sample points. Several data sets of sample 

1305 points sharing the same x-coordinates can be fitted at once by 

1306 passing in a 2D-array that contains one dataset per column. 

1307 deg : int or 1-D array_like 

1308 Degree(s) of the fitting polynomials. If `deg` is a single integer 

1309 all terms up to and including the `deg`'th term are included in the 

1310 fit. For NumPy versions >= 1.11.0 a list of integers specifying the 

1311 degrees of the terms to include may be used instead. 

1312 rcond : float, optional 

1313 Relative condition number of the fit. Singular values smaller than 

1314 this relative to the largest singular value will be ignored. The 

1315 default value is len(x)*eps, where eps is the relative precision of 

1316 the float type, about 2e-16 in most cases. 

1317 full : bool, optional 

1318 Switch determining nature of return value. When it is False (the 

1319 default) just the coefficients are returned, when True diagnostic 

1320 information from the singular value decomposition is also returned. 

1321 w : array_like, shape (`M`,), optional 

1322 Weights. If not None, the contribution of each point 

1323 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the 

1324 weights are chosen so that the errors of the products ``w[i]*y[i]`` 

1325 all have the same variance. The default value is None. 

1326 

1327 .. versionadded:: 1.5.0 

1328 

1329 Returns 

1330 ------- 

1331 coef : ndarray, shape (M,) or (M, K) 

1332 Legendre coefficients ordered from low to high. If `y` was 

1333 2-D, the coefficients for the data in column k of `y` are in 

1334 column `k`. If `deg` is specified as a list, coefficients for 

1335 terms not included in the fit are set equal to zero in the 

1336 returned `coef`. 

1337 

1338 [residuals, rank, singular_values, rcond] : list 

1339 These values are only returned if `full` = True 

1340 

1341 resid -- sum of squared residuals of the least squares fit 

1342 rank -- the numerical rank of the scaled Vandermonde matrix 

1343 sv -- singular values of the scaled Vandermonde matrix 

1344 rcond -- value of `rcond`. 

1345 

1346 For more details, see `linalg.lstsq`. 

1347 

1348 Warns 

1349 ----- 

1350 RankWarning 

1351 The rank of the coefficient matrix in the least-squares fit is 

1352 deficient. The warning is only raised if `full` = False. The 

1353 warnings can be turned off by 

1354 

1355 >>> import warnings 

1356 >>> warnings.simplefilter('ignore', np.RankWarning) 

1357 

1358 See Also 

1359 -------- 

1360 chebfit, polyfit, lagfit, hermfit, hermefit 

1361 legval : Evaluates a Legendre series. 

1362 legvander : Vandermonde matrix of Legendre series. 

1363 legweight : Legendre weight function (= 1). 

1364 linalg.lstsq : Computes a least-squares fit from the matrix. 

1365 scipy.interpolate.UnivariateSpline : Computes spline fits. 

1366 

1367 Notes 

1368 ----- 

1369 The solution is the coefficients of the Legendre series `p` that 

1370 minimizes the sum of the weighted squared errors 

1371 

1372 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2, 

1373 

1374 where :math:`w_j` are the weights. This problem is solved by setting up 

1375 as the (typically) overdetermined matrix equation 

1376 

1377 .. math:: V(x) * c = w * y, 

1378 

1379 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the 

1380 coefficients to be solved for, `w` are the weights, and `y` are the 

1381 observed values. This equation is then solved using the singular value 

1382 decomposition of `V`. 

1383 

1384 If some of the singular values of `V` are so small that they are 

1385 neglected, then a `RankWarning` will be issued. This means that the 

1386 coefficient values may be poorly determined. Using a lower order fit 

1387 will usually get rid of the warning. The `rcond` parameter can also be 

1388 set to a value smaller than its default, but the resulting fit may be 

1389 spurious and have large contributions from roundoff error. 

1390 

1391 Fits using Legendre series are usually better conditioned than fits 

1392 using power series, but much can depend on the distribution of the 

1393 sample points and the smoothness of the data. If the quality of the fit 

1394 is inadequate splines may be a good alternative. 

1395 

1396 References 

1397 ---------- 

1398 .. [1] Wikipedia, "Curve fitting", 

1399 https://en.wikipedia.org/wiki/Curve_fitting 

1400 

1401 Examples 

1402 -------- 

1403 

1404 """ 

1405 return pu._fit(legvander, x, y, deg, rcond, full, w) 

1406 

1407 

1408def legcompanion(c): 

1409 """Return the scaled companion matrix of c. 

1410 

1411 The basis polynomials are scaled so that the companion matrix is 

1412 symmetric when `c` is an Legendre basis polynomial. This provides 

1413 better eigenvalue estimates than the unscaled case and for basis 

1414 polynomials the eigenvalues are guaranteed to be real if 

1415 `numpy.linalg.eigvalsh` is used to obtain them. 

1416 

1417 Parameters 

1418 ---------- 

1419 c : array_like 

1420 1-D array of Legendre series coefficients ordered from low to high 

1421 degree. 

1422 

1423 Returns 

1424 ------- 

1425 mat : ndarray 

1426 Scaled companion matrix of dimensions (deg, deg). 

1427 

1428 Notes 

1429 ----- 

1430 

1431 .. versionadded:: 1.7.0 

1432 

1433 """ 

1434 # c is a trimmed copy 

1435 [c] = pu.as_series([c]) 

1436 if len(c) < 2: 

1437 raise ValueError('Series must have maximum degree of at least 1.') 

1438 if len(c) == 2: 

1439 return np.array([[-c[0]/c[1]]]) 

1440 

1441 n = len(c) - 1 

1442 mat = np.zeros((n, n), dtype=c.dtype) 

1443 scl = 1./np.sqrt(2*np.arange(n) + 1) 

1444 top = mat.reshape(-1)[1::n+1] 

1445 bot = mat.reshape(-1)[n::n+1] 

1446 top[...] = np.arange(1, n)*scl[:n-1]*scl[1:n] 

1447 bot[...] = top 

1448 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*(n/(2*n - 1)) 

1449 return mat 

1450 

1451 

1452def legroots(c): 

1453 """ 

1454 Compute the roots of a Legendre series. 

1455 

1456 Return the roots (a.k.a. "zeros") of the polynomial 

1457 

1458 .. math:: p(x) = \\sum_i c[i] * L_i(x). 

1459 

1460 Parameters 

1461 ---------- 

1462 c : 1-D array_like 

1463 1-D array of coefficients. 

1464 

1465 Returns 

1466 ------- 

1467 out : ndarray 

1468 Array of the roots of the series. If all the roots are real, 

1469 then `out` is also real, otherwise it is complex. 

1470 

1471 See Also 

1472 -------- 

1473 polyroots, chebroots, lagroots, hermroots, hermeroots 

1474 

1475 Notes 

1476 ----- 

1477 The root estimates are obtained as the eigenvalues of the companion 

1478 matrix, Roots far from the origin of the complex plane may have large 

1479 errors due to the numerical instability of the series for such values. 

1480 Roots with multiplicity greater than 1 will also show larger errors as 

1481 the value of the series near such points is relatively insensitive to 

1482 errors in the roots. Isolated roots near the origin can be improved by 

1483 a few iterations of Newton's method. 

1484 

1485 The Legendre series basis polynomials aren't powers of ``x`` so the 

1486 results of this function may seem unintuitive. 

1487 

1488 Examples 

1489 -------- 

1490 >>> import numpy.polynomial.legendre as leg 

1491 >>> leg.legroots((1, 2, 3, 4)) # 4L_3 + 3L_2 + 2L_1 + 1L_0, all real roots 

1492 array([-0.85099543, -0.11407192, 0.51506735]) # may vary 

1493 

1494 """ 

1495 # c is a trimmed copy 

1496 [c] = pu.as_series([c]) 

1497 if len(c) < 2: 

1498 return np.array([], dtype=c.dtype) 

1499 if len(c) == 2: 

1500 return np.array([-c[0]/c[1]]) 

1501 

1502 # rotated companion matrix reduces error 

1503 m = legcompanion(c)[::-1,::-1] 

1504 r = la.eigvals(m) 

1505 r.sort() 

1506 return r 

1507 

1508 

1509def leggauss(deg): 

1510 """ 

1511 Gauss-Legendre quadrature. 

1512 

1513 Computes the sample points and weights for Gauss-Legendre quadrature. 

1514 These sample points and weights will correctly integrate polynomials of 

1515 degree :math:`2*deg - 1` or less over the interval :math:`[-1, 1]` with 

1516 the weight function :math:`f(x) = 1`. 

1517 

1518 Parameters 

1519 ---------- 

1520 deg : int 

1521 Number of sample points and weights. It must be >= 1. 

1522 

1523 Returns 

1524 ------- 

1525 x : ndarray 

1526 1-D ndarray containing the sample points. 

1527 y : ndarray 

1528 1-D ndarray containing the weights. 

1529 

1530 Notes 

1531 ----- 

1532 

1533 .. versionadded:: 1.7.0 

1534 

1535 The results have only been tested up to degree 100, higher degrees may 

1536 be problematic. The weights are determined by using the fact that 

1537 

1538 .. math:: w_k = c / (L'_n(x_k) * L_{n-1}(x_k)) 

1539 

1540 where :math:`c` is a constant independent of :math:`k` and :math:`x_k` 

1541 is the k'th root of :math:`L_n`, and then scaling the results to get 

1542 the right value when integrating 1. 

1543 

1544 """ 

1545 ideg = pu._deprecate_as_int(deg, "deg") 

1546 if ideg <= 0: 

1547 raise ValueError("deg must be a positive integer") 

1548 

1549 # first approximation of roots. We use the fact that the companion 

1550 # matrix is symmetric in this case in order to obtain better zeros. 

1551 c = np.array([0]*deg + [1]) 

1552 m = legcompanion(c) 

1553 x = la.eigvalsh(m) 

1554 

1555 # improve roots by one application of Newton 

1556 dy = legval(x, c) 

1557 df = legval(x, legder(c)) 

1558 x -= dy/df 

1559 

1560 # compute the weights. We scale the factor to avoid possible numerical 

1561 # overflow. 

1562 fm = legval(x, c[1:]) 

1563 fm /= np.abs(fm).max() 

1564 df /= np.abs(df).max() 

1565 w = 1/(fm * df) 

1566 

1567 # for Legendre we can also symmetrize 

1568 w = (w + w[::-1])/2 

1569 x = (x - x[::-1])/2 

1570 

1571 # scale w to get the right value 

1572 w *= 2. / w.sum() 

1573 

1574 return x, w 

1575 

1576 

1577def legweight(x): 

1578 """ 

1579 Weight function of the Legendre polynomials. 

1580 

1581 The weight function is :math:`1` and the interval of integration is 

1582 :math:`[-1, 1]`. The Legendre polynomials are orthogonal, but not 

1583 normalized, with respect to this weight function. 

1584 

1585 Parameters 

1586 ---------- 

1587 x : array_like 

1588 Values at which the weight function will be computed. 

1589 

1590 Returns 

1591 ------- 

1592 w : ndarray 

1593 The weight function at `x`. 

1594 

1595 Notes 

1596 ----- 

1597 

1598 .. versionadded:: 1.7.0 

1599 

1600 """ 

1601 w = x*0.0 + 1.0 

1602 return w 

1603 

1604# 

1605# Legendre series class 

1606# 

1607 

1608class Legendre(ABCPolyBase): 

1609 """A Legendre series class. 

1610 

1611 The Legendre class provides the standard Python numerical methods 

1612 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the 

1613 attributes and methods listed in the `ABCPolyBase` documentation. 

1614 

1615 Parameters 

1616 ---------- 

1617 coef : array_like 

1618 Legendre coefficients in order of increasing degree, i.e., 

1619 ``(1, 2, 3)`` gives ``1*P_0(x) + 2*P_1(x) + 3*P_2(x)``. 

1620 domain : (2,) array_like, optional 

1621 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped 

1622 to the interval ``[window[0], window[1]]`` by shifting and scaling. 

1623 The default value is [-1, 1]. 

1624 window : (2,) array_like, optional 

1625 Window, see `domain` for its use. The default value is [-1, 1]. 

1626 

1627 .. versionadded:: 1.6.0 

1628 

1629 """ 

1630 # Virtual Functions 

1631 _add = staticmethod(legadd) 

1632 _sub = staticmethod(legsub) 

1633 _mul = staticmethod(legmul) 

1634 _div = staticmethod(legdiv) 

1635 _pow = staticmethod(legpow) 

1636 _val = staticmethod(legval) 

1637 _int = staticmethod(legint) 

1638 _der = staticmethod(legder) 

1639 _fit = staticmethod(legfit) 

1640 _line = staticmethod(legline) 

1641 _roots = staticmethod(legroots) 

1642 _fromroots = staticmethod(legfromroots) 

1643 

1644 # Virtual properties 

1645 nickname = 'leg' 

1646 domain = np.array(legdomain) 

1647 window = np.array(legdomain) 

1648 basis_name = 'P'