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

Volker Lorrmann lorrmann at physik.uni-wuerzburg.de
Sat Oct 7 09:04:00 CDT 2006


That looks really amazing. Exactly what i´ve been looking for. Thanks!!

Volker
> 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



More information about the SciPy-user mailing list