[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