Coverage for /usr/lib/python3/dist-packages/mpmath/calculus/optimization.py: 12%

1from __future__ import print_function 

2 

3from copy import copy 

4 

5from ..libmp.backend import xrange 

6 

7class OptimizationMethods(object): 

8 def __init__(ctx): 

9 pass 

10 

11############## 

12# 1D-SOLVERS # 

13############## 

14 

15class Newton: 

16 """ 

17 1d-solver generating pairs of approximative root and error. 

18 

19 Needs starting points x0 close to the root. 

20 

21 Pro: 

22 

23 * converges fast 

24 * sometimes more robust than secant with bad second starting point 

25 

26 Contra: 

27 

28 * converges slowly for multiple roots 

29 * needs first derivative 

30 * 2 function evaluations per iteration 

31 """ 

32 maxsteps = 20 

33 

34 def __init__(self, ctx, f, x0, **kwargs): 

35 self.ctx = ctx 

36 if len(x0) == 1: 

37 self.x0 = x0[0] 

38 else: 

39 raise ValueError('expected 1 starting point, got %i' % len(x0)) 

40 self.f = f 

41 if not 'df' in kwargs: 

42 def df(x): 

43 return self.ctx.diff(f, x) 

44 else: 

45 df = kwargs['df'] 

46 self.df = df 

47 

48 def __iter__(self): 

49 f = self.f 

50 df = self.df 

51 x0 = self.x0 

52 while True: 

53 x1 = x0 - f(x0) / df(x0) 

54 error = abs(x1 - x0) 

55 x0 = x1 

56 yield (x1, error) 

57 

58class Secant: 

59 """ 

60 1d-solver generating pairs of approximative root and error. 

61 

62 Needs starting points x0 and x1 close to the root. 

63 x1 defaults to x0 + 0.25. 

64 

65 Pro: 

66 

67 * converges fast 

68 

69 Contra: 

70 

71 * converges slowly for multiple roots 

72 """ 

73 maxsteps = 30 

74 

75 def __init__(self, ctx, f, x0, **kwargs): 

76 self.ctx = ctx 

77 if len(x0) == 1: 

78 self.x0 = x0[0] 

79 self.x1 = self.x0 + 0.25 

80 elif len(x0) == 2: 

81 self.x0 = x0[0] 

82 self.x1 = x0[1] 

83 else: 

84 raise ValueError('expected 1 or 2 starting points, got %i' % len(x0)) 

85 self.f = f 

86 

87 def __iter__(self): 

88 f = self.f 

89 x0 = self.x0 

90 x1 = self.x1 

91 f0 = f(x0) 

92 while True: 

93 f1 = f(x1) 

94 l = x1 - x0 

95 if not l: 

96 break 

97 s = (f1 - f0) / l 

98 if not s: 

99 break 

100 x0, x1 = x1, x1 - f1/s 

101 f0 = f1 

102 yield x1, abs(l) 

103 

104class MNewton: 

105 """ 

106 1d-solver generating pairs of approximative root and error. 

107 

108 Needs starting point x0 close to the root. 

109 Uses modified Newton's method that converges fast regardless of the 

110 multiplicity of the root. 

111 

112 Pro: 

113 

114 * converges fast for multiple roots 

115 

116 Contra: 

117 

118 * needs first and second derivative of f 

119 * 3 function evaluations per iteration 

120 """ 

121 maxsteps = 20 

122 

123 def __init__(self, ctx, f, x0, **kwargs): 

124 self.ctx = ctx 

125 if not len(x0) == 1: 

126 raise ValueError('expected 1 starting point, got %i' % len(x0)) 

127 self.x0 = x0[0] 

128 self.f = f 

129 if not 'df' in kwargs: 

130 def df(x): 

131 return self.ctx.diff(f, x) 

132 else: 

133 df = kwargs['df'] 

134 self.df = df 

135 if not 'd2f' in kwargs: 

136 def d2f(x): 

137 return self.ctx.diff(df, x) 

138 else: 

139 d2f = kwargs['df'] 

