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# Authors: Pearu Peterson, Pauli Virtanen, John Travers 

2""" 

3First-order ODE integrators. 

4 

5User-friendly interface to various numerical integrators for solving a 

6system of first order ODEs with prescribed initial conditions:: 

7 

8 d y(t)[i] 

9 --------- = f(t,y(t))[i], 

10 d t 

11 

12 y(t=0)[i] = y0[i], 

13 

14where:: 

15 

16 i = 0, ..., len(y0) - 1 

17 

18class ode 

19--------- 

20 

21A generic interface class to numeric integrators. It has the following 

22methods:: 

23 

24 integrator = ode(f, jac=None) 

25 integrator = integrator.set_integrator(name, **params) 

26 integrator = integrator.set_initial_value(y0, t0=0.0) 

27 integrator = integrator.set_f_params(*args) 

28 integrator = integrator.set_jac_params(*args) 

29 y1 = integrator.integrate(t1, step=False, relax=False) 

30 flag = integrator.successful() 

31 

32class complex_ode 

33----------------- 

34 

35This class has the same generic interface as ode, except it can handle complex 

36f, y and Jacobians by transparently translating them into the equivalent 

37real-valued system. It supports the real-valued solvers (i.e., not zvode) and is 

38an alternative to ode with the zvode solver, sometimes performing better. 

39""" 

40# XXX: Integrators must have: 

41# =========================== 

42# cvode - C version of vode and vodpk with many improvements. 

43# Get it from http://www.netlib.org/ode/cvode.tar.gz. 

44# To wrap cvode to Python, one must write the extension module by 

45# hand. Its interface is too much 'advanced C' that using f2py 

46# would be too complicated (or impossible). 

47# 

48# How to define a new integrator: 

49# =============================== 

50# 

51# class myodeint(IntegratorBase): 

52# 

53# runner = <odeint function> or None 

54# 

55# def __init__(self,...): # required 

56# <initialize> 

57# 

58# def reset(self,n,has_jac): # optional 

59# # n - the size of the problem (number of equations) 

60# # has_jac - whether user has supplied its own routine for Jacobian 

61# <allocate memory,initialize further> 

62# 

63# def run(self,f,jac,y0,t0,t1,f_params,jac_params): # required 

64# # this method is called to integrate from t=t0 to t=t1 

65# # with initial condition y0. f and jac are user-supplied functions 

66# # that define the problem. f_params,jac_params are additional 

67# # arguments 

68# # to these functions. 

69# <calculate y1> 

70# if <calculation was unsuccessful>: 

71# self.success = 0 

72# return t1,y1 

73# 

74# # In addition, one can define step() and run_relax() methods (they 

75# # take the same arguments as run()) if the integrator can support 

76# # these features (see IntegratorBase doc strings). 

77# 

78# if myodeint.runner: 

79# IntegratorBase.integrator_classes.append(myodeint) 

80 

81__all__ = ['ode', 'complex_ode'] 

82__version__ = "$Id$" 

83__docformat__ = "restructuredtext en" 

84 

85import re 

86import warnings 

87 

88from numpy import asarray, array, zeros, isscalar, real, imag, vstack 

89 

90from . import vode as _vode 

91from . import _dop 

92from . import lsoda as _lsoda 

93 

94 

95_dop_int_dtype = _dop.types.intvar.dtype 

96_vode_int_dtype = _vode.types.intvar.dtype 

97_lsoda_int_dtype = _lsoda.types.intvar.dtype 

98 

99 

100# ------------------------------------------------------------------------------ 

101# User interface 

102# ------------------------------------------------------------------------------ 

103 

104 

105class ode(object): 

