Theory

The first order ODE-IVP solver is from this paper:

• F. Stenger, B. Keys, M. O'Reilly, K.Parker, and S.-A Gustafson. “ODE-IVP-PACK via Sinc indefinite integration and Newton's method,” NUMERICAL ALGORITHMS, 20 (1999), PP. 241-268
  

We have used the same sinc discretization ideas, but apply a few optimizations where available and solve the resulting nonlinear equations through other means besides Newton's Method.
The ODE-IVP solver for a single equation solves this problem:
We can rewrite this in terms of an integral equation:
This has the following Frechet derivative:
Here, we can apply sinc indefinite integration to get the following vector equation:

The corresponding matrix equation for the Frechet derivative is:

We can use a nonlinear solver to solve these last two equations.  An iteration of Newton's method takes the form:

Implementation

The first-order, nonlinear ODE IVP solver is implemented in sinc.ode_system.DiffEqOptv2Lib with the DiffEqOptV2 class (it is in the ode_system module because it forms the base of the ODE IVP system solver).

The implementation relies on an underlying optimization packages provided by SciPy (http://scipy.org).  Given a nonlinear vector-valued function, t(f) and its Jacobian, t'(f), it will attempt to minimize the norm of t(f).  The following methods have been tested:

• Newton's Method
• BFGS
• L-BFGS (reduced-memory BFGS).  http://en.wikipedia.org/wiki/L-BFGS. 

See http://www.scipy.org/doc/api_docs/scipy.optimize.optimize.html for a full list of routines available; there are many other methods available we have not touched.

The following routine implements t(f):

  def T( self, x ):
    """ The integral operator corresponding to the differential equation. """
    x = array( x ) 
    # T(X) = X - I0 - ( h I^{-1} D( 1/phi' ) ) K
    retVal = x - self.init - self.h * dot(cSincLib.mdMultiply( self.I_n1(), self.df2 ), self.F( x ) )
    return retVal.reshape( (prod(retVal.shape),) )

This is a fairly straightforward implementation of the formula.

The function t'(f) is implemented in Tp:

  def Tp( self, x ):
    """ Computes the sinc approx. of the Frechet derivative of the operator T.
        From Stenger's paper:
        T'_{F} = I_1 x I_2 - I_1 x { h * I^{-1} * D(1/phi') } [ B_ij ]
    """ 
    x = array(x)
    I_n1xdf2 = -self.h * cSincLib.mdMultiply( self.I_n1() , self.df2 )
    if self.const_jacobian:
      j = self.j( self.nodes, x )
      return atleast_2d(j), I_n1xdf2
    else:
      j = ones( (self.m,), float64 ); j *= self.j( self.nodes, x )
      tmp2 = cSincLib.mdMultiply( I_n1xdf2, j )
      return cSincLib.mdAdd( tmp2, ones( self.m, float64 ) )