140 self.d2f = d2f 

141 

142 def __iter__(self): 

143 x = self.x0 

144 f = self.f 

145 df = self.df 

146 d2f = self.d2f 

147 while True: 

148 prevx = x 

149 fx = f(x) 

150 if fx == 0: 

151 break 

152 dfx = df(x) 

153 d2fx = d2f(x) 

154 # x = x - F(x)/F'(x) with F(x) = f(x)/f'(x) 

155 x -= fx / (dfx - fx * d2fx / dfx) 

156 error = abs(x - prevx) 

157 yield x, error 

158 

159class Halley: 

160 """ 

161 1d-solver generating pairs of approximative root and error. 

162 

163 Needs a starting point x0 close to the root. 

164 Uses Halley's method with cubic convergence rate. 

165 

166 Pro: 

167 

168 * converges even faster the Newton's method 

169 * useful when computing with *many* digits 

170 

171 Contra: 

172 

173 * needs first and second derivative of f 

174 * 3 function evaluations per iteration 

175 * converges slowly for multiple roots 

176 """ 

177 

178 maxsteps = 20 

179 

180 def __init__(self, ctx, f, x0, **kwargs): 

181 self.ctx = ctx 

182 if not len(x0) == 1: 

183 raise ValueError('expected 1 starting point, got %i' % len(x0)) 

184 self.x0 = x0[0] 

185 self.f = f 

186 if not 'df' in kwargs: 

187 def df(x): 

188 return self.ctx.diff(f, x) 

189 else: 

190 df = kwargs['df'] 

191 self.df = df 

192 if not 'd2f' in kwargs: 

193 def d2f(x): 

194 return self.ctx.diff(df, x) 

195 else: 

196 d2f = kwargs['df'] 

197 self.d2f = d2f 

198 

199 def __iter__(self): 

200 x = self.x0 

201 f = self.f 

202 df = self.df 

203 d2f = self.d2f 

204 while True: 

205 prevx = x 

206 fx = f(x) 

207 dfx = df(x) 

208 d2fx = d2f(x) 

209 x -= 2*fx*dfx / (2*dfx**2 - fx*d2fx) 

210 error = abs(x - prevx) 

211 yield x, error 

212 

213class Muller: 

214 """ 

215 1d-solver generating pairs of approximative root and error. 

216 

217 Needs starting points x0, x1 and x2 close to the root. 

218 x1 defaults to x0 + 0.25; x2 to x1 + 0.25. 

219 Uses Muller's method that converges towards complex roots. 

220 

221 Pro: 

222 

223 * converges fast (somewhat faster than secant) 

224 * can find complex roots 

225 

226 Contra: 

227 

228 * converges slowly for multiple roots 

229 * may have complex values for real starting points and real roots 

230 

231 http://en.wikipedia.org/wiki/Muller's_method 

232 """ 

233 maxsteps = 30 

234 

235 def __init__(self, ctx, f, x0, **kwargs): 

236 self.ctx = ctx 

237 if len(x0) == 1: 

238 self.x0 = x0[0] 

239 self.x1 = self.x0 + 0.25 

240 self.x2 = self.x1 + 0.25 

241 elif len(x0) == 2: 

242 self.x0 = x0[0] 

243 self.x1 = x0[1] 

244 self.x2 = self.x1 + 0.25 

245 elif len(x0) == 3: 

246 self.x0 = x0[0] 

247 self.x1 = x0[1] 

248 self.x2 = x0[2] 

249 else: 

250 raise ValueError('expected 1, 2 or 3 starting points, got %i' 

251 % len(x0)) 

252 self.f = f 

253 self.verbose = kwargs['verbose'] 

254 

255 def __iter__(self): 

256 f = self.f 

257 x0 = self.x0 

258 x1 = self.x1 

259 x2 = self.x2 

260 fx0 = f(x0) 

261 fx1 = f(x1) 

262 fx2 = f(x2) 

263 while True: 

264 # TODO: maybe refactoring with function for divided differences 

265 # calculate divided differences 