106 """ 

107 A generic interface class to numeric integrators. 

108 

109 Solve an equation system :math:`y'(t) = f(t,y)` with (optional) ``jac = df/dy``. 

110 

111 *Note*: The first two arguments of ``f(t, y, ...)`` are in the 

112 opposite order of the arguments in the system definition function used 

113 by `scipy.integrate.odeint`. 

114 

115 Parameters 

116 ---------- 

117 f : callable ``f(t, y, *f_args)`` 

118 Right-hand side of the differential equation. t is a scalar, 

119 ``y.shape == (n,)``. 

120 ``f_args`` is set by calling ``set_f_params(*args)``. 

121 `f` should return a scalar, array or list (not a tuple). 

122 jac : callable ``jac(t, y, *jac_args)``, optional 

123 Jacobian of the right-hand side, ``jac[i,j] = d f[i] / d y[j]``. 

124 ``jac_args`` is set by calling ``set_jac_params(*args)``. 

125 

126 Attributes 

127 ---------- 

128 t : float 

129 Current time. 

130 y : ndarray 

131 Current variable values. 

132 

133 See also 

134 -------- 

135 odeint : an integrator with a simpler interface based on lsoda from ODEPACK 

136 quad : for finding the area under a curve 

137 

138 Notes 

139 ----- 

140 Available integrators are listed below. They can be selected using 

141 the `set_integrator` method. 

142 

143 "vode" 

144 

145 Real-valued Variable-coefficient Ordinary Differential Equation 

146 solver, with fixed-leading-coefficient implementation. It provides 

147 implicit Adams method (for non-stiff problems) and a method based on 

148 backward differentiation formulas (BDF) (for stiff problems). 

149 

150 Source: http://www.netlib.org/ode/vode.f 

151 

152 .. warning:: 

153 

154 This integrator is not re-entrant. You cannot have two `ode` 

155 instances using the "vode" integrator at the same time. 

156 

157 This integrator accepts the following parameters in `set_integrator` 

158 method of the `ode` class: 

159 

160 - atol : float or sequence 

161 absolute tolerance for solution 

162 - rtol : float or sequence 

163 relative tolerance for solution 

164 - lband : None or int 

165 - uband : None or int 

166 Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband. 

167 Setting these requires your jac routine to return the jacobian 

168 in packed format, jac_packed[i-j+uband, j] = jac[i,j]. The 

169 dimension of the matrix must be (lband+uband+1, len(y)). 

170 - method: 'adams' or 'bdf' 

171 Which solver to use, Adams (non-stiff) or BDF (stiff) 

172 - with_jacobian : bool 

173 This option is only considered when the user has not supplied a 

174 Jacobian function and has not indicated (by setting either band) 

175 that the Jacobian is banded. In this case, `with_jacobian` specifies 

176 whether the iteration method of the ODE solver's correction step is 

177 chord iteration with an internally generated full Jacobian or 

178 functional iteration with no Jacobian. 

179 - nsteps : int 

180 Maximum number of (internally defined) steps allowed during one 

181 call to the solver. 

182 - first_step : float 

183 - min_step : float 

184 - max_step : float 

185 Limits for the step sizes used by the integrator. 

186 - order : int 

187 Maximum order used by the integrator, 

188 order <= 12 for Adams, <= 5 for BDF. 

189 

190 "zvode" 

191 

192 Complex-valued Variable-coefficient Ordinary Differential Equation 

193 solver, with fixed-leading-coefficient implementation. It provides 

194 implicit Adams method (for non-stiff problems) and a method based on 

195 backward differentiation formulas (BDF) (for stiff problems). 

196 

197 Source: http://www.netlib.org/ode/zvode.f 

198 

199 .. warning:: 

200 

201 This integrator is not re-entrant. You cannot have two `ode` 

202 instances using the "zvode" integrator at the same time. 

203 

204 This integrator accepts the same parameters in `set_integrator` 

205 as the "vode" solver. 

206 

207 .. note:: 

208 

209 When using ZVODE for a stiff system, it should only be used for 

210 the case in which the function f is analytic, that is, when each f(i) 

211 is an analytic function of each y(j). Analyticity means that the 

212 partial derivative df(i)/dy(j) is a unique complex number, and this 

213 fact is critical in the way ZVODE solves the dense or banded linear 

214 systems that arise in the stiff case. For a complex stiff ODE system 

215 in which f is not analytic, ZVODE is likely to have convergence 

216 failures, and for this problem one should instead use DVODE on the 

217 equivalent real system (in the real and imaginary parts of y). 

218 

219 "lsoda" 

220 

221 Real-valued Variable-coefficient Ordinary Differential Equation 

222 solver, with fixed-leading-coefficient implementation. It provides 

223 automatic method switching between implicit Adams method (for non-stiff 

224 problems) and a method based on backward differentiation formulas (BDF) 

225 (for stiff problems). 

226 

227 Source: http://www.netlib.org/odepack 

228 

229 .. warning:: 

230 

231 This integrator is not re-entrant. You cannot have two `ode` 

232 instances using the "lsoda" integrator at the same time. 

233 

234 This integrator accepts the following parameters in `set_integrator` 

235 method of the `ode` class: 

236 

237 - atol : float or sequence 

238 absolute tolerance for solution 

239 - rtol : float or sequence 

240 relative tolerance for solution 

241 - lband : None or int 

242 - uband : None or int 

243 Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband. 

244 Setting these requires your jac routine to return the jacobian 

245 in packed format, jac_packed[i-j+uband, j] = jac[i,j]. 

246 - with_jacobian : bool 

247 *Not used.* 

248 - nsteps : int 

249 Maximum number of (internally defined) steps allowed during one 

250 call to the solver. 

251 - first_step : float 

252 - min_step : float 

253 - max_step : float 

254 Limits for the step sizes used by the integrator. 

255 - max_order_ns : int 

256 Maximum order used in the nonstiff case (default 12). 

257 - max_order_s : int 

258 Maximum order used in the stiff case (default 5). 

259 - max_hnil : int 

260 Maximum number of messages reporting too small step size (t + h = t) 

261 (default 0) 

262 - ixpr : int 

263 Whether to generate extra printing at method switches (default False). 

264 

265 "dopri5" 

266 

267 This is an explicit runge-kutta method of order (4)5 due to Dormand & 

268 Prince (with stepsize control and dense output). 

269 

270 Authors: 

271 

272 E. Hairer and G. Wanner 

273 Universite de Geneve, Dept. de Mathematiques 

274 CH-1211 Geneve 24, Switzerland 

275 e-mail: ernst.hairer@math.unige.ch, gerhard.wanner@math.unige.ch 

276 

277 This code is described in [HNW93]_. 

278 

279 This integrator accepts the following parameters in set_integrator() 

280 method of the ode class: 

281 

282 - atol : float or sequence 

283 absolute tolerance for solution 

284 - rtol : float or sequence 

285 relative tolerance for solution 

286 - nsteps : int 

287 Maximum number of (internally defined) steps allowed during one 

288 call to the solver. 

289 - first_step : float 

290 - max_step : float 

291 - safety : float 

292 Safety factor on new step selection (default 0.9) 

293 - ifactor : float 

294 - dfactor : float 

295 Maximum factor to increase/decrease step size by in one step 

296 - beta : float 

297 Beta parameter for stabilised step size control. 

298 - verbosity : int 

299 Switch for printing messages (< 0 for no messages). 

300 

301 "dop853" 

302 

303 This is an explicit runge-kutta method of order 8(5,3) due to Dormand 

304 & Prince (with stepsize control and dense output). 

305 

306 Options and references the same as "dopri5". 

307 

308 Examples 

309 -------- 

310 

311 A problem to integrate and the corresponding jacobian: 

312 

313 >>> from scipy.integrate import ode 

314 >>> 

315 >>> y0, t0 = [1.0j, 2.0], 0 

316 >>> 

317 >>> def f(t, y, arg1): 

318 ... return [1j*arg1*y[0] + y[1], -arg1*y[1]**2] 

319 >>> def jac(t, y, arg1): 

320 ... return [[1j*arg1, 1], [0, -arg1*2*y[1]]] 

321 

322 The integration: 

323 

324 >>> r = ode(f, jac).set_integrator('zvode', method='bdf') 

325 >>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0) 

326 >>> t1 = 10 

327 >>> dt = 1 

328 >>> while r.successful() and r.t < t1: 

329 ... print(r.t+dt, r.integrate(r.t+dt)) 

330 1 [-0.71038232+0.23749653j 0.40000271+0.j ] 

331 2.0 [0.19098503-0.52359246j 0.22222356+0.j ] 

332 3.0 [0.47153208+0.52701229j 0.15384681+0.j ] 

333 4.0 [-0.61905937+0.30726255j 0.11764744+0.j ] 

334 5.0 [0.02340997-0.61418799j 0.09523835+0.j ] 

335 6.0 [0.58643071+0.339819j 0.08000018+0.j ] 

336 7.0 [-0.52070105+0.44525141j 0.06896565+0.j ] 

337 8.0 [-0.15986733-0.61234476j 0.06060616+0.j ] 

338 9.0 [0.64850462+0.15048982j 0.05405414+0.j ] 

339 10.0 [-0.38404699+0.56382299j 0.04878055+0.j ] 

340 

341 References 

342 ---------- 

343 .. [HNW93] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary 

344 Differential Equations i. Nonstiff Problems. 2nd edition. 

345 Springer Series in Computational Mathematics, 

346 Springer-Verlag (1993) 

347 

348 """ 

