[Scipy-svn] r3950 - in trunk/scipy/integrate: . tests
scipy-svn@scip...
scipy-svn@scip...
Tue Feb 19 23:27:20 CST 2008
Author: oliphant
Date: 2008-02-19 23:27:16 -0600 (Tue, 19 Feb 2008)
New Revision: 3950
Modified:
trunk/scipy/integrate/info.py
trunk/scipy/integrate/ode.py
trunk/scipy/integrate/tests/test_integrate.py
trunk/scipy/integrate/vode.pyf
Log:
Apply patch for ticket #334 which adds zvode to scipy's ode functionality.
Modified: trunk/scipy/integrate/info.py
===================================================================
--- trunk/scipy/integrate/info.py 2008-02-20 01:35:04 UTC (rev 3949)
+++ trunk/scipy/integrate/info.py 2008-02-20 05:27:16 UTC (rev 3950)
@@ -27,7 +27,7 @@
Interface to numerical integrators of ODE systems.
odeint -- General integration of ordinary differential equations.
- ode -- Integrate ODE using vode routine.
+ ode -- Integrate ODE using VODE and ZVODE routines.
"""
Modified: trunk/scipy/integrate/ode.py
===================================================================
--- trunk/scipy/integrate/ode.py 2008-02-20 01:35:04 UTC (rev 3949)
+++ trunk/scipy/integrate/ode.py 2008-02-20 05:27:16 UTC (rev 3950)
@@ -1,140 +1,208 @@
-## Automatically adapted for scipy Oct 21, 2005 by
+# Authors: Pearu Peterson, Pauli Virtanen
+"""
+First-order ODE integrators
-#!/usr/bin/env python
-#Author: Pearu Peterson
-#Date: 3 Feb 2002
-#$Revision$
-"""
User-friendly interface to various numerical integrators for solving a
-system of first order ODEs with prescribed initial conditions:
+system of first order ODEs with prescribed initial conditions::
- d y(t)[i]
- --------- = f(t,y(t))[i],
- d t
+ d y(t)[i]
+ --------- = f(t,y(t))[i],
+ d t
+
+ y(t=0)[i] = y0[i],
- y(t=0)[i] = y0[i],
+where::
-where i = 0, ..., len(y0) - 1
+ i = 0, ..., len(y0) - 1
-Provides:
- ode - a generic interface class to numeric integrators. It has the
- following methods:
- integrator = ode(f,jac=None)
- integrator = integrator.set_integrator(name,**params)
- integrator = integrator.set_initial_value(y0,t0=0.0)
- integrator = integrator.set_f_params(*args)
- integrator = integrator.set_jac_params(*args)
- y1 = integrator.integrate(t1,step=0,relax=0)
- flag = integrator.successful()
+class ode
+---------
-Supported integrators:
- vode - Variable-coefficient Ordinary Differential Equation solver,
- with fixed-leading-coefficient implementation.
- It provides implicit Adams method (for non-stiff problems)
- and a method based on backward differentiation formulas (BDF)
- (for stiff problems).
- Source: http://www.netlib.org/ode/vode.f
- This integrator accepts the following parameters in
- set_integrator() method of the ode class:
- atol=float|seq
- rtol=float|seq
- lband=None|int
- rband=None|int
- method='adams'|'bdf'
- with_jacobian=0|1
- nsteps = int
- (first|min|max)_step = float
- order = int # <=12 for adams, <=5 for bdf
+A generic interface class to numeric integrators. It has the following
+methods::
+
+ integrator = ode(f,jac=None)
+ integrator = integrator.set_integrator(name,**params)
+ integrator = integrator.set_initial_value(y0,t0=0.0)
+ integrator = integrator.set_f_params(*args)
+ integrator = integrator.set_jac_params(*args)
+ y1 = integrator.integrate(t1,step=0,relax=0)
+ flag = integrator.successful()
+
"""
+
+integrator_info = \
"""
-XXX: Integrators must have:
-===========================
-cvode - C version of vode and vodpk with many improvements.
- Get it from http://www.netlib.org/ode/cvode.tar.gz
- To wrap cvode to Python, one must write extension module by
- hand. Its interface is too much 'advanced C' that using f2py
- would be too complicated (or impossible).
+Available integrators
+---------------------
-How to define a new integrator:
-===============================
+vode
+~~~~
-class myodeint(IntegratorBase):
+Real-valued Variable-coefficient Ordinary Differential Equation
+solver, with fixed-leading-coefficient implementation. It provides
+implicit Adams method (for non-stiff problems) and a method based on
+backward differentiation formulas (BDF) (for stiff problems).
- runner = <odeint function> or None
+Source: http://www.netlib.org/ode/vode.f
- def __init__(self,...): # required
- <initialize>
+This integrator accepts the following parameters in set_integrator()
+method of the ode class:
- def reset(self,n,has_jac): # optional
- # n - the size of the problem (number of equations)
- # has_jac - whether user has supplied its own routine for Jacobian
- <allocate memory,initialize further>
+- atol : float or sequence
+ absolute tolerance for solution
+- rtol : float or sequence
+ relative tolerance for solution
+- lband : None or int
+- rband : None or int
+ Jacobian band width, jac[i,j] != 0 for i-lband <= j <= i+rband.
+ Setting these requires your jac routine to return the jacobian
+ in packed format, jac_packed[i-j+lband, j] = jac[i,j].
+- method: 'adams' or 'bdf'
+ Which solver to use, Adams (non-stiff) or BDF (stiff)
+- with_jacobian : bool
+ Whether to use the jacobian
+- nsteps : int
+ Maximum number of (internally defined) steps allowed during one
+ call to the solver.
+- first_step : float
+- min_step : float
+- max_step : float
+ Limits for the step sizes used by the integrator.
+- order : int
+ Maximum order used by the integrator,
+ order <= 12 for Adams, <= 5 for BDF.
- def run(self,f,jac,y0,t0,t1,f_params,jac_params): # required
- # this method is called to integrate from t=t0 to t=t1
- # with initial condition y0. f and jac are user-supplied functions
- # that define the problem. f_params,jac_params are additional arguments
- # to these functions.
- <calculate y1>
- if <calculation was unsuccesful>:
- self.success = 0
- return t1,y1
+zvode
+~~~~~
- # In addition, one can define step() and run_relax() methods (they
- # take the same arguments as run()) if the integrator can support
- # these features (see IntegratorBase doc strings).
+Complex-valued Variable-coefficient Ordinary Differential Equation
+solver, with fixed-leading-coefficient implementation. It provides
+implicit Adams method (for non-stiff problems) and a method based on
+backward differentiation formulas (BDF) (for stiff problems).
-if myodeint.runner:
- IntegratorBase.integrator_classes.append(myodeint)
+Source: http://www.netlib.org/ode/zvode.f
+
+This integrator accepts the same parameters in set_integrator()
+as the "vode" solver.
+
+:Note:
+ When using ZVODE for a stiff system, it should only be used for
+ the case in which the function f is analytic, that is, when each f(i)
+ is an analytic function of each y(j). Analyticity means that the
+ partial derivative df(i)/dy(j) is a unique complex number, and this
+ fact is critical in the way ZVODE solves the dense or banded linear
+ systems that arise in the stiff case. For a complex stiff ODE system
+ in which f is not analytic, ZVODE is likely to have convergence
+ failures, and for this problem one should instead use DVODE on the
+ equivalent real system (in the real and imaginary parts of y).
+
"""
+__doc__ += integrator_info
+
+# XXX: Integrators must have:
+# ===========================
+# cvode - C version of vode and vodpk with many improvements.
+# Get it from http://www.netlib.org/ode/cvode.tar.gz
+# To wrap cvode to Python, one must write extension module by
+# hand. Its interface is too much 'advanced C' that using f2py
+# would be too complicated (or impossible).
+#
+# How to define a new integrator:
+# ===============================
+#
+# class myodeint(IntegratorBase):
+#
+# runner = <odeint function> or None
+#
+# def __init__(self,...): # required
+# <initialize>
+#
+# def reset(self,n,has_jac): # optional
+# # n - the size of the problem (number of equations)
+# # has_jac - whether user has supplied its own routine for Jacobian
+# <allocate memory,initialize further>
+#
+# def run(self,f,jac,y0,t0,t1,f_params,jac_params): # required
+# # this method is called to integrate from t=t0 to t=t1
+# # with initial condition y0. f and jac are user-supplied functions
+# # that define the problem. f_params,jac_params are additional
+# # arguments
+# # to these functions.
+# <calculate y1>
+# if <calculation was unsuccesful>:
+# self.success = 0
+# return t1,y1
+#
+# # In addition, one can define step() and run_relax() methods (they
+# # take the same arguments as run()) if the integrator can support
+# # these features (see IntegratorBase doc strings).
+#
+# if myodeint.runner:
+# IntegratorBase.integrator_classes.append(myodeint)
+
__all__ = ['ode']
__version__ = "$Id$"
+__docformat__ = "restructuredtext en"
from numpy import asarray, array, zeros, sin, int32, isscalar
import re, sys
+#------------------------------------------------------------------------------
+# User interface
+#------------------------------------------------------------------------------
+
class ode(object):
"""\
-ode - a generic interface class to numeric integrators. It has the
- following methods:
- integrator = ode(f,jac=None)
- integrator = integrator.set_integrator(name,**params)
- integrator = integrator.set_initial_value(y0,t0=0.0)
- integrator = integrator.set_f_params(*args)
- integrator = integrator.set_jac_params(*args)
- y1 = integrator.integrate(t1,step=0,relax=0)
- flag = integrator.successful()
+A generic interface class to numeric integrators.
- Typical usage:
- r = ode(f,jac).set_integrator('vode').set_initial_value(y0,t0)
- t1 = <final t>
- dt = <step>
- while r.successful() and r.t < t1:
- r.integrate(r.t+dt)
- print r.t, r.y
- where f and jac have the following signatures:
- def f(t,y[,arg1,..]):
- return <f(t,y)>
- def jac(t,y[,arg1,..]):
- return <df/dy(t,y)>
+See also
+--------
+odeint : an integrator with a simpler interface based on lsoda from ODEPACK
+quad : for finding the area under a curve
-See also:
- odeint - an integrator with a simpler interface based on lsoda from ODEPACK
- quad - for finding the area under a curve
- """
+Examples
+--------
+A problem to integrate and the corresponding jacobian:
- def __init__(self,f,jac=None):
- """Define equation y' = f(y,t) where (optional) jac = df/dy.
- User-supplied functions must have the following signatures:
- def f(t,y,...):
- return <f(t,y)>
- def jac(t,y,...):
- return <jac(t,y)>
- where ... means extra parameters that can be set with
- set_(f|jac)_params(*args)
- methods.
+>>> from scipy import eye
+>>> from scipy.integrate import ode
+>>>
+>>> y0, t0 = [1.0j, 2.0], 0
+>>>
+>>> def f(t, y, arg1):
+>>> return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
+>>> def jac(t, y, arg1):
+>>> return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
+
+The integration:
+
+>>> r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
+>>> r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
+>>> t1 = 10
+>>> dt = 1
+>>> while r.successful() and r.t < t1:
+>>> r.integrate(r.t+dt)
+>>> print r.t, r.y
+
+"""
+
+ __doc__ += integrator_info
+
+ def __init__(self, f, jac=None):
"""
+ Define equation y' = f(y,t) where (optional) jac = df/dy.
+
+ Parameters
+ ----------
+ f : f(t, y, *f_args)
+ Rhs of the equation. t is a scalar, y.shape == (n,).
+ f_args is set by calling set_f_params(*args)
+ jac : jac(t, y, *jac_args)
+ Jacobian of the rhs, jac[i,j] = d f[i] / d y[j]
+ jac_args is set by calling set_f_params(*args)
+ """
self.stiff = 0
self.f = f
self.jac = jac
@@ -142,20 +210,29 @@
self.jac_params = ()
self.y = []
- def set_initial_value(self,y,t=0.0):
+ def set_initial_value(self, y, t=0.0):
"""Set initial conditions y(t) = y."""
if isscalar(y):
y = [y]
n_prev = len(self.y)
- self.y = asarray(y, float)
- self.t = t
if not n_prev:
self.set_integrator('') # find first available integrator
+ self.y = asarray(y, self._integrator.scalar)
+ self.t = t
self._integrator.reset(len(self.y),self.jac is not None)
return self
- def set_integrator(self,name,**integrator_params):
- """Set integrator by name."""
+ def set_integrator(self, name, **integrator_params):
+ """
+ Set integrator by name.
+
+ Parameters
+ ----------
+ name : str
+ Name of the integrator
+ integrator_params
+ Additional parameters for the integrator.
+ """
integrator = find_integrator(name)
if integrator is None:
print 'No integrator name match with %s or is not available.'\
@@ -164,11 +241,11 @@
self._integrator = integrator(**integrator_params)
if not len(self.y):
self.t = 0.0
- self.y = array([0.0], float)
+ self.y = array([0.0], self._integrator.scalar)
self._integrator.reset(len(self.y),self.jac is not None)
return self
- def integrate(self,t,step=0,relax=0):
+ def integrate(self, t, step=0, relax=0):
"""Find y=y(t), set y as an initial condition, and return y."""
if step and self._integrator.supports_step:
mth = self._integrator.step
@@ -188,23 +265,22 @@
return self._integrator.success==1
def set_f_params(self,*args):
- """Set extra-parameters for user-supplied function f."""
+ """Set extra parameters for user-supplied function f."""
self.f_params = args
return self
def set_jac_params(self,*args):
- """Set extra-parameters for user-supplied function jac."""
+ """Set extra parameters for user-supplied function jac."""
self.jac_params = args
return self
-#############################################################
-#### Nothing interesting for an end-user in what follows ####
-#############################################################
+#------------------------------------------------------------------------------
+# ODE integrators
+#------------------------------------------------------------------------------
def find_integrator(name):
for cl in IntegratorBase.integrator_classes:
if re.match(name,cl.__name__,re.I):
- print 'Found integrator',cl.__name__
return cl
return
@@ -215,6 +291,7 @@
supports_run_relax = None
supports_step = None
integrator_classes = []
+ scalar = float
def reset(self,n,has_jac):
"""Prepare integrator for call: allocate memory, set flags, etc.
@@ -379,47 +456,107 @@
IntegratorBase.integrator_classes.append(vode)
-def test1():
- def f(t,y):
- a = sin(6*t)
- return y*y-a+y
+class zvode(vode):
+ try:
+ import vode as _vode
+ except ImportError:
+ print sys.exc_value
+ _vode = None
+ runner = getattr(_vode,'zvode',None)
- ode_runner = ode(f)
- ode_runner.set_integrator('vode')
- ode_runner.set_initial_value([0.1,0.11,.1]*10)
+ supports_run_relax = 1
+ supports_step = 1
+ scalar = complex
- while ode_runner.successful() and ode_runner.t < 50:
- y1 = ode_runner.integrate(ode_runner.t+2)
- print ode_runner.t,y1[:3]
+ def reset(self, n, has_jac):
+ # Calculate parameters for Fortran subroutine dvode.
+ if has_jac:
+ if self.mu is None and self.ml is None:
+ miter = 1
+ else:
+ if self.mu is None: self.mu = 0
+ if self.ml is None: self.ml = 0
+ miter = 4
+ else:
+ if self.mu is None and self.ml is None:
+ if self.with_jacobian:
+ miter = 2
+ else:
+ miter = 0
+ else:
+ if self.mu is None: self.mu = 0
+ if self.ml is None: self.ml = 0
+ if self.ml==self.mu==0:
+ miter = 3
+ else:
+ miter = 5
-def test2():
- # Stiff problem. Requires analytic Jacobian.
- def f(t,y):
- ydot0 = -0.04*y[0] + 1e4*y[1]*y[2]
- ydot2 = 3e7*y[1]*y[1]
- ydot1 = -ydot0-ydot2
- return [ydot0,ydot1,ydot2]
- def jac(t,y):
- jc = [[-0.04,1e4*y[2] ,1e4*y[1]],
- [0.04 ,-1e4*y[2]-6e7*y[1],-1e4*y[1]],
- [0.0 ,6e7*y[1] ,0.0]]
- return jc
- r = ode(f,jac).set_integrator('vode',
- rtol=1e-4,
- atol=[1e-8,1e-14,1e-6],
- method='bdf',
- )
- r.set_initial_value([1,0,0])
- print 'At t=%s y=%s'%(r.t,r.y)
- tout = 0.4
- for i in range(12):
- r.integrate(tout)
- print 'At t=%s y=%s'%(r.t,r.y)
- tout *= 10
+ mf = 10*self.meth + miter
-if __name__ == "__main__":
- print 'Integrators available:',\
- ', '.join(map(lambda c:c.__name__,
- IntegratorBase.integrator_classes))
- test1()
- test2()
+ if mf in (10,):
+ lzw = 15*n
+ elif mf in (11, 12):
+ lzw = 15*n + 2*n**2
+ elif mf in (-11, -12):
+ lzw = 15*n + n**2
+ elif mf in (13,):
+ lzw = 16*n
+ elif mf in (14,15):
+ lzw = 17*n + (3*self.ml + 2*self.mu)*n
+ elif mf in (-14,-15):
+ lzw = 16*n + (2*self.ml + self.mu)*n
+ elif mf in (20,):
+ lzw = 8*n
+ elif mf in (21, 22):
+ lzw = 8*n + 2*n**2
+ elif mf in (-21,-22):
+ lzw = 8*n + n**2
+ elif mf in (23,):
+ lzw = 9*n
+ elif mf in (24, 25):
+ lzw = 10*n + (3*self.ml + 2*self.mu)*n
+ elif mf in (-24, -25):
+ lzw = 9*n + (2*self.ml + self.mu)*n
+
+ lrw = 20 + n
+
+ if miter in (0, 3):
+ liw = 30
+ else:
+ liw = 30 + n
+
+ zwork = zeros((lzw,), complex)
+ self.zwork = zwork
+
+ rwork = zeros((lrw,), float)
+ rwork[4] = self.first_step
+ rwork[5] = self.max_step
+ rwork[6] = self.min_step
+ self.rwork = rwork
+
+ iwork = zeros((liw,), int32)
+ if self.ml is not None:
+ iwork[0] = self.ml
+ if self.mu is not None:
+ iwork[1] = self.mu
+ iwork[4] = self.order
+ iwork[5] = self.nsteps
+ iwork[6] = 2 # mxhnil
+ self.iwork = iwork
+
+ self.call_args = [self.rtol,self.atol,1,1,
+ self.zwork,self.rwork,self.iwork,mf]
+ self.success = 1
+
+ def run(self,*args):
+ y1,t,istate = self.runner(*(args[:5]+tuple(self.call_args)+args[5:]))
+ if istate < 0:
+ print 'zvode:', self.messages.get(istate,
+ 'Unexpected istate=%s'%istate)
+ self.success = 0
+ else:
+ self.call_args[3] = 2 # upgrade istate from 1 to 2
+ return y1, t
+
+if zvode.runner:
+ IntegratorBase.integrator_classes.append(zvode)
Modified: trunk/scipy/integrate/tests/test_integrate.py
===================================================================
--- trunk/scipy/integrate/tests/test_integrate.py 2008-02-20 01:35:04 UTC (rev 3949)
+++ trunk/scipy/integrate/tests/test_integrate.py 2008-02-20 05:27:16 UTC (rev 3950)
@@ -1,52 +1,146 @@
-#!/usr/bin/env python
-
-# Test provided by Nils Wagner.
-# File created by Ed Schofield on Nov 16.
-
-""" Tests for numerical integration.
+# Authors: Nils Wagner, Ed Schofield, Pauli Virtanen
"""
+Tests for numerical integration.
+"""
import numpy
-from numpy import arange, zeros, array, dot, sqrt, cos, sin
+from numpy import (arange, zeros, array, dot, sqrt, cos, sin, absolute,
+ eye, pi, exp, allclose)
from scipy.linalg import norm
+
from scipy.testing import *
-from scipy.integrate import odeint
+from scipy.integrate import odeint, ode
-class TestODEInt(TestCase):
- """ Test odeint: free vibration of a simple oscillator
+#------------------------------------------------------------------------------
+# Test ODE integrators
+#------------------------------------------------------------------------------
+
+class TestOdeint(TestCase):
+ """
+ Check integrate.odeint
+ """
+ def _do_problem(self, problem):
+ t = arange(0.0, problem.stop_t, 0.05)
+ z, infodict = odeint(problem.f, problem.z0, t, full_output=True)
+ assert problem.verify(z, t)
+
+ def test_odeint(self):
+ for problem_cls in PROBLEMS:
+ problem = problem_cls()
+ if problem.cmplx: continue
+ self._do_problem(problem)
+
+class TestOde(TestCase):
+ """
+ Check integrate.ode
+ """
+ def _do_problem(self, problem, integrator, method='adams'):
+
+ # ode has callback arguments in different order than odeint
+ f = lambda t, z: problem.f(z, t)
+ jac = None
+ if hasattr(problem, 'jac'):
+ jac = lambda t, z: problem.jac(z, t)
+
+ ig = ode(f, jac)
+ ig.set_integrator(integrator,
+ atol=problem.atol/10,
+ rtol=problem.rtol/10,
+ method=method)
+ ig.set_initial_value(problem.z0, t=0.0)
+ z = ig.integrate(problem.stop_t)
+
+ assert ig.successful(), (problem, method)
+ assert problem.verify(array([z]), problem.stop_t), (problem, method)
+
+ def test_vode(self):
+ """Check the vode solver"""
+ for problem_cls in PROBLEMS:
+ problem = problem_cls()
+ if problem.cmplx: continue
+ if not problem.stiff:
+ self._do_problem(problem, 'vode', 'adams')
+ self._do_problem(problem, 'vode', 'bdf')
+
+ def test_zvode(self):
+ """Check the zvode solver"""
+ for problem_cls in PROBLEMS:
+ problem = problem_cls()
+ if not problem.stiff:
+ self._do_problem(problem, 'zvode', 'adams')
+ self._do_problem(problem, 'zvode', 'bdf')
+
+#------------------------------------------------------------------------------
+# Test problems
+#------------------------------------------------------------------------------
+
+class ODE:
+ """
+ ODE problem
+ """
+ stiff = False
+ cmplx = False
+ stop_t = 1
+ z0 = []
+
+ atol = 1e-6
+ rtol = 1e-5
+
+class SimpleOscillator(ODE):
+ r"""
+ Free vibration of a simple oscillator::
m \ddot{u} + k u = 0, u(0) = u_0 \dot{u}(0) \dot{u}_0
-
- Solution:
+ Solution::
u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m)
"""
+ stop_t = 1 + 0.09
+ z0 = array([1.0, 0.1], float)
- def setUp(self):
- self.k = 4.0
- self.m = 1.0
+ k = 4.0
+ m = 1.0
- def F(self, z, t):
+ def f(self, z, t):
tmp = zeros((2,2), float)
tmp[0,1] = 1.0
tmp[1,0] = -self.k / self.m
- return dot(tmp,z)
+ return dot(tmp, z)
- def test_odeint1(self):
+ def verify(self, zs, t):
omega = sqrt(self.k / self.m)
- z0 = zeros(2, float)
- z0[0] = 1.0 # initial displacement
- z0[1] = 0.1 # initial velocity
- t = arange(0.0, 1+0.09, 0.1)
+ u = self.z0[0]*cos(omega*t)+self.z0[1]*sin(omega*t)/omega
+ return allclose(u, zs[:,0], atol=self.atol, rtol=self.rtol)
- # Analytical solution
- #
- u = z0[0]*cos(omega*t)+z0[1]*sin(omega*t)/omega
+class ComplexExp(ODE):
+ r"""The equation :lm:`\dot u = i u`"""
+ stop_t = 1.23*pi
+ z0 = exp([1j,2j,3j,4j,5j])
+ cmplx = True
- # Numerical solution
- z, infodict = odeint(self.F, z0, t, full_output=True)
+ def f(self, z, t):
+ return 1j*z
- res = norm(u - z[:,0])
- print 'Residual:', res
- assert res < 1.0e-6
+ def jac(self, z, t):
+ return 1j*eye(5)
+ def verify(self, zs, t):
+ u = self.z0 * exp(1j*t)
+ return allclose(u, zs, atol=self.atol, rtol=self.rtol)
+
+class Pi(ODE):
+ r"""Integrate 1/(t + 1j) from t=-10 to t=10"""
+ stop_t = 20
+ z0 = [0]
+ cmplx = True
+
+ def f(self, z, t):
+ return array([1./(t - 10 + 1j)])
+ def verify(self, zs, t):
+ u = -2j*numpy.arctan(10)
+ return allclose(u, zs[-1,:], atol=self.atol, rtol=self.rtol)
+
+PROBLEMS = [SimpleOscillator, ComplexExp, Pi]
+
+#------------------------------------------------------------------------------
+
if __name__ == "__main__":
nose.run(argv=['', __file__])
Modified: trunk/scipy/integrate/vode.pyf
===================================================================
--- trunk/scipy/integrate/vode.pyf 2008-02-20 01:35:04 UTC (rev 3949)
+++ trunk/scipy/integrate/vode.pyf 2008-02-20 05:27:16 UTC (rev 3950)
@@ -25,12 +25,35 @@
end subroutine jac
end interface
end python module dvode__user__routines
+
+python module zvode__user__routines
+ interface zvode_user_interface
+ subroutine f(n,t,y,ydot,rpar,ipar)
+ integer intent(hide) :: n
+ double precision intent(in) :: t
+ double complex dimension(n),intent(in,c) :: y
+ double complex dimension(n),intent(out,c) :: ydot
+ double precision intent(hide) :: rpar
+ integer intent(hide) :: ipar
+ end subroutine f
+ subroutine jac(n,t,y,ml,mu,jac,nrowpd,rpar,ipar)
+ integer intent(hide) :: n
+ double precision :: t
+ double complex dimension(n),intent(c,in) :: y
+ integer intent(hide) :: ml,mu
+ integer intent(hide):: nrowpd
+ double complex intent(out) :: jac(nrowpd, n)
+ double precision intent(hide) :: rpar
+ integer intent(hide) :: ipar
+ end subroutine jac
+ end interface
+end python module zvode__user__routines
python module vode
interface
subroutine dvode(f,jac,neq,y,t,tout,itol,rtol,atol,itask,istate,iopt,rwork,lrw,iwork,liw,mf,rpar,ipar)
! y1,t,istate = dvode(f,jac,y0,t0,t1,rtol,atol,itask,istate,rwork,iwork,mf)
- callstatement (*f2py_func)(cb_f_in_dvode__user__routines,&neq,y,&t,&tout,&itol,&rtol,atol,&itask,&istate,&iopt,rwork,&lrw,iwork,&liw,cb_jac_in_dvode__user__routines,&mf,&rpar,&ipar)
+ callstatement (*f2py_func)(cb_f_in_dvode__user__routines,&neq,y,&t,&tout,&itol,rtol,atol,&itask,&istate,&iopt,rwork,&lrw,iwork,&liw,cb_jac_in_dvode__user__routines,&mf,&rpar,&ipar)
use dvode__user__routines
external f
external jac
@@ -56,4 +79,36 @@
integer intent(hide) :: ipar = 0
end subroutine dvode
end interface
+
+ interface
+ subroutine zvode(f,jac,neq,y,t,tout,itol,rtol,atol,itask,istate,iopt,zwork,lzw,rwork,lrw,iwork,liw,mf,rpar,ipar)
+ ! y1,t,istate = zvode(f,jac,y0,t0,t1,rtol,atol,itask,istate,rwork,iwork,mf)
+ callstatement (*f2py_func)(cb_f_in_zvode__user__routines,&neq,y,&t,&tout,&itol,rtol,atol,&itask,&istate,&iopt,zwork,&lzw,rwork,&lrw,iwork,&liw,cb_jac_in_zvode__user__routines,&mf,&rpar,&ipar)
+ use zvode__user__routines
+ external f
+ external jac
+
+ integer intent(hide),depend(y) :: neq = len(y)
+ double complex dimension(neq),intent(in,out,copy) :: y
+ double precision intent(in,out):: t
+ double precision intent(in):: tout
+ integer intent(hide),depend(atol) :: itol = (len(atol)<=1 && len(rtol)<=1?1:(len(rtol)<=1?2:(len(atol)<=1?3:4)))
+ double precision dimension(*),intent(in),check(len(atol)<&
+ &=1||len(atol)>=neq),depend(neq) :: atol
+ double precision dimension(*),intent(in),check(len(rtol)<&
+ &=1||len(rtol)>=neq),depend(neq) :: rtol
+ integer intent(in),check(itask>0 && itask<6) :: itask
+ integer intent(in,out),check(istate>0 && istate<4) :: istate
+ integer intent(hide) :: iopt = 1
+ double complex dimension(lzw),intent(in,cache) :: zwork
+ integer intent(hide),check(len(zwork)>=lzw),depend(zwork) :: lzw=len(zwork)
+ double precision dimension(lrw),intent(in,cache) :: rwork
+ integer intent(hide),check(len(rwork)>=lrw),depend(rwork) :: lrw=len(rwork)
+ integer dimension(liw),intent(in,cache) :: iwork
+ integer intent(hide),check(len(iwork)>=liw),depend(iwork) :: liw=len(iwork)
+ integer intent(in) :: mf
+ double precision intent(hide) :: rpar = 0.0
+ integer intent(hide) :: ipar = 0
+ end subroutine zvode
+ end interface
end python module vode
More information about the Scipy-svn
mailing list