266 fx2x1 = (fx1 - fx2) / (x1 - x2) 

267 fx2x0 = (fx0 - fx2) / (x0 - x2) 

268 fx1x0 = (fx0 - fx1) / (x0 - x1) 

269 w = fx2x1 + fx2x0 - fx1x0 

270 fx2x1x0 = (fx1x0 - fx2x1) / (x0 - x2) 

271 if w == 0 and fx2x1x0 == 0: 

272 if self.verbose: 

273 print('canceled with') 

274 print('x0 =', x0, ', x1 =', x1, 'and x2 =', x2) 

275 break 

276 x0 = x1 

277 fx0 = fx1 

278 x1 = x2 

279 fx1 = fx2 

280 # denominator should be as large as possible => choose sign 

281 r = self.ctx.sqrt(w**2 - 4*fx2*fx2x1x0) 

282 if abs(w - r) > abs(w + r): 

283 r = -r 

284 x2 -= 2*fx2 / (w + r) 

285 fx2 = f(x2) 

286 error = abs(x2 - x1) 

287 yield x2, error 

288 

289# TODO: consider raising a ValueError when there's no sign change in a and b 

290class Bisection: 

291 """ 

292 1d-solver generating pairs of approximative root and error. 

293 

294 Uses bisection method to find a root of f in [a, b]. 

295 Might fail for multiple roots (needs sign change). 

296 

297 Pro: 

298 

299 * robust and reliable 

300 

301 Contra: 

302 

303 * converges slowly 

304 * needs sign change 

305 """ 

306 maxsteps = 100 

307 

308 def __init__(self, ctx, f, x0, **kwargs): 

309 self.ctx = ctx 

310 if len(x0) != 2: 

311 raise ValueError('expected interval of 2 points, got %i' % len(x0)) 

312 self.f = f 

313 self.a = x0[0] 

314 self.b = x0[1] 

315 

316 def __iter__(self): 

317 f = self.f 

318 a = self.a 

319 b = self.b 

320 l = b - a 

321 fb = f(b) 

322 while True: 

323 m = self.ctx.ldexp(a + b, -1) 

324 fm = f(m) 

325 sign = fm * fb 

326 if sign < 0: 

327 a = m 

328 elif sign > 0: 

329 b = m 

330 fb = fm 

331 else: 

332 yield m, self.ctx.zero 

333 l /= 2 

334 yield (a + b)/2, abs(l) 

335 

336def _getm(method): 

337 """ 

338 Return a function to calculate m for Illinois-like methods. 

339 """ 

340 if method == 'illinois': 

341 def getm(fz, fb): 

342 return 0.5 

343 elif method == 'pegasus': 

344 def getm(fz, fb): 

345 return fb/(fb + fz) 

346 elif method == 'anderson': 

347 def getm(fz, fb): 

348 m = 1 - fz/fb 

349 if m > 0: 

350 return m 

351 else: 

352 return 0.5 

353 else: 

354 raise ValueError("method '%s' not recognized" % method) 

355 return getm 

356 

357class Illinois: 

358 """ 

359 1d-solver generating pairs of approximative root and error. 

360 

361 Uses Illinois method or similar to find a root of f in [a, b]. 

362 Might fail for multiple roots (needs sign change). 

363 Combines bisect with secant (improved regula falsi). 

364 

365 The only difference between the methods is the scaling factor m, which is 

366 used to ensure convergence (you can choose one using the 'method' keyword): 

367 

368 Illinois method ('illinois'): 

369 m = 0.5 

370 

371 Pegasus method ('pegasus'): 

372 m = fb/(fb + fz) 

373 

374 Anderson-Bjoerk method ('anderson'): 

375 m = 1 - fz/fb if positive else 0.5 

376 

377 Pro: 

378 

379 * converges very fast 

380 

381 Contra: 

382 

383 * has problems with multiple roots 

384 * needs sign change 

385 """ 

386 maxsteps = 30 

387 

388 def __init__(self, ctx, f, x0, **kwargs): 

389 self.ctx = ctx 