349 

350 def __init__(self, f, jac=None): 

351 self.stiff = 0 

352 self.f = f 

353 self.jac = jac 

354 self.f_params = () 

355 self.jac_params = () 

356 self._y = [] 

357 

358 @property 

359 def y(self): 

360 return self._y 

361 

362 def set_initial_value(self, y, t=0.0): 

363 """Set initial conditions y(t) = y.""" 

364 if isscalar(y): 

365 y = [y] 

366 n_prev = len(self._y) 

367 if not n_prev: 

368 self.set_integrator('') # find first available integrator 

369 self._y = asarray(y, self._integrator.scalar) 

370 self.t = t 

371 self._integrator.reset(len(self._y), self.jac is not None) 

372 return self 

373 

374 def set_integrator(self, name, **integrator_params): 

375 """ 

376 Set integrator by name. 

377 

378 Parameters 

379 ---------- 

380 name : str 

381 Name of the integrator. 

382 integrator_params 

383 Additional parameters for the integrator. 

384 """ 

385 integrator = find_integrator(name) 

386 if integrator is None: 

387 # FIXME: this really should be raise an exception. Will that break 

388 # any code? 

389 warnings.warn('No integrator name match with %r or is not ' 

390 'available.' % name) 

391 else: 

392 self._integrator = integrator(**integrator_params) 

393 if not len(self._y): 

394 self.t = 0.0 

395 self._y = array([0.0], self._integrator.scalar) 

396 self._integrator.reset(len(self._y), self.jac is not None) 

397 return self 

398 

399 def integrate(self, t, step=False, relax=False): 

400 """Find y=y(t), set y as an initial condition, and return y. 

401 

402 Parameters 

403 ---------- 

404 t : float 

405 The endpoint of the integration step. 

406 step : bool 

407 If True, and if the integrator supports the step method, 

408 then perform a single integration step and return. 

409 This parameter is provided in order to expose internals of 

410 the implementation, and should not be changed from its default 

411 value in most cases. 

412 relax : bool 

413 If True and if the integrator supports the run_relax method, 

414 then integrate until t_1 >= t and return. ``relax`` is not 

415 referenced if ``step=True``. 

416 This parameter is provided in order to expose internals of 

417 the implementation, and should not be changed from its default 

418 value in most cases. 

419 

420 Returns 

421 ------- 

422 y : float 

423 The integrated value at t 

424 """ 

425 if step and self._integrator.supports_step: 

426 mth = self._integrator.step 

427 elif relax and self._integrator.supports_run_relax: 

428 mth = self._integrator.run_relax 

429 else: 

430 mth = self._integrator.run 

431 

432 try: 

433 self._y, self.t = mth(self.f, self.jac or (lambda: None), 

434 self._y, self.t, t, 

435 self.f_params, self.jac_params) 

436 except SystemError: 

437 # f2py issue with tuple returns, see ticket 1187. 

438 raise ValueError('Function to integrate must not return a tuple.') 

439 

440 return self._y 

441 

442 def successful(self): 

443 """Check if integration was successful.""" 

444 try: 

445 self._integrator 

446 except AttributeError: 

447 self.set_integrator('') 

448 return self._integrator.success == 1 

449 

450 def get_return_code(self): 

451 """Extracts the return code for the integration to enable better control 

452 if the integration fails. 

453 

454 In general, a return code > 0 implies success, while a return code < 0 

455 implies failure. 

456 

457 Notes 

458 ----- 

459 This section describes possible return codes and their meaning, for available 

460 integrators that can be selected by `set_integrator` method. 

461 

462 "vode" 

463 

464 =========== ======= 

465 Return Code Message 

466 =========== ======= 

467 2 Integration successful. 

468 -1 Excess work done on this call. (Perhaps wrong MF.) 

469 -2 Excess accuracy requested. (Tolerances too small.) 

470 -3 Illegal input detected. (See printed message.) 

471 -4 Repeated error test failures. (Check all input.) 

472 -5 Repeated convergence failures. (Perhaps bad Jacobian 

473 supplied or wrong choice of MF or tolerances.) 

474 -6 Error weight became zero during problem. (Solution 

475 component i vanished, and ATOL or ATOL(i) = 0.) 

476 =========== ======= 

477 

478 "zvode" 

479 

480 =========== ======= 

481 Return Code Message 

482 =========== ======= 

483 2 Integration successful. 

484 -1 Excess work done on this call. (Perhaps wrong MF.) 

485 -2 Excess accuracy requested. (Tolerances too small.) 

486 -3 Illegal input detected. (See printed message.) 

487 -4 Repeated error test failures. (Check all input.) 

488 -5 Repeated convergence failures. (Perhaps bad Jacobian 

489 supplied or wrong choice of MF or tolerances.) 

490 -6 Error weight became zero during problem. (Solution 

491 component i vanished, and ATOL or ATOL(i) = 0.) 

492 =========== ======= 

493 

494 "dopri5" 

495 

496 =========== ======= 

497 Return Code Message 

498 =========== ======= 

499 1 Integration successful. 

500 2 Integration successful (interrupted by solout). 

501 -1 Input is not consistent. 

502 -2 Larger nsteps is needed. 

503 -3 Step size becomes too small. 

504 -4 Problem is probably stiff (interrupted). 

505 =========== ======= 

506 

507 "dop853" 

508 

509 =========== ======= 

510 Return Code Message 

511 =========== ======= 

512 1 Integration successful. 

513 2 Integration successful (interrupted by solout). 

514 -1 Input is not consistent. 

515 -2 Larger nsteps is needed. 

516 -3 Step size becomes too small. 

517 -4 Problem is probably stiff (interrupted). 

518 =========== ======= 

519 

520 "lsoda" 

521 

522 =========== ======= 

523 Return Code Message 

524 =========== ======= 

525 2 Integration successful. 

526 -1 Excess work done on this call (perhaps wrong Dfun type). 

527 -2 Excess accuracy requested (tolerances too small). 

528 -3 Illegal input detected (internal error). 

529 -4 Repeated error test failures (internal error). 

530 -5 Repeated convergence failures (perhaps bad Jacobian or tolerances). 

531 -6 Error weight became zero during problem. 

532 -7 Internal workspace insufficient to finish (internal error). 

533 =========== ======= 

534 """ 

