[SciPy-user] COLNEW the BVP solver for scipy [Was: Solving Two boundary value problems ...]

[Pauli Virtanen]

John Hassler wrote:
> I don't think that there's a function specifically for boundary value
> ode problems in SciPy, but it is pretty simple to combine an ode solver
> with a zero finder to use the "shooting method."  I'm not familiar with
> how rlab or Matlab does it, but the Scilab "bvode" uses a finite
> difference method.  If the problem isn't too complex, the shooting
> method is easy to set up, and works as well as the finite difference
> method.  (In my vast experience ... well, ok, two applications ... the
> shooting method was more forgiving and less brittle than the finite
> difference method.)

Hi all,

I've recently written a Python wrapper for a modified version of the COLNEW
bvp solver (Scilab's bvode is also based on COLNEW). The code can be found
from::

    http://www.elisanet.fi/ptvirtan/bvp/

In addition to wrapping COLNEW, the package has the following features

* I vectorized COLNEW over mesh points. (Vectorized code is approx 10x
  faster than non-vectorized, at least in my small tests.)
* There are simple estimators for the Jacobians so that they need not be
  explicitly provided. They appear to work well despite their simplicity.

This is still in development -- the code has so far not been in real use.
However, I made a test suite of the test problems for COLSYS, COLNEW and
MUS plus some others, and the code solves them correctly.

About usage and syntax: For example solving Bessel's equation plus another
one with some boundary conditions,

    u''(x) = -u'(x) / x + (nu**2/x**2 - 1) * u(x)        for 1 <= x <= 10
    v'(x) = x**(nu+1) * u(x)                             for 1 <= x <= 10

    u(1)  = J_{nu}(1)
    u(10) = J_{nu}(10)
    v(5)  = 5**(nu+1) * J_{nu+1}(5)

looks like this (yes, it is vectorized)

    import scipy as N
    N.pkgload('special')
    import bvp

    nu = 3.4123
    degrees = [2, 1]

    def fsub(x, z):
        """The equations"""
        u, du, v = z     # it's neat to name the variables
        return N.array([-du/x + (nu**2/x**2 - 1)*u, x**(nu+1) * u])

    def gsub(z):
        """The boundary conditions"""
        u, du, v = z
        return N.array([u[0] - N.special.jv(nu,   1),
                        v[1] - 5**(nu+1) * N.special.jv(nu+1, 5),
                        u[2] - N.special.jv(nu,   10)])

    boundary_points = [1, 5, 10]
    tol = [1e-5, 0, 1e-5]

    solution = bvp.colnew.solve(
        boundary_points, degrees, fsub, gsub,
        is_linear=True, tolerances=tol,
        vectorized=True, maximum_mesh_size=300)

    solution_values = solution([2,3,4,5,6])

I hope this package will be of some use to people. I'd be interested to hear
if it proves (un)useful or any other suggestions the people on the list
might have.

        Pauli Virtanen