390 if len(x0) != 2: 

391 raise ValueError('expected interval of 2 points, got %i' % len(x0)) 

392 self.a = x0[0] 

393 self.b = x0[1] 

394 self.f = f 

395 self.tol = kwargs['tol'] 

396 self.verbose = kwargs['verbose'] 

397 self.method = kwargs.get('method', 'illinois') 

398 self.getm = _getm(self.method) 

399 if self.verbose: 

400 print('using %s method' % self.method) 

401 

402 def __iter__(self): 

403 method = self.method 

404 f = self.f 

405 a = self.a 

406 b = self.b 

407 fa = f(a) 

408 fb = f(b) 

409 m = None 

410 while True: 

411 l = b - a 

412 if l == 0: 

413 break 

414 s = (fb - fa) / l 

415 z = a - fa/s 

416 fz = f(z) 

417 if abs(fz) < self.tol: 

418 # TODO: better condition (when f is very flat) 

419 if self.verbose: 

420 print('canceled with z =', z) 

421 yield z, l 

422 break 

423 if fz * fb < 0: # root in [z, b] 

424 a = b 

425 fa = fb 

426 b = z 

427 fb = fz 

428 else: # root in [a, z] 

429 m = self.getm(fz, fb) 

430 b = z 

431 fb = fz 

432 fa = m*fa # scale down to ensure convergence 

433 if self.verbose and m and not method == 'illinois': 

434 print('m:', m) 

435 yield (a + b)/2, abs(l) 

436 

437def Pegasus(*args, **kwargs): 

438 """ 

439 1d-solver generating pairs of approximative root and error. 

440 

441 Uses Pegasus method to find a root of f in [a, b]. 

442 Wrapper for illinois to use method='pegasus'. 

443 """ 

444 kwargs['method'] = 'pegasus' 

445 return Illinois(*args, **kwargs) 

446 

447def Anderson(*args, **kwargs): 

448 """ 

449 1d-solver generating pairs of approximative root and error. 

450 

451 Uses Anderson-Bjoerk method to find a root of f in [a, b]. 

452 Wrapper for illinois to use method='pegasus'. 

453 """ 

454 kwargs['method'] = 'anderson' 

455 return Illinois(*args, **kwargs) 

456 

457# TODO: check whether it's possible to combine it with Illinois stuff 

458class Ridder: 

459 """ 

460 1d-solver generating pairs of approximative root and error. 

461 

462 Ridders' method to find a root of f in [a, b]. 

463 Is told to perform as well as Brent's method while being simpler. 

464 

465 Pro: 

466 

467 * very fast 

468 * simpler than Brent's method 

469 

470 Contra: 

471 

472 * two function evaluations per step 

473 * has problems with multiple roots 

474 * needs sign change 

475 

476 http://en.wikipedia.org/wiki/Ridders'_method 

477 """ 

478 maxsteps = 30 

479 

480 def __init__(self, ctx, f, x0, **kwargs): 

481 self.ctx = ctx 

482 self.f = f 

483 if len(x0) != 2: 

484 raise ValueError('expected interval of 2 points, got %i' % len(x0)) 

485 self.x1 = x0[0] 

486 self.x2 = x0[1] 

487 self.verbose = kwargs['verbose'] 

488 self.tol = kwargs['tol'] 

489 

490 def __iter__(self): 

491 ctx = self.ctx 

492 f = self.f 

493 x1 = self.x1 

494 fx1 = f(x1) 

495 x2 = self.x2 

496 fx2 = f(x2) 

497 while True: 

498 x3 = 0.5*(x1 + x2) 

499 fx3 = f(x3) 

500 x4 = x3 + (x3 - x1) * ctx.sign(fx1 - fx2) * fx3 / ctx.sqrt(fx3**2 - fx1*fx2) 

501 fx4 = f(x4) 

502 if abs(fx4) < self.tol: 

503 # TODO: better condition (when f is very flat) 

504 if self.verbose: 

505 print('canceled with f(x4) =', fx4) 

