Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/integrate/_ode.py : 21%

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.
5User-friendly interface to various numerical integrators for solving a
6system of first order ODEs with prescribed initial conditions::
8 d y(t)[i]
9 --------- = f(t,y(t))[i],
10 d t
12 y(t=0)[i] = y0[i],
14where::
16 i = 0, ..., len(y0) - 1
18class ode
19---------
21A generic interface class to numeric integrators. It has the following
22methods::
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()
32class complex_ode
33-----------------
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)
81__all__ = ['ode', 'complex_ode']
82__version__ = "$Id$"
83__docformat__ = "restructuredtext en"
85import re
86import warnings
88from numpy import asarray, array, zeros, isscalar, real, imag, vstack
90from . import vode as _vode
91from . import _dop
92from . import lsoda as _lsoda
95_dop_int_dtype = _dop.types.intvar.dtype
96_vode_int_dtype = _vode.types.intvar.dtype
97_lsoda_int_dtype = _lsoda.types.intvar.dtype
100# ------------------------------------------------------------------------------
101# User interface
102# ------------------------------------------------------------------------------
105class ode(object):
106 """
107 A generic interface class to numeric integrators.
109 Solve an equation system :math:`y'(t) = f(t,y)` with (optional) ``jac = df/dy``.
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`.
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)``.
126 Attributes
127 ----------
128 t : float
129 Current time.
130 y : ndarray
131 Current variable values.
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
138 Notes
139 -----
140 Available integrators are listed below. They can be selected using
141 the `set_integrator` method.
143 "vode"
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).
150 Source: http://www.netlib.org/ode/vode.f
152 .. warning::
154 This integrator is not re-entrant. You cannot have two `ode`
155 instances using the "vode" integrator at the same time.
157 This integrator accepts the following parameters in `set_integrator`
158 method of the `ode` class:
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.
190 "zvode"
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).
197 Source: http://www.netlib.org/ode/zvode.f
199 .. warning::
201 This integrator is not re-entrant. You cannot have two `ode`
202 instances using the "zvode" integrator at the same time.
204 This integrator accepts the same parameters in `set_integrator`
205 as the "vode" solver.
207 .. note::
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).
219 "lsoda"
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).
227 Source: http://www.netlib.org/odepack
229 .. warning::
231 This integrator is not re-entrant. You cannot have two `ode`
232 instances using the "lsoda" integrator at the same time.
234 This integrator accepts the following parameters in `set_integrator`
235 method of the `ode` class:
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).
265 "dopri5"
267 This is an explicit runge-kutta method of order (4)5 due to Dormand &
268 Prince (with stepsize control and dense output).
270 Authors:
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
277 This code is described in [HNW93]_.
279 This integrator accepts the following parameters in set_integrator()
280 method of the ode class:
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).
301 "dop853"
303 This is an explicit runge-kutta method of order 8(5,3) due to Dormand
304 & Prince (with stepsize control and dense output).
306 Options and references the same as "dopri5".
308 Examples
309 --------
311 A problem to integrate and the corresponding jacobian:
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]]]
322 The integration:
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 ]
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)
348 """
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 = []
358 @property
359 def y(self):
360 return self._y
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
374 def set_integrator(self, name, **integrator_params):
375 """
376 Set integrator by name.
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
399 def integrate(self, t, step=False, relax=False):
400 """Find y=y(t), set y as an initial condition, and return y.
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.
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
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.')
440 return self._y
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
450 def get_return_code(self):
451 """Extracts the return code for the integration to enable better control
452 if the integration fails.
454 In general, a return code > 0 implies success, while a return code < 0
455 implies failure.
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.
462 "vode"
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 =========== =======
478 "zvode"
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 =========== =======
494 "dopri5"
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 =========== =======
507 "dop853"
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 =========== =======
520 "lsoda"
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
541 def set_f_params(self, *args):
542 """Set extra parameters for user-supplied function f."""
543 self.f_params = args
544 return self
546 def set_jac_params(self, *args):
547 """Set extra parameters for user-supplied function jac."""
548 self.jac_params = args
549 return self
551 def set_solout(self, solout):
552 """
553 Set callable to be called at every successful integration step.
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
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")
574def _transform_banded_jac(bjac):
575 """
576 Convert a real matrix of the form (for example)
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]
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
593class complex_ode(ode):
594 """
595 A wrapper of ode for complex systems.
597 This functions similarly as `ode`, but re-maps a complex-valued
598 equation system to a real-valued one before using the integrators.
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)``.
609 Attributes
610 ----------
611 t : float
612 Current time.
613 y : ndarray
614 Current variable values.
616 Examples
617 --------
618 For usage examples, see `ode`.
620 """
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)
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
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))
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]
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)
659 return jac_tmp
661 @property
662 def y(self):
663 return self._y[::2] + 1j * self._y[1::2]
665 def set_integrator(self, name, **integrator_params):
666 """
667 Set integrator by name.
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")
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
689 return ode.set_integrator(self, name, **integrator_params)
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)
699 def integrate(self, t, step=False, relax=False):
700 """Find y=y(t), set y as an initial condition, and return y.
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.
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]
728 def set_solout(self, solout):
729 """
730 Set callable to be called at every successful integration step.
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
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")
749# ------------------------------------------------------------------------------
750# ODE integrators
751# ------------------------------------------------------------------------------
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
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.
765 """
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)
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
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
792 def check_handle(self):
793 if self.handle is not self.__class__.active_global_handle:
794 raise IntegratorConcurrencyError(self.__class__.__name__)
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 """
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)')
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__)
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__)
820 # XXX: __str__ method for getting visual state of the integrator
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 """
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
834 return jac_wrapper
837class vode(IntegratorBase):
838 runner = getattr(_vode, 'dvode', None)
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
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 ):
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
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
884 self.initialized = False
886 def _determine_mf_and_set_bands(self, has_jac):
887 """
888 Determine the `MF` parameter (Method Flag) for the Fortran subroutine `dvode`.
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.
895 Returns
897 mf = 10*self.meth + miter
899 self.meth is the linear multistep method:
900 self.meth == 1: method="adams"
901 self.meth == 2: method="bdf"
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.
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 """
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
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.
941 mf = 10 * self.meth + miter
942 return mf
944 def reset(self, n, has_jac):
945 mf = self._determine_mf_and_set_bands(has_jac)
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)
966 if mf % 10 in [0, 3]:
967 liw = 30
968 else:
969 liw = 30 + n
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
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
987 self.call_args = [self.rtol, self.atol, 1, 1,
988 self.rwork, self.iwork, mf]
989 self.success = 1
990 self.initialized = False
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()
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)
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
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
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
1034if vode.runner is not None:
1035 IntegratorBase.integrator_classes.append(vode)
1038class zvode(vode):
1039 runner = getattr(_vode, 'zvode', None)
1041 supports_run_relax = 1
1042 supports_step = 1
1043 scalar = complex
1044 active_global_handle = 0
1046 def reset(self, n, has_jac):
1047 mf = self._determine_mf_and_set_bands(has_jac)
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
1074 lrw = 20 + n
1076 if mf % 10 in (0, 3):
1077 liw = 30
1078 else:
1079 liw = 30 + n
1081 zwork = zeros((lzw,), complex)
1082 self.zwork = zwork
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
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
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
1106if zvode.runner is not None:
1107 IntegratorBase.integrator_classes.append(zvode)
1110class dopri5(IntegratorBase):
1111 runner = getattr(_dop, 'dopri5', None)
1112 name = 'dopri5'
1113 supports_solout = True
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 }
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)
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
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
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
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
1193if dopri5.runner is not None:
1194 IntegratorBase.integrator_classes.append(dopri5)
1197class dop853(dopri5):
1198 runner = getattr(_dop, 'dop853', None)
1199 name = 'dop853'
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)
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
1236if dop853.runner is not None:
1237 IntegratorBase.integrator_classes.append(dop853)
1240class lsoda(IntegratorBase):
1241 runner = getattr(_lsoda, 'lsoda', None)
1242 active_global_handle = 0
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 }
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 ):
1270 self.with_jacobian = with_jacobian
1271 self.rtol = rtol
1272 self.atol = atol
1273 self.mu = uband
1274 self.ml = lband
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
1286 self.initialized = False
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
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
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
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
1373if lsoda.runner:
1374 IntegratorBase.integrator_classes.append(lsoda)