535 try: 

536 self._integrator 

537 except AttributeError: 

538 self.set_integrator('') 

539 return self._integrator.istate 

540 

541 def set_f_params(self, *args): 

542 """Set extra parameters for user-supplied function f.""" 

543 self.f_params = args 

544 return self 

545 

546 def set_jac_params(self, *args): 

547 """Set extra parameters for user-supplied function jac.""" 

548 self.jac_params = args 

549 return self 

550 

551 def set_solout(self, solout): 

552 """ 

553 Set callable to be called at every successful integration step. 

554 

555 Parameters 

556 ---------- 

557 solout : callable 

558 ``solout(t, y)`` is called at each internal integrator step, 

559 t is a scalar providing the current independent position 

560 y is the current soloution ``y.shape == (n,)`` 

561 solout should return -1 to stop integration 

562 otherwise it should return None or 0 

563 

564 """ 

565 if self._integrator.supports_solout: 

566 self._integrator.set_solout(solout) 

567 if self._y is not None: 

568 self._integrator.reset(len(self._y), self.jac is not None) 

569 else: 

570 raise ValueError("selected integrator does not support solout," 

571 " choose another one") 

572 

573 

574def _transform_banded_jac(bjac): 

575 """ 

576 Convert a real matrix of the form (for example) 

577 

578 [0 0 A B] [0 0 0 B] 

579 [0 0 C D] [0 0 A D] 

580 [E F G H] to [0 F C H] 

581 [I J K L] [E J G L] 

582 [I 0 K 0] 

583 

584 That is, every other column is shifted up one. 

585 """ 

586 # Shift every other column. 

587 newjac = zeros((bjac.shape[0] + 1, bjac.shape[1])) 

588 newjac[1:, ::2] = bjac[:, ::2] 

589 newjac[:-1, 1::2] = bjac[:, 1::2] 

590 return newjac 

591 

592 

593class complex_ode(ode): 

594 """ 

595 A wrapper of ode for complex systems. 

596 

597 This functions similarly as `ode`, but re-maps a complex-valued 

598 equation system to a real-valued one before using the integrators. 

599 

600 Parameters 

601 ---------- 

602 f : callable ``f(t, y, *f_args)`` 

603 Rhs of the equation. t is a scalar, ``y.shape == (n,)``. 

604 ``f_args`` is set by calling ``set_f_params(*args)``. 

605 jac : callable ``jac(t, y, *jac_args)`` 

606 Jacobian of the rhs, ``jac[i,j] = d f[i] / d y[j]``. 

607 ``jac_args`` is set by calling ``set_f_params(*args)``. 

608 

609 Attributes 

610 ---------- 

611 t : float 

612 Current time. 

613 y : ndarray 

614 Current variable values. 

615 

616 Examples 

617 -------- 

618 For usage examples, see `ode`. 

619 

620 """ 

621 

622 def __init__(self, f, jac=None): 

623 self.cf = f 

624 self.cjac = jac 

625 if jac is None: 

626 ode.__init__(self, self._wrap, None) 

627 else: 

628 ode.__init__(self, self._wrap, self._wrap_jac) 

629 

630 def _wrap(self, t, y, *f_args): 

631 f = self.cf(*((t, y[::2] + 1j * y[1::2]) + f_args)) 

632 # self.tmp is a real-valued array containing the interleaved 

633 # real and imaginary parts of f. 

634 self.tmp[::2] = real(f) 

635 self.tmp[1::2] = imag(f) 

636 return self.tmp 

637 

638 def _wrap_jac(self, t, y, *jac_args): 

639 # jac is the complex Jacobian computed by the user-defined function. 

640 jac = self.cjac(*((t, y[::2] + 1j * y[1::2]) + jac_args)) 

641 

642 # jac_tmp is the real version of the complex Jacobian. Each complex 

643 # entry in jac, say 2+3j, becomes a 2x2 block of the form 

644 # [2 -3] 

645 # [3 2] 

646 jac_tmp = zeros((2 * jac.shape[0], 2 * jac.shape[1])) 

647 jac_tmp[1::2, 1::2] = jac_tmp[::2, ::2] = real(jac) 

648 jac_tmp[1::2, ::2] = imag(jac) 

649 jac_tmp[::2, 1::2] = -jac_tmp[1::2, ::2] 

650 

651 ml = getattr(self._integrator, 'ml', None) 

652 mu = getattr(self._integrator, 'mu', None) 

653 if ml is not None or mu is not None: 

654 # Jacobian is banded. The user's Jacobian function has computed 

655 # the complex Jacobian in packed format. The corresponding 

656 # real-valued version has every other column shifted up. 

657 jac_tmp = _transform_banded_jac(jac_tmp) 

658 

659 return jac_tmp 

660 

661 @property 

662 def y(self): 

663 return self._y[::2] + 1j * self._y[1::2] 

664 

665 def set_integrator(self, name, **integrator_params): 

666 """ 

667 Set integrator by name. 

668 

669 Parameters 

670 ---------- 

671 name : str 

672 Name of the integrator 

673 integrator_params 

674 Additional parameters for the integrator. 

675 """ 

676 if name == 'zvode': 

677 raise ValueError("zvode must be used with ode, not complex_ode") 

678 

679 lband = integrator_params.get('lband') 

680 uband = integrator_params.get('uband') 

681 if lband is not None or uband is not None: 

682 # The Jacobian is banded. Override the user-supplied bandwidths 

683 # (which are for the complex Jacobian) with the bandwidths of 

684 # the corresponding real-valued Jacobian wrapper of the complex 

685 # Jacobian. 

686 integrator_params['lband'] = 2 * (lband or 0) + 1 

687 integrator_params['uband'] = 2 * (uband or 0) + 1 

688 

689 return ode.set_integrator(self, name, **integrator_params) 

690 

691 def set_initial_value(self, y, t=0.0): 

692 """Set initial conditions y(t) = y.""" 

693 y = asarray(y) 

694 self.tmp = zeros(y.size * 2, 'float') 

695 self.tmp[::2] = real(y) 

696 self.tmp[1::2] = imag(y) 

697 return ode.set_initial_value(self, self.tmp, t) 

698 

699 def integrate(self, t, step=False, relax=False): 