506 yield x4, abs(x1 - x2) 

507 break 

508 if fx4 * fx2 < 0: # root in [x4, x2] 

509 x1 = x4 

510 fx1 = fx4 

511 else: # root in [x1, x4] 

512 x2 = x4 

513 fx2 = fx4 

514 error = abs(x1 - x2) 

515 yield (x1 + x2)/2, error 

516 

517class ANewton: 

518 """ 

519 EXPERIMENTAL 1d-solver generating pairs of approximative root and error. 

520 

521 Uses Newton's method modified to use Steffensens method when convergence is 

522 slow. (I.e. for multiple roots.) 

523 """ 

524 maxsteps = 20 

525 

526 def __init__(self, ctx, f, x0, **kwargs): 

527 self.ctx = ctx 

528 if not len(x0) == 1: 

529 raise ValueError('expected 1 starting point, got %i' % len(x0)) 

530 self.x0 = x0[0] 

531 self.f = f 

532 if not 'df' in kwargs: 

533 def df(x): 

534 return self.ctx.diff(f, x) 

535 else: 

536 df = kwargs['df'] 

537 self.df = df 

538 def phi(x): 

539 return x - f(x) / df(x) 

540 self.phi = phi 

541 self.verbose = kwargs['verbose'] 

542 

543 def __iter__(self): 

544 x0 = self.x0 

545 f = self.f 

546 df = self.df 

547 phi = self.phi 

548 error = 0 

549 counter = 0 

550 while True: 

551 prevx = x0 

552 try: 

553 x0 = phi(x0) 

554 except ZeroDivisionError: 

555 if self.verbose: 

556 print('ZeroDivisionError: canceled with x =', x0) 

557 break 

558 preverror = error 

559 error = abs(prevx - x0) 

560 # TODO: decide not to use convergence acceleration 

561 if error and abs(error - preverror) / error < 1: 

562 if self.verbose: 

563 print('converging slowly') 

564 counter += 1 

565 if counter >= 3: 

566 # accelerate convergence 

567 phi = steffensen(phi) 

568 counter = 0 

569 if self.verbose: 

570 print('accelerating convergence') 

571 yield x0, error 

572 

573# TODO: add Brent 

574 

575############################ 

576# MULTIDIMENSIONAL SOLVERS # 

577############################ 

578 

579def jacobian(ctx, f, x): 

580 """ 

581 Calculate the Jacobian matrix of a function at the point x0. 

582 

583 This is the first derivative of a vectorial function: 

584 

585 f : R^m -> R^n with m >= n 

586 """ 

587 x = ctx.matrix(x) 

588 h = ctx.sqrt(ctx.eps) 

589 fx = ctx.matrix(f(*x)) 

590 m = len(fx) 

591 n = len(x) 

592 J = ctx.matrix(m, n) 

593 for j in xrange(n): 

594 xj = x.copy() 

595 xj[j] += h 

596 Jj = (ctx.matrix(f(*xj)) - fx) / h 

597 for i in xrange(m): 

598 J[i,j] = Jj[i] 

599 return J 

600 

601# TODO: test with user-specified jacobian matrix 

602class MDNewton: 

603 """ 

604 Find the root of a vector function numerically using Newton's method. 

605 

606 f is a vector function representing a nonlinear equation system. 

607 

608 x0 is the starting point close to the root. 

609 

610 J is a function returning the Jacobian matrix for a point. 

611 

612 Supports overdetermined systems. 

613 

614 Use the 'norm' keyword to specify which norm to use. Defaults to max-norm. 

615 The function to calculate the Jacobian matrix can be given using the 

616 keyword 'J'. Otherwise it will be calculated numerically. 

617 

618 Please note that this method converges only locally. Especially for high- 

619 dimensional systems it is not trivial to find a good starting point being 

620 close enough to the root. 

621 

622 It is recommended to use a faster, low-precision solver from SciPy [1] or 

623 OpenOpt [2] to get an initial guess. Afterwards you can use this method for 

624 root-polishing to any precision. 

625 

626 [1] http://scipy.org 

627 

628 [2] http://openopt.org/Welcome 

629 """ 

