Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_ode.py: 21%
507 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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']
83import re
84import warnings
86from numpy import asarray, array, zeros, isscalar, real, imag, vstack
88from . import _vode
89from . import _dop
90from . import _lsoda
93_dop_int_dtype = _dop.types.intvar.dtype
94_vode_int_dtype = _vode.types.intvar.dtype
95_lsoda_int_dtype = _lsoda.types.intvar.dtype
98# ------------------------------------------------------------------------------
99# User interface
100# ------------------------------------------------------------------------------
103class ode:
104 """
105 A generic interface class to numeric integrators.
107 Solve an equation system :math:`y'(t) = f(t,y)` with (optional) ``jac = df/dy``.
109 *Note*: The first two arguments of ``f(t, y, ...)`` are in the
110 opposite order of the arguments in the system definition function used
111 by `scipy.integrate.odeint`.
113 Parameters
114 ----------
115 f : callable ``f(t, y, *f_args)``
116 Right-hand side of the differential equation. t is a scalar,
117 ``y.shape == (n,)``.
118 ``f_args`` is set by calling ``set_f_params(*args)``.
119 `f` should return a scalar, array or list (not a tuple).
120 jac : callable ``jac(t, y, *jac_args)``, optional
121 Jacobian of the right-hand side, ``jac[i,j] = d f[i] / d y[j]``.
122 ``jac_args`` is set by calling ``set_jac_params(*args)``.
124 Attributes
125 ----------
126 t : float
127 Current time.
128 y : ndarray
129 Current variable values.
131 See also
132 --------
133 odeint : an integrator with a simpler interface based on lsoda from ODEPACK
134 quad : for finding the area under a curve
136 Notes
137 -----
138 Available integrators are listed below. They can be selected using
139 the `set_integrator` method.
141 "vode"
143 Real-valued Variable-coefficient Ordinary Differential Equation
144 solver, with fixed-leading-coefficient implementation. It provides
145 implicit Adams method (for non-stiff problems) and a method based on
146 backward differentiation formulas (BDF) (for stiff problems).
148 Source: http://www.netlib.org/ode/vode.f
150 .. warning::
152 This integrator is not re-entrant. You cannot have two `ode`
153 instances using the "vode" integrator at the same time.
155 This integrator accepts the following parameters in `set_integrator`
156 method of the `ode` class:
158 - atol : float or sequence
159 absolute tolerance for solution
160 - rtol : float or sequence
161 relative tolerance for solution
162 - lband : None or int
163 - uband : None or int
164 Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
165 Setting these requires your jac routine to return the jacobian
166 in packed format, jac_packed[i-j+uband, j] = jac[i,j]. The
167 dimension of the matrix must be (lband+uband+1, len(y)).
168 - method: 'adams' or 'bdf'
169 Which solver to use, Adams (non-stiff) or BDF (stiff)
170 - with_jacobian : bool
171 This option is only considered when the user has not supplied a
172 Jacobian function and has not indicated (by setting either band)
173 that the Jacobian is banded. In this case, `with_jacobian` specifies
174 whether the iteration method of the ODE solver's correction step is
175 chord iteration with an internally generated full Jacobian or
176 functional iteration with no Jacobian.
177 - nsteps : int
178 Maximum number of (internally defined) steps allowed during one
179 call to the solver.
180 - first_step : float
181 - min_step : float
182 - max_step : float
183 Limits for the step sizes used by the integrator.
184 - order : int
185 Maximum order used by the integrator,
186 order <= 12 for Adams, <= 5 for BDF.
188 "zvode"
190 Complex-valued Variable-coefficient Ordinary Differential Equation
191 solver, with fixed-leading-coefficient implementation. It provides
192 implicit Adams method (for non-stiff problems) and a method based on
193 backward differentiation formulas (BDF) (for stiff problems).
195 Source: http://www.netlib.org/ode/zvode.f
197 .. warning::
199 This integrator is not re-entrant. You cannot have two `ode`
200 instances using the "zvode" integrator at the same time.
202 This integrator accepts the same parameters in `set_integrator`
203 as the "vode" solver.
205 .. note::
207 When using ZVODE for a stiff system, it should only be used for
208 the case in which the function f is analytic, that is, when each f(i)
209 is an analytic function of each y(j). Analyticity means that the
210 partial derivative df(i)/dy(j) is a unique complex number, and this
211 fact is critical in the way ZVODE solves the dense or banded linear
212 systems that arise in the stiff case. For a complex stiff ODE system
213 in which f is not analytic, ZVODE is likely to have convergence
214 failures, and for this problem one should instead use DVODE on the
215 equivalent real system (in the real and imaginary parts of y).
217 "lsoda"
219 Real-valued Variable-coefficient Ordinary Differential Equation
220 solver, with fixed-leading-coefficient implementation. It provides
221 automatic method switching between implicit Adams method (for non-stiff
222 problems) and a method based on backward differentiation formulas (BDF)
223 (for stiff problems).
225 Source: http://www.netlib.org/odepack
227 .. warning::
229 This integrator is not re-entrant. You cannot have two `ode`
230 instances using the "lsoda" integrator at the same time.
232 This integrator accepts the following parameters in `set_integrator`
233 method of the `ode` class:
235 - atol : float or sequence
236 absolute tolerance for solution
237 - rtol : float or sequence
238 relative tolerance for solution
239 - lband : None or int
240 - uband : None or int
241 Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+uband.
242 Setting these requires your jac routine to return the jacobian
243 in packed format, jac_packed[i-j+uband, j] = jac[i,j].
244 - with_jacobian : bool
245 *Not used.*
246 - nsteps : int
247 Maximum number of (internally defined) steps allowed during one
248 call to the solver.
249 - first_step : float
250 - min_step : float
251 - max_step : float
252 Limits for the step sizes used by the integrator.
253 - max_order_ns : int
254 Maximum order used in the nonstiff case (default 12).
255 - max_order_s : int
256 Maximum order used in the stiff case (default 5).
257 - max_hnil : int
258 Maximum number of messages reporting too small step size (t + h = t)
259 (default 0)
260 - ixpr : int
261 Whether to generate extra printing at method switches (default False).
263 "dopri5"
265 This is an explicit runge-kutta method of order (4)5 due to Dormand &
266 Prince (with stepsize control and dense output).
268 Authors:
270 E. Hairer and G. Wanner
271 Universite de Geneve, Dept. de Mathematiques
272 CH-1211 Geneve 24, Switzerland
273 e-mail: ernst.hairer@math.unige.ch, gerhard.wanner@math.unige.ch
275 This code is described in [HNW93]_.
277 This integrator accepts the following parameters in set_integrator()
278 method of the ode class:
280 - atol : float or sequence
281 absolute tolerance for solution
282 - rtol : float or sequence
283 relative tolerance for solution
284 - nsteps : int
285 Maximum number of (internally defined) steps allowed during one
286 call to the solver.
287 - first_step : float
288 - max_step : float
289 - safety : float
290 Safety factor on new step selection (default 0.9)
291 - ifactor : float
292 - dfactor : float
293 Maximum factor to increase/decrease step size by in one step
294 - beta : float
295 Beta parameter for stabilised step size control.
296 - verbosity : int
297 Switch for printing messages (< 0 for no messages).
299 "dop853"
301 This is an explicit runge-kutta method of order 8(5,3) due to Dormand
302 & Prince (with stepsize control and dense output).
304 Options and references the same as "dopri5".
306 Examples
307 --------
309 A problem to integrate and the corresponding jacobian:
311 >>> from scipy.integrate import ode
312 >>>
313 >>> y0, t0 = [1.0j, 2.0], 0
314 >>>
315 >>> def f(t, y, arg1):
316 ... return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
317 >>> def jac(t, y, arg1):
318 ... return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
320 The integration:
322 >>> r = ode(f, jac).set_integrator('zvode', method='bdf')
323 >>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
324 >>> t1 = 10
325 >>> dt = 1
326 >>> while r.successful() and r.t < t1:
327 ... print(r.t+dt, r.integrate(r.t+dt))
328 1 [-0.71038232+0.23749653j 0.40000271+0.j ]
329 2.0 [0.19098503-0.52359246j 0.22222356+0.j ]
330 3.0 [0.47153208+0.52701229j 0.15384681+0.j ]
331 4.0 [-0.61905937+0.30726255j 0.11764744+0.j ]
332 5.0 [0.02340997-0.61418799j 0.09523835+0.j ]
333 6.0 [0.58643071+0.339819j 0.08000018+0.j ]
334 7.0 [-0.52070105+0.44525141j 0.06896565+0.j ]
335 8.0 [-0.15986733-0.61234476j 0.06060616+0.j ]
336 9.0 [0.64850462+0.15048982j 0.05405414+0.j ]
337 10.0 [-0.38404699+0.56382299j 0.04878055+0.j ]
339 References
340 ----------
341 .. [HNW93] E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary
342 Differential Equations i. Nonstiff Problems. 2nd edition.
343 Springer Series in Computational Mathematics,
344 Springer-Verlag (1993)
346 """
348 def __init__(self, f, jac=None):
349 self.stiff = 0
350 self.f = f
351 self.jac = jac
352 self.f_params = ()
353 self.jac_params = ()
354 self._y = []
356 @property
357 def y(self):
358 return self._y
360 def set_initial_value(self, y, t=0.0):
361 """Set initial conditions y(t) = y."""
362 if isscalar(y):
363 y = [y]
364 n_prev = len(self._y)
365 if not n_prev:
366 self.set_integrator('') # find first available integrator
367 self._y = asarray(y, self._integrator.scalar)
368 self.t = t
369 self._integrator.reset(len(self._y), self.jac is not None)
370 return self
372 def set_integrator(self, name, **integrator_params):
373 """
374 Set integrator by name.
376 Parameters
377 ----------
378 name : str
379 Name of the integrator.
380 **integrator_params
381 Additional parameters for the integrator.
382 """
383 integrator = find_integrator(name)
384 if integrator is None:
385 # FIXME: this really should be raise an exception. Will that break
386 # any code?
387 warnings.warn('No integrator name match with %r or is not '
388 'available.' % name)
389 else:
390 self._integrator = integrator(**integrator_params)
391 if not len(self._y):
392 self.t = 0.0
393 self._y = array([0.0], self._integrator.scalar)
394 self._integrator.reset(len(self._y), self.jac is not None)
395 return self
397 def integrate(self, t, step=False, relax=False):
398 """Find y=y(t), set y as an initial condition, and return y.
400 Parameters
401 ----------
402 t : float
403 The endpoint of the integration step.
404 step : bool
405 If True, and if the integrator supports the step method,
406 then perform a single integration step and return.
407 This parameter is provided in order to expose internals of
408 the implementation, and should not be changed from its default
409 value in most cases.
410 relax : bool
411 If True and if the integrator supports the run_relax method,
412 then integrate until t_1 >= t and return. ``relax`` is not
413 referenced if ``step=True``.
414 This parameter is provided in order to expose internals of
415 the implementation, and should not be changed from its default
416 value in most cases.
418 Returns
419 -------
420 y : float
421 The integrated value at t
422 """
423 if step and self._integrator.supports_step:
424 mth = self._integrator.step
425 elif relax and self._integrator.supports_run_relax:
426 mth = self._integrator.run_relax
427 else:
428 mth = self._integrator.run
430 try:
431 self._y, self.t = mth(self.f, self.jac or (lambda: None),
432 self._y, self.t, t,
433 self.f_params, self.jac_params)
434 except SystemError as e:
435 # f2py issue with tuple returns, see ticket 1187.
436 raise ValueError(
437 'Function to integrate must not return a tuple.'
438 ) from e
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:
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().__init__(rtol, atol, nsteps, max_step, first_step, safety,
1214 ifactor, dfactor, beta, method, verbosity)
1216 def reset(self, n, has_jac):
1217 work = zeros((11 * n + 21,), float)
1218 work[1] = self.safety
1219 work[2] = self.dfactor
1220 work[3] = self.ifactor
1221 work[4] = self.beta
1222 work[5] = self.max_step
1223 work[6] = self.first_step
1224 self.work = work
1225 iwork = zeros((21,), _dop_int_dtype)
1226 iwork[0] = self.nsteps
1227 iwork[2] = self.verbosity
1228 self.iwork = iwork
1229 self.call_args = [self.rtol, self.atol, self._solout,
1230 self.iout, self.work, self.iwork]
1231 self.success = 1
1234if dop853.runner is not None:
1235 IntegratorBase.integrator_classes.append(dop853)
1238class lsoda(IntegratorBase):
1239 runner = getattr(_lsoda, 'lsoda', None)
1240 active_global_handle = 0
1242 messages = {
1243 2: "Integration successful.",
1244 -1: "Excess work done on this call (perhaps wrong Dfun type).",
1245 -2: "Excess accuracy requested (tolerances too small).",
1246 -3: "Illegal input detected (internal error).",
1247 -4: "Repeated error test failures (internal error).",
1248 -5: "Repeated convergence failures (perhaps bad Jacobian or tolerances).",
1249 -6: "Error weight became zero during problem.",
1250 -7: "Internal workspace insufficient to finish (internal error)."
1251 }
1253 def __init__(self,
1254 with_jacobian=False,
1255 rtol=1e-6, atol=1e-12,
1256 lband=None, uband=None,
1257 nsteps=500,
1258 max_step=0.0, # corresponds to infinite
1259 min_step=0.0,
1260 first_step=0.0, # determined by solver
1261 ixpr=0,
1262 max_hnil=0,
1263 max_order_ns=12,
1264 max_order_s=5,
1265 method=None
1266 ):
1268 self.with_jacobian = with_jacobian
1269 self.rtol = rtol
1270 self.atol = atol
1271 self.mu = uband
1272 self.ml = lband
1274 self.max_order_ns = max_order_ns
1275 self.max_order_s = max_order_s
1276 self.nsteps = nsteps
1277 self.max_step = max_step
1278 self.min_step = min_step
1279 self.first_step = first_step
1280 self.ixpr = ixpr
1281 self.max_hnil = max_hnil
1282 self.success = 1
1284 self.initialized = False
1286 def reset(self, n, has_jac):
1287 # Calculate parameters for Fortran subroutine dvode.
1288 if has_jac:
1289 if self.mu is None and self.ml is None:
1290 jt = 1
1291 else:
1292 if self.mu is None:
1293 self.mu = 0
1294 if self.ml is None:
1295 self.ml = 0
1296 jt = 4
1297 else:
1298 if self.mu is None and self.ml is None:
1299 jt = 2
1300 else:
1301 if self.mu is None:
1302 self.mu = 0
1303 if self.ml is None:
1304 self.ml = 0
1305 jt = 5
1306 lrn = 20 + (self.max_order_ns + 4) * n
1307 if jt in [1, 2]:
1308 lrs = 22 + (self.max_order_s + 4) * n + n * n
1309 elif jt in [4, 5]:
1310 lrs = 22 + (self.max_order_s + 5 + 2 * self.ml + self.mu) * n
1311 else:
1312 raise ValueError('Unexpected jt=%s' % jt)
1313 lrw = max(lrn, lrs)
1314 liw = 20 + n
1315 rwork = zeros((lrw,), float)
1316 rwork[4] = self.first_step
1317 rwork[5] = self.max_step
1318 rwork[6] = self.min_step
1319 self.rwork = rwork
1320 iwork = zeros((liw,), _lsoda_int_dtype)
1321 if self.ml is not None:
1322 iwork[0] = self.ml
1323 if self.mu is not None:
1324 iwork[1] = self.mu
1325 iwork[4] = self.ixpr
1326 iwork[5] = self.nsteps
1327 iwork[6] = self.max_hnil
1328 iwork[7] = self.max_order_ns
1329 iwork[8] = self.max_order_s
1330 self.iwork = iwork
1331 self.call_args = [self.rtol, self.atol, 1, 1,
1332 self.rwork, self.iwork, jt]
1333 self.success = 1
1334 self.initialized = False
1336 def run(self, f, jac, y0, t0, t1, f_params, jac_params):
1337 if self.initialized:
1338 self.check_handle()
1339 else:
1340 self.initialized = True
1341 self.acquire_new_handle()
1342 args = [f, y0, t0, t1] + self.call_args[:-1] + \
1343 [jac, self.call_args[-1], f_params, 0, jac_params]
1344 y1, t, istate = self.runner(*args)
1345 self.istate = istate
1346 if istate < 0:
1347 unexpected_istate_msg = 'Unexpected istate={:d}'.format(istate)
1348 warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,
1349 self.messages.get(istate, unexpected_istate_msg)))
1350 self.success = 0
1351 else:
1352 self.call_args[3] = 2 # upgrade istate from 1 to 2
1353 self.istate = 2
1354 return y1, t
1356 def step(self, *args):
1357 itask = self.call_args[2]
1358 self.call_args[2] = 2
1359 r = self.run(*args)
1360 self.call_args[2] = itask
1361 return r
1363 def run_relax(self, *args):
1364 itask = self.call_args[2]
1365 self.call_args[2] = 3
1366 r = self.run(*args)
1367 self.call_args[2] = itask
1368 return r
1371if lsoda.runner:
1372 IntegratorBase.integrator_classes.append(lsoda)