700 """Find y=y(t), set y as an initial condition, and return y. 

701 

702 Parameters 

703 ---------- 

704 t : float 

705 The endpoint of the integration step. 

706 step : bool 

707 If True, and if the integrator supports the step method, 

708 then perform a single integration step and return. 

709 This parameter is provided in order to expose internals of 

710 the implementation, and should not be changed from its default 

711 value in most cases. 

712 relax : bool 

713 If True and if the integrator supports the run_relax method, 

714 then integrate until t_1 >= t and return. ``relax`` is not 

715 referenced if ``step=True``. 

716 This parameter is provided in order to expose internals of 

717 the implementation, and should not be changed from its default 

718 value in most cases. 

719 

720 Returns 

721 ------- 

722 y : float 

723 The integrated value at t 

724 """ 

725 y = ode.integrate(self, t, step, relax) 

726 return y[::2] + 1j * y[1::2] 

727 

728 def set_solout(self, solout): 

729 """ 

730 Set callable to be called at every successful integration step. 

731 

732 Parameters 

733 ---------- 

734 solout : callable 

735 ``solout(t, y)`` is called at each internal integrator step, 

736 t is a scalar providing the current independent position 

737 y is the current soloution ``y.shape == (n,)`` 

738 solout should return -1 to stop integration 

739 otherwise it should return None or 0 

740 

741 """ 

742 if self._integrator.supports_solout: 

743 self._integrator.set_solout(solout, complex=True) 

744 else: 

745 raise TypeError("selected integrator does not support solouta," 

746 + "choose another one") 

747 

748 

749# ------------------------------------------------------------------------------ 

750# ODE integrators 

751# ------------------------------------------------------------------------------ 

752 

753def find_integrator(name): 

754 for cl in IntegratorBase.integrator_classes: 

755 if re.match(name, cl.__name__, re.I): 

756 return cl 

757 return None 

758 

759 

760class IntegratorConcurrencyError(RuntimeError): 

761 """ 

762 Failure due to concurrent usage of an integrator that can be used 

763 only for a single problem at a time. 

764 

765 """ 

766 

767 def __init__(self, name): 

768 msg = ("Integrator `%s` can be used to solve only a single problem " 

769 "at a time. If you want to integrate multiple problems, " 

770 "consider using a different integrator " 

771 "(see `ode.set_integrator`)") % name 

772 RuntimeError.__init__(self, msg) 

773 

774 

775class IntegratorBase(object): 

776 runner = None # runner is None => integrator is not available 

777 success = None # success==1 if integrator was called successfully 

778 istate = None # istate > 0 means success, istate < 0 means failure 

779 supports_run_relax = None 

780 supports_step = None 

781 supports_solout = False 

782 integrator_classes = [] 

783 scalar = float 

784 

785 def acquire_new_handle(self): 

786 # Some of the integrators have internal state (ancient 

787 # Fortran...), and so only one instance can use them at a time. 

788 # We keep track of this, and fail when concurrent usage is tried. 

789 self.__class__.active_global_handle += 1 

790 self.handle = self.__class__.active_global_handle 

791 

792 def check_handle(self): 

793 if self.handle is not self.__class__.active_global_handle: 

794 raise IntegratorConcurrencyError(self.__class__.__name__) 

795 

796 def reset(self, n, has_jac): 

797 """Prepare integrator for call: allocate memory, set flags, etc. 

798 n - number of equations. 

799 has_jac - if user has supplied function for evaluating Jacobian. 

800 """ 

801 

802 def run(self, f, jac, y0, t0, t1, f_params, jac_params): 

803 """Integrate from t=t0 to t=t1 using y0 as an initial condition. 

804 Return 2-tuple (y1,t1) where y1 is the result and t=t1 

805 defines the stoppage coordinate of the result. 

806 """ 

807 raise NotImplementedError('all integrators must define ' 

808 'run(f, jac, t0, t1, y0, f_params, jac_params)') 

809 

810 def step(self, f, jac, y0, t0, t1, f_params, jac_params): 

811 """Make one integration step and return (y1,t1).""" 

812 raise NotImplementedError('%s does not support step() method' % 

813 self.__class__.__name__) 

814 

815 def run_relax(self, f, jac, y0, t0, t1, f_params, jac_params): 

816 """Integrate from t=t0 to t>=t1 and return (y1,t).""" 

817 raise NotImplementedError('%s does not support run_relax() method' % 

818 self.__class__.__name__) 

819 

820 # XXX: __str__ method for getting visual state of the integrator 

821 

822 

823def _vode_banded_jac_wrapper(jacfunc, ml, jac_params): 

824 """ 

825 Wrap a banded Jacobian function with a function that pads 

826 the Jacobian with `ml` rows of zeros. 

827 """ 

828 

829 def jac_wrapper(t, y): 

830 jac = asarray(jacfunc(t, y, *jac_params)) 

831 padded_jac = vstack((jac, zeros((ml, jac.shape[1])))) 

832 return padded_jac 

833 

834 return jac_wrapper 

835 

836 

837class vode(IntegratorBase): 

838 runner = getattr(_vode, 'dvode', None) 

839 

840 messages = {-1: 'Excess work done on this call. (Perhaps wrong MF.)', 

841 -2: 'Excess accuracy requested. (Tolerances too small.)', 

842 -3: 'Illegal input detected. (See printed message.)', 

843 -4: 'Repeated error test failures. (Check all input.)', 

844 -5: 'Repeated convergence failures. (Perhaps bad' 

845 ' Jacobian supplied or wrong choice of MF or tolerances.)', 

846 -6: 'Error weight became zero during problem. (Solution' 

847 ' component i vanished, and ATOL or ATOL(i) = 0.)' 

848 } 

849 supports_run_relax = 1 

850 supports_step = 1 

851 active_global_handle = 0 

852 

853 def __init__(self, 

854 method='adams', 

855 with_jacobian=False, 

856 rtol=1e-6, atol=1e-12, 

857 lband=None, uband=None, 

858 order=12, 

859 nsteps=500, 

860 max_step=0.0, # corresponds to infinite 

861 min_step=0.0, 

862 first_step=0.0, # determined by solver 

863 ): 

864 

865 if re.match(method, r'adams', re.I): 

866 self.meth = 1 

867 elif re.match(method, r'bdf', re.I): 

868 self.meth = 2 

869 else: 

870 raise ValueError('Unknown integration method %s' % method) 