630 maxsteps = 10 

631 

632 def __init__(self, ctx, f, x0, **kwargs): 

633 self.ctx = ctx 

634 self.f = f 

635 if isinstance(x0, (tuple, list)): 

636 x0 = ctx.matrix(x0) 

637 assert x0.cols == 1, 'need a vector' 

638 self.x0 = x0 

639 if 'J' in kwargs: 

640 self.J = kwargs['J'] 

641 else: 

642 def J(*x): 

643 return ctx.jacobian(f, x) 

644 self.J = J 

645 self.norm = kwargs['norm'] 

646 self.verbose = kwargs['verbose'] 

647 

648 def __iter__(self): 

649 f = self.f 

650 x0 = self.x0 

651 norm = self.norm 

652 J = self.J 

653 fx = self.ctx.matrix(f(*x0)) 

654 fxnorm = norm(fx) 

655 cancel = False 

656 while not cancel: 

657 # get direction of descent 

658 fxn = -fx 

659 Jx = J(*x0) 

660 s = self.ctx.lu_solve(Jx, fxn) 

661 if self.verbose: 

662 print('Jx:') 

663 print(Jx) 

664 print('s:', s) 

665 # damping step size TODO: better strategy (hard task) 

666 l = self.ctx.one 

667 x1 = x0 + s 

668 while True: 

669 if x1 == x0: 

670 if self.verbose: 

671 print("canceled, won't get more excact") 

672 cancel = True 

673 break 

674 fx = self.ctx.matrix(f(*x1)) 

675 newnorm = norm(fx) 

676 if newnorm < fxnorm: 

677 # new x accepted 

678 fxnorm = newnorm 

679 x0 = x1 

680 break 

681 l /= 2 

682 x1 = x0 + l*s 

683 yield (x0, fxnorm) 

684 

685############# 

686# UTILITIES # 

687############# 

688 

689str2solver = {'newton':Newton, 'secant':Secant, 'mnewton':MNewton, 

690 'halley':Halley, 'muller':Muller, 'bisect':Bisection, 

691 'illinois':Illinois, 'pegasus':Pegasus, 'anderson':Anderson, 

692 'ridder':Ridder, 'anewton':ANewton, 'mdnewton':MDNewton} 

693 

694def findroot(ctx, f, x0, solver='secant', tol=None, verbose=False, verify=True, **kwargs): 