871 self.with_jacobian = with_jacobian 

872 self.rtol = rtol 

873 self.atol = atol 

874 self.mu = uband 

875 self.ml = lband 

876 

877 self.order = order 

878 self.nsteps = nsteps 

879 self.max_step = max_step 

880 self.min_step = min_step 

881 self.first_step = first_step 

882 self.success = 1 

883 

884 self.initialized = False 

885 

886 def _determine_mf_and_set_bands(self, has_jac): 

887 """ 

888 Determine the `MF` parameter (Method Flag) for the Fortran subroutine `dvode`. 

889 

890 In the Fortran code, the legal values of `MF` are: 

891 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, 

892 -11, -12, -14, -15, -21, -22, -24, -25 

893 but this Python wrapper does not use negative values. 

894 

895 Returns 

896 

897 mf = 10*self.meth + miter 

898 

899 self.meth is the linear multistep method: 

900 self.meth == 1: method="adams" 

901 self.meth == 2: method="bdf" 

902 

903 miter is the correction iteration method: 

904 miter == 0: Functional iteraton; no Jacobian involved. 

905 miter == 1: Chord iteration with user-supplied full Jacobian. 

906 miter == 2: Chord iteration with internally computed full Jacobian. 

907 miter == 3: Chord iteration with internally computed diagonal Jacobian. 

908 miter == 4: Chord iteration with user-supplied banded Jacobian. 

909 miter == 5: Chord iteration with internally computed banded Jacobian. 

910 

911 Side effects: If either self.mu or self.ml is not None and the other is None, 

912 then the one that is None is set to 0. 

913 """ 

914 

915 jac_is_banded = self.mu is not None or self.ml is not None 

916 if jac_is_banded: 

917 if self.mu is None: 

918 self.mu = 0 

919 if self.ml is None: 

920 self.ml = 0 

921 

922 # has_jac is True if the user provided a Jacobian function. 

923 if has_jac: 

924 if jac_is_banded: 

925 miter = 4 

926 else: 

927 miter = 1 

928 else: 

929 if jac_is_banded: 

930 if self.ml == self.mu == 0: 

931 miter = 3 # Chord iteration with internal diagonal Jacobian. 

932 else: 

933 miter = 5 # Chord iteration with internal banded Jacobian. 

934 else: 

935 # self.with_jacobian is set by the user in the call to ode.set_integrator. 

936 if self.with_jacobian: 

937 miter = 2 # Chord iteration with internal full Jacobian. 

938 else: 

939 miter = 0 # Functional iteraton; no Jacobian involved. 

940 

941 mf = 10 * self.meth + miter 

942 return mf 

943 

944 def reset(self, n, has_jac): 

945 mf = self._determine_mf_and_set_bands(has_jac) 

946 

947 if mf == 10: 

948 lrw = 20 + 16 * n 

949 elif mf in [11, 12]: 

950 lrw = 22 + 16 * n + 2 * n * n 

951 elif mf == 13: 

952 lrw = 22 + 17 * n 

953 elif mf in [14, 15]: 

954 lrw = 22 + 18 * n + (3 * self.ml + 2 * self.mu) * n 

955 elif mf == 20: 

956 lrw = 20 + 9 * n 

957 elif mf in [21, 22]: 

958 lrw = 22 + 9 * n + 2 * n * n 

959 elif mf == 23: 

960 lrw = 22 + 10 * n 

961 elif mf in [24, 25]: 

962 lrw = 22 + 11 * n + (3 * self.ml + 2 * self.mu) * n 

963 else: 

964 raise ValueError('Unexpected mf=%s' % mf) 

965 

966 if mf % 10 in [0, 3]: 

967 liw = 30 

968 else: 

969 liw = 30 + n 

970 

971 rwork = zeros((lrw,), float) 

972 rwork[4] = self.first_step 

973 rwork[5] = self.max_step 

974 rwork[6] = self.min_step 

975 self.rwork = rwork 

976 

977 iwork = zeros((liw,), _vode_int_dtype) 

978 if self.ml is not None: 

979 iwork[0] = self.ml 

980 if self.mu is not None: 

981 iwork[1] = self.mu 

982 iwork[4] = self.order 

983 iwork[5] = self.nsteps 

984 iwork[6] = 2 # mxhnil 

985 self.iwork = iwork 

986 

987 self.call_args = [self.rtol, self.atol, 1, 1, 

988 self.rwork, self.iwork, mf] 

989 self.success = 1 

990 self.initialized = False 

991 

992 def run(self, f, jac, y0, t0, t1, f_params, jac_params): 

993 if self.initialized: 

994 self.check_handle() 

995 else: 

996 self.initialized = True 

997 self.acquire_new_handle() 

998 

999 if self.ml is not None and self.ml > 0: 

1000 # Banded Jacobian. Wrap the user-provided function with one 

1001 # that pads the Jacobian array with the extra `self.ml` rows 

1002 # required by the f2py-generated wrapper. 

1003 jac = _vode_banded_jac_wrapper(jac, self.ml, jac_params) 

1004 

1005 args = ((f, jac, y0, t0, t1) + tuple(self.call_args) + 

1006 (f_params, jac_params)) 

1007 y1, t, istate = self.runner(*args) 

1008 self.istate = istate 

1009 if istate < 0: 

1010 unexpected_istate_msg = 'Unexpected istate={:d}'.format(istate) 

1011 warnings.warn('{:s}: {:s}'.format(self.__class__.__name__, 

1012 self.messages.get(istate, unexpected_istate_msg))) 

1013 self.success = 0 

1014 else: 

1015 self.call_args[3] = 2 # upgrade istate from 1 to 2 

1016 self.istate = 2 

1017 return y1, t 

1018 

1019 def step(self, *args): 

1020 itask = self.call_args[2] 

1021 self.call_args[2] = 2 

1022 r = self.run(*args) 

1023 self.call_args[2] = itask 

1024 return r 

1025 

1026 def run_relax(self, *args): 

1027 itask = self.call_args[2] 

1028 self.call_args[2] = 3 

1029 r = self.run(*args) 

1030 self.call_args[2] = itask 

1031 return r 

1032 

1033 

1034if vode.runner is not None: 

1035 IntegratorBase.integrator_classes.append(vode) 

1036 

1037 

1038class zvode(vode): 

1039 runner = getattr(_vode, 'zvode', None) 

1040 

1041 supports_run_relax = 1 