695 r""" 

696 Find an approximate solution to `f(x) = 0`, using *x0* as starting point or 

697 interval for *x*. 

698 

699 Multidimensional overdetermined systems are supported. 

700 You can specify them using a function or a list of functions. 

701 

702 Mathematically speaking, this function returns `x` such that 

703 `|f(x)|^2 \leq \mathrm{tol}` is true within the current working precision. 

704 If the computed value does not meet this criterion, an exception is raised. 

705 This exception can be disabled with *verify=False*. 

706 

707 For interval arithmetic (``iv.findroot()``), please note that 

708 the returned interval ``x`` is not guaranteed to contain `f(x)=0`! 

709 It is only some `x` for which `|f(x)|^2 \leq \mathrm{tol}` certainly holds 

710 regardless of numerical error. This may be improved in the future. 

711 

712 **Arguments** 

713 

714 *f* 

715 one dimensional function 

716 *x0* 

717 starting point, several starting points or interval (depends on solver) 

718 *tol* 

719 the returned solution has an error smaller than this 

720 *verbose* 

721 print additional information for each iteration if true 

722 *verify* 

723 verify the solution and raise a ValueError if `|f(x)|^2 > \mathrm{tol}` 

724 *solver* 

725 a generator for *f* and *x0* returning approximative solution and error 

726 *maxsteps* 

727 after how many steps the solver will cancel 

728 *df* 

729 first derivative of *f* (used by some solvers) 

730 *d2f* 

731 second derivative of *f* (used by some solvers) 

732 *multidimensional* 

733 force multidimensional solving 

734 *J* 

735 Jacobian matrix of *f* (used by multidimensional solvers) 

736 *norm* 

737 used vector norm (used by multidimensional solvers) 

738 

739 solver has to be callable with ``(f, x0, **kwargs)`` and return an generator 

740 yielding pairs of approximative solution and estimated error (which is 

741 expected to be positive). 

742 You can use the following string aliases: 

743 'secant', 'mnewton', 'halley', 'muller', 'illinois', 'pegasus', 'anderson', 

744 'ridder', 'anewton', 'bisect' 

745 

746 See mpmath.calculus.optimization for their documentation. 

747 

748 **Examples** 

749 

750 The function :func:`~mpmath.findroot` locates a root of a given function using the 

751 secant method by default. A simple example use of the secant method is to 

752 compute `\pi` as the root of `\sin x` closest to `x_0 = 3`:: 

753 

754 >>> from mpmath import * 

755 >>> mp.dps = 30; mp.pretty = True 

756 >>> findroot(sin, 3) 

757 3.14159265358979323846264338328 

758 

759 The secant method can be used to find complex roots of analytic functions, 

760 although it must in that case generally be given a nonreal starting value 

761 (or else it will never leave the real line):: 

762 

763 >>> mp.dps = 15 

764 >>> findroot(lambda x: x**3 + 2*x + 1, j) 

765 (0.226698825758202 + 1.46771150871022j) 

766 

767 A nice application is to compute nontrivial roots of the Riemann zeta 

768 function with many digits (good initial values are needed for convergence):: 

769 

770 >>> mp.dps = 30 

771 >>> findroot(zeta, 0.5+14j) 

772 (0.5 + 14.1347251417346937904572519836j) 

773 

774 The secant method can also be used as an optimization algorithm, by passing 

775 it a derivative of a function. The following example locates the positive 

776 minimum of the gamma function:: 

777 

778 >>> mp.dps = 20 

779 >>> findroot(lambda x: diff(gamma, x), 1) 

780 1.4616321449683623413 

781 

782 Finally, a useful application is to compute inverse functions, such as the 

783 Lambert W function which is the inverse of `w e^w`, given the first 

784 term of the solution's asymptotic expansion as the initial value. In basic 

785 cases, this gives identical results to mpmath's built-in ``lambertw`` 

786 function:: 

787 

788 >>> def lambert(x): 

789 ... return findroot(lambda w: w*exp(w) - x, log(1+x)) 

790 ... 

791 >>> mp.dps = 15 

792 >>> lambert(1); lambertw(1) 

793 0.567143290409784 

794 0.567143290409784 

795 >>> lambert(1000); lambert(1000) 

796 5.2496028524016 

797 5.2496028524016 

798 

799 Multidimensional functions are also supported:: 

800 

801 >>> f = [lambda x1, x2: x1**2 + x2, 

802 ... lambda x1, x2: 5*x1**2 - 3*x1 + 2*x2 - 3] 

803 >>> findroot(f, (0, 0)) 

804 [-0.618033988749895] 

805 [-0.381966011250105] 

806 >>> findroot(f, (10, 10)) 

807 [ 1.61803398874989] 

808 [-2.61803398874989] 

809 

810 You can verify this by solving the system manually. 

811 

812 Please note that the following (more general) syntax also works:: 

813 

814 >>> def f(x1, x2): 

815 ... return x1**2 + x2, 5*x1**2 - 3*x1 + 2*x2 - 3 

816 ... 

817 >>> findroot(f, (0, 0)) 

818 [-0.618033988749895] 

819 [-0.381966011250105] 

820 

821 

822 **Multiple roots** 

823 

824 For multiple roots all methods of the Newtonian family (including secant) 

825 converge slowly. Consider this example:: 

826 

827 >>> f = lambda x: (x - 1)**99 

828 >>> findroot(f, 0.9, verify=False) 

829 0.918073542444929 

830 

831 Even for a very close starting point the secant method converges very 

832 slowly. Use ``verbose=True`` to illustrate this. 

833 

834 It is possible to modify Newton's method to make it converge regardless of 

835 the root's multiplicity:: 

836 

837 >>> findroot(f, -10, solver='mnewton') 

838 1.0 

839 

840 This variant uses the first and second derivative of the function, which is 

841 not very efficient. 

842 

843 Alternatively you can use an experimental Newtonian solver that keeps track 

844 of the speed of convergence and accelerates it using Steffensen's method if 

845 necessary:: 

846 

847 >>> findroot(f, -10, solver='anewton', verbose=True) 

848 x: -9.88888888888888888889 

849 error: 0.111111111111111111111 

850 converging slowly 

851 x: -9.77890011223344556678 

852 error: 0.10998877665544332211 

853 converging slowly 

854 x: -9.67002233332199662166 

855 error: 0.108877778911448945119 

856 converging slowly 

857 accelerating convergence 

858 x: -9.5622443299551077669 

859 error: 0.107778003366888854764 

860 converging slowly 

861 x: 0.99999999999999999214 

862 error: 10.562244329955107759 

863 x: 1.0 

864 error: 7.8598304758094664213e-18 

865 ZeroDivisionError: canceled with x = 1.0 

866 1.0 

867 

868 **Complex roots** 

869 

870 For complex roots it's recommended to use Muller's method as it converges 

871 even for real starting points very fast:: 

872 

873 >>> findroot(lambda x: x**4 + x + 1, (0, 1, 2), solver='muller') 

874 (0.727136084491197 + 0.934099289460529j) 

875 

876 

877 **Intersection methods** 

878 

879 When you need to find a root in a known interval, it's highly recommended to 

880 use an intersection-based solver like ``'anderson'`` or ``'ridder'``. 

881 Usually they converge faster and more reliable. They have however problems 

882 with multiple roots and usually need a sign change to find a root:: 

883 

884 >>> findroot(lambda x: x**3, (-1, 1), solver='anderson') 

885 0.0 

886 

887 Be careful with symmetric functions:: 

888 

889 >>> findroot(lambda x: x**2, (-1, 1), solver='anderson') #doctest:+ELLIPSIS 

890 Traceback (most recent call last): 

891 ... 

892 ZeroDivisionError 

893 

894 It fails even for better starting points, because there is no sign change:: 

895 

896 >>> findroot(lambda x: x**2, (-1, .5), solver='anderson') 

897 Traceback (most recent call last): 

898 ... 

899 ValueError: Could not find root within given tolerance. (1.0 > 2.16840434497100886801e-19) 

900 Try another starting point or tweak arguments. 

901 

902 """ 

903 prec = ctx.prec 

904 try: 

905 ctx.prec += 20 

906 

907 # initialize arguments 

908 if tol is None: 

909 tol = ctx.eps * 2**10 

910 

911 kwargs['verbose'] = kwargs.get('verbose', verbose) 

912 

913 if 'd1f' in kwargs: 

914 kwargs['df'] = kwargs['d1f'] 

915 

916 kwargs['tol'] = tol 

917 if isinstance(x0, (list, tuple)): 

918 x0 = [ctx.convert(x) for x in x0] 

919 else: 

920 x0 = [ctx.convert(x0)] 

921 

922 if isinstance(solver, str): 

923 try: 

924 solver = str2solver[solver] 

925 except KeyError: 

926 raise ValueError('could not recognize solver') 

927 

928 # accept list of functions 

929 if isinstance(f, (list, tuple)): 

930 f2 = copy(f) 

931 def tmp(*args): 

932 return [fn(*args) for fn in f2] 

933 f = tmp 

934 

935 # detect multidimensional functions 

936 try: 

937 fx = f(*x0) 

938 multidimensional = isinstance(fx, (list, tuple, ctx.matrix)) 

939 except TypeError: 

940 fx = f(x0[0]) 

941 multidimensional = False 

942 if 'multidimensional' in kwargs: 

943 multidimensional = kwargs['multidimensional'] 

944 if multidimensional: 

945 # only one multidimensional solver available at the moment