1042 supports_step = 1 

1043 scalar = complex 

1044 active_global_handle = 0 

1045 

1046 def reset(self, n, has_jac): 

1047 mf = self._determine_mf_and_set_bands(has_jac) 

1048 

1049 if mf in (10,): 

1050 lzw = 15 * n 

1051 elif mf in (11, 12): 

1052 lzw = 15 * n + 2 * n ** 2 

1053 elif mf in (-11, -12): 

1054 lzw = 15 * n + n ** 2 

1055 elif mf in (13,): 

1056 lzw = 16 * n 

1057 elif mf in (14, 15): 

1058 lzw = 17 * n + (3 * self.ml + 2 * self.mu) * n 

1059 elif mf in (-14, -15): 

1060 lzw = 16 * n + (2 * self.ml + self.mu) * n 

1061 elif mf in (20,): 

1062 lzw = 8 * n 

1063 elif mf in (21, 22): 

1064 lzw = 8 * n + 2 * n ** 2 

1065 elif mf in (-21, -22): 

1066 lzw = 8 * n + n ** 2 

1067 elif mf in (23,): 

1068 lzw = 9 * n 

1069 elif mf in (24, 25): 

1070 lzw = 10 * n + (3 * self.ml + 2 * self.mu) * n 

1071 elif mf in (-24, -25): 

1072 lzw = 9 * n + (2 * self.ml + self.mu) * n 

1073 

1074 lrw = 20 + n 

1075 

1076 if mf % 10 in (0, 3): 

1077 liw = 30 

1078 else: 

1079 liw = 30 + n 

1080 

1081 zwork = zeros((lzw,), complex) 

1082 self.zwork = zwork 

1083 

1084 rwork = zeros((lrw,), float) 

1085 rwork[4] = self.first_step 

1086 rwork[5] = self.max_step 

1087 rwork[6] = self.min_step 

1088 self.rwork = rwork 

1089 

1090 iwork = zeros((liw,), _vode_int_dtype) 

1091 if self.ml is not None: 

1092 iwork[0] = self.ml 

1093 if self.mu is not None: 

1094 iwork[1] = self.mu 

1095 iwork[4] = self.order 

1096 iwork[5] = self.nsteps 

1097 iwork[6] = 2 # mxhnil 

1098 self.iwork = iwork 

1099 

1100 self.call_args = [self.rtol, self.atol, 1, 1, 

1101 self.zwork, self.rwork, self.iwork, mf] 

1102 self.success = 1 

1103 self.initialized = False 

1104 

1105 

1106if zvode.runner is not None: 

1107 IntegratorBase.integrator_classes.append(zvode) 

1108 

1109 

1110class dopri5(IntegratorBase): 

1111 runner = getattr(_dop, 'dopri5', None) 

1112 name = 'dopri5' 

1113 supports_solout = True 

1114 

1115 messages = {1: 'computation successful', 

1116 2: 'computation successful (interrupted by solout)', 

1117 -1: 'input is not consistent', 

1118 -2: 'larger nsteps is needed', 

1119 -3: 'step size becomes too small', 

1120 -4: 'problem is probably stiff (interrupted)', 

1121 } 

1122 

1123 def __init__(self, 

1124 rtol=1e-6, atol=1e-12, 

1125 nsteps=500, 

1126 max_step=0.0, 

1127 first_step=0.0, # determined by solver 

1128 safety=0.9, 

1129 ifactor=10.0, 

1130 dfactor=0.2, 

1131 beta=0.0, 

1132 method=None, 

1133 verbosity=-1, # no messages if negative 

1134 ): 

1135 self.rtol = rtol 

1136 self.atol = atol 

1137 self.nsteps = nsteps 

1138 self.max_step = max_step 

1139 self.first_step = first_step 

1140 self.safety = safety 

1141 self.ifactor = ifactor 

1142 self.dfactor = dfactor 

1143 self.beta = beta 

1144 self.verbosity = verbosity 

1145 self.success = 1 

1146 self.set_solout(None) 

1147 

1148 def set_solout(self, solout, complex=False): 

1149 self.solout = solout 

1150 self.solout_cmplx = complex 

1151 if solout is None: 

1152 self.iout = 0 

1153 else: 

1154 self.iout = 1 

1155 

1156 def reset(self, n, has_jac): 

1157 work = zeros((8 * n + 21,), float) 

1158 work[1] = self.safety 

1159 work[2] = self.dfactor 

1160 work[3] = self.ifactor 

1161 work[4] = self.beta 

1162 work[5] = self.max_step 

1163 work[6] = self.first_step 

1164 self.work = work 

1165 iwork = zeros((21,), _dop_int_dtype) 

1166 iwork[0] = self.nsteps 

1167 iwork[2] = self.verbosity 

1168 self.iwork = iwork 

1169 self.call_args = [self.rtol, self.atol, self._solout, 

1170 self.iout, self.work, self.iwork] 

1171 self.success = 1 

1172 

1173 def run(self, f, jac, y0, t0, t1, f_params, jac_params): 

1174 x, y, iwork, istate = self.runner(*((f, t0, y0, t1) + 

1175 tuple(self.call_args) + (f_params,))) 

1176 self.istate = istate 

1177 if istate < 0: 

1178 unexpected_istate_msg = 'Unexpected istate={:d}'.format(istate) 

1179 warnings.warn('{:s}: {:s}'.format(self.__class__.__name__, 

1180 self.messages.get(istate, unexpected_istate_msg))) 

1181 self.success = 0 

1182 return y, x 

1183 

1184 def _solout(self, nr, xold, x, y, nd, icomp, con): 

1185 if self.solout is not None: 

1186 if self.solout_cmplx: 

1187 y = y[::2] + 1j * y[1::2] 

1188 return self.solout(x, y) 

1189 else: 

1190 return 1 

1191 

1192 

1193if dopri5.runner is not None: 

1194 IntegratorBase.integrator_classes.append(dopri5) 

1195 

1196 

1197class dop853(dopri5): 

1198 runner = getattr(_dop, 'dop853', None) 

1199 name = 'dop853' 

1200 

1201 def __init__(self, 

1202 rtol=1e-6, atol=1e-12, 

1203 nsteps=500, 

1204 max_step=0.0, 

1205 first_step=0.0, # determined by solver 

1206 safety=0.9, 

1207 ifactor=6.0, 

1208 dfactor=0.3, 

1209 beta=0.0, 

1210 method=None, 

1211 verbosity=-1, # no messages if negative 

1212 ): 

1213 super(self.__class__, self).__init__(rtol, atol, nsteps, max_step, 

1214 first_step, safety, ifactor, 

1215 dfactor, beta, method, 

1216 verbosity) 

1217 

1218 def reset(self, n, has_jac): 

1219 work = zeros((11 * n + 21,), float) 

1220 work[1] = self.safety 

1221 work[2] = self.dfactor 

1222 work[3] = self.ifactor 

1223 work[4] = self.beta 

1224 work[5] = self.max_step 

1225 work[6] = self.first_step 

1226 self.work = work 

1227 iwork = zeros((21,), _dop_int_dtype) 

1228 iwork[0] = self.nsteps 

1229 iwork[2] = self.verbosity 

1230 self.iwork = iwork 

1231 self.call_args = [self.rtol, self.atol, self._solout, 

1232 self.iout, self.work, self.iwork] 

1233 self.success = 1 

1234 

1235 

1236if dop853.runner is not None: 

1237 IntegratorBase.integrator_classes.append(dop853) 

1238 

1239 

1240class lsoda(IntegratorBase): 

1241 runner = getattr(_lsoda, 'lsoda', None) 

1242 active_global_handle = 0 

1243 

1244 messages = { 

1245 2: "Integration successful.", 

1246 -1: "Excess work done on this call (perhaps wrong Dfun type).", 

1247 -2: "Excess accuracy requested (tolerances too small).", 

1248 -3: "Illegal input detected (internal error).", 

1249 -4: "Repeated error test failures (internal error).", 

1250 -5: "Repeated convergence failures (perhaps bad Jacobian or tolerances).", 

1251 -6: "Error weight became zero during problem.", 

1252 -7: "Internal workspace insufficient to finish (internal error)." 

1253 } 

1254 

1255 def __init__(self, 

1256 with_jacobian=False, 

1257 rtol=1e-6, atol=1e-12, 

1258 lband=None, uband=None, 

1259 nsteps=500, 

1260 max_step=0.0, # corresponds to infinite 

1261 min_step=0.0, 

1262 first_step=0.0, # determined by solver 

1263 ixpr=0, 

1264 max_hnil=0, 

1265 max_order_ns=12, 

1266 max_order_s=5, 

1267 method=None 

1268 ): 

1269 

1270 self.with_jacobian = with_jacobian 

1271 self.rtol = rtol 

1272 self.atol = atol 

1273 self.mu = uband 

1274 self.ml = lband 

1275 

1276 self.max_order_ns = max_order_ns 

1277 self.max_order_s = max_order_s 

1278 self.nsteps = nsteps 

1279 self.max_step = max_step 

1280 self.min_step = min_step 

1281 self.first_step = first_step 

1282 self.ixpr = ixpr 

1283 self.max_hnil = max_hnil 

1284 self.success = 1 

1285 

1286 self.initialized = False 

1287 

1288 def reset(self, n, has_jac): 

1289 # Calculate parameters for Fortran subroutine dvode. 

1290 if has_jac: 

1291 if self.mu is None and self.ml is None: 

1292 jt = 1 

1293 else: 

1294 if self.mu is None: 

1295 self.mu = 0 

1296 if self.ml is None: 

1297 self.ml = 0 

1298 jt = 4 

1299 else: 

1300 if self.mu is None and self.ml is None: 

1301 jt = 2 

1302 else: 

1303 if self.mu is None: 

1304 self.mu = 0 

1305 if self.ml is None: 

1306 self.ml = 0 

1307 jt = 5 

1308 lrn = 20 + (self.max_order_ns + 4) * n 

1309 if jt in [1, 2]: 

1310 lrs = 22 + (self.max_order_s + 4) * n + n * n 

1311 elif jt in [4, 5]: 

1312 lrs = 22 + (self.max_order_s + 5 + 2 * self.ml + self.mu) * n 

1313 else: 

1314 raise ValueError('Unexpected jt=%s' % jt) 

1315 lrw = max(lrn, lrs) 

1316 liw = 20 + n 

1317 rwork = zeros((lrw,), float) 

1318 rwork[4] = self.first_step 

1319 rwork[5] = self.max_step 

1320 rwork[6] = self.min_step 

1321 self.rwork = rwork 

1322 iwork = zeros((liw,), _lsoda_int_dtype) 

1323 if self.ml is not None: 

1324 iwork[0] = self.ml 

1325 if self.mu is not None: 

1326 iwork[1] = self.mu 

1327 iwork[4] = self.ixpr 

1328 iwork[5] = self.nsteps 

1329 iwork[6] = self.max_hnil 

1330 iwork[7] = self.max_order_ns 

1331 iwork[8] = self.max_order_s 

1332 self.iwork = iwork 

1333 self.call_args = [self.rtol, self.atol, 1, 1, 

1334 self.rwork, self.iwork, jt] 

1335 self.success = 1 

1336 self.initialized = False 

1337 

1338 def run(self, f, jac, y0, t0, t1, f_params, jac_params): 

1339 if self.initialized: 

1340 self.check_handle() 

1341 else: 

1342 self.initialized = True 

1343 self.acquire_new_handle() 

1344 args = [f, y0, t0, t1] + self.call_args[:-1] + \ 

1345 [jac, self.call_args[-1], f_params, 0, jac_params] 

1346 y1, t, istate = self.runner(*args) 

1347 self.istate = istate 

1348 if istate < 0: 

1349 unexpected_istate_msg = 'Unexpected istate={:d}'.format(istate) 

1350 warnings.warn('{:s}: {:s}'.format(self.__class__.__name__, 

1351 self.messages.get(istate, unexpected_istate_msg))) 

1352 self.success = 0 

1353 else: 

1354 self.call_args[3] = 2 # upgrade istate from 1 to 2 

1355 self.istate = 2 

1356 return y1, t 

1357 

1358 def step(self, *args): 

1359 itask = self.call_args[2] 

1360 self.call_args[2] = 2 

1361 r = self.run(*args) 

1362 self.call_args[2] = itask 

1363 return r 

1364 

1365 def run_relax(self, *args): 

1366 itask = self.call_args[2] 

1367 self.call_args[2] = 3 

1368 r = self.run(*args) 

1369 self.call_args[2] = itask 

1370 return r 

1371 

1372 

1373if lsoda.runner: 

1374 IntegratorBase.integrator_classes.append(lsoda)