sampledoc

scikits.bvp1lg

Abstract:Multi-point boundary value problem solver for Python

This is a Python wrapper for a modified version of the COLNEW boundary value problem solver by U. Ascher and G. Bader. Modifications made include vectorization over mesh points and better compatibility with Python.

It is written in the form of a Scikit. You should be able to compile it by unpacking it and typing the usual

cd scikits.bvp1lg-0.2.5
python setup.py build

The package requires recent versions of numpy and scipy.

Note

This package is still in development, and although it passes its test suite (consisting of the test problems for COLSYS and COLNEW), it has otherwise not been in regular use. Use with due caution and double-check your results.

Note

The licensing conditions for the COLNEW and MUS codes are rather unclear. So ask the authors of the codes listed in LICENSE.txt at least before using them commercially. Some additional information:

  • COLNEW is based on COLSYS, which is algorithm #569 in the Collected Algorithms of the ACM catalogue. These are covered by a restrictive non-commercial-use-only license. Nonetheless, COLNEW is used for example in Scilab, but which still has a restrictive license.

Files


Example 1

Consider the following problem:

../../_images/bessel-eq.png

The following Python code solves this:

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

import scikits.bvp1lg as 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)

# Analytical solution is of course known here, so check it

x = solution.mesh
assert N.allclose(solution(x)[:,0], N.special.jv(nu, x),
                  rtol=2e-3, atol=1e-8)
assert N.allclose(solution(x)[:,2], x**(nu+1)*N.special.jv(nu+1, x),
                  rtol=2e-3, atol=1e-8)

# Plot it

import pylab

pylab.plot(x, N.special.jv(nu, x), 'g.',
           x, solution(x)[:,0], 'b-')
pylab.title('$u$')
pylab.savefig('bessel-solution.png')
pylab.show()

pylab.plot(x, N.special.jv(nu+1, x)*x**(nu+1), 'g.',
           x, solution(x)[:,2], 'b-',)
pylab.title('$v$')
pylab.savefig('bessel-solution-v.png')
pylab.show()

The solutions look like this: u and v compared to the exact results

../../_images/bessel-solution.png ../../_images/bessel-solution-v.png

Example 2

../../_images/nonlin-eq.png

Second, a non-linear problem with a free parameter C For large C the solution becomes less and less smooth (see below).

The following Python program solves this using simple continuation, up to C = 500. I suspect shooting methods might not work here. Now, also define the jacobians of the functions to speed the code up:

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

C = None

# The exact solution is known

def exact_solution(x):
    return N.arccosh(N.sqrt((1.+C)/C)*N.cosh(N.sqrt(C)*(x - .5)))

def f(x, z):
    """Right-hand side"""
    u, du = z
    return N.array([N.cosh(u)/N.sinh(u)**3])

def df(x, z):
    """Jacobian of f"""
    u, du = z
    return N.array([[-(1 + 2*N.cosh(u)**2)/N.sinh(u)**4, 0*x]])

def g(z):
    """Boundary conditions"""
    u, du = z
    v0 = exact_solution(0)
    v1 = exact_solution(1)
    return N.array([u[0] - v0,
                    u[1] - v1])

def dg(z):
    """Jacobian of g"""
    return N.array([[1, 0],
                    [1, 0]])

def guess(x):
    z = N.zeros((2,) + N.asarray(x).shape)
    dm = N.zeros((1,) + N.asarray(x).shape)
    z[0] = 1
    return z, dm

# Solve the problem for multiple values of C using simple continuation.

Cvals = N.r_[N.linspace(1, 1.69, 5),
             N.linspace(1.69, 1.72, 10), # branch point here?
             N.linspace(1.8, 500, 40)]

x = N.linspace(0, 1, 501)

solution = guess
for Cval in Cvals:
    C = Cval
    print C
    solution = bvp.colnew.solve([0, 1], [2],
                                f, g,
                                dfsub=df, dgsub=dg,
                                initial_guess=solution,
                                tolerances=[1e-5, 1e-5],
                                maximum_mesh_size=1500,
                                vectorized=True,
                                coarsen_initial_guess_mesh=True)
    assert N.allclose(exact_solution(x), solution(x)[:,0], rtol=1e-5)

# Plot it

import pylab

x = solution.mesh

pylab.plot(x, exact_solution(x), 'g.',
           x, solution(x)[:,0], 'b-')
pylab.title('$u$')
pylab.savefig('nonlin-solution.png')
pylab.show()

Typical run time for the above on my machine (Athlon XP 2500) is 1.5s, when the plotting commands are commented out. Without the jacobians, the run time is 1.8s – the advantage of specifying them here is quite small. But note that this is not at all reliable as benchmark, the time includes Python startup.

The solution looks like the following. It approaches a function u(x) = A |x - 1/2| as C increases.

../../_images/nonlin-solution.png

COLNEW’s adaptive mesh refinement can be seen clearly in this picture.

More examples and a test suite can be found in the package.

Example 3

Example 3 of the Matlab BVP Tutorial: finding the fourth eigenvalue (lmd) of a Sturm-Liouville eigenvalue problem:

y'' - (lmd-2*q*cos(4*x))*y = 0

This shows how to solve a bvp-problem with an unknown parameter. In particular the parameter lmd is interpreted as an additional function with vanishing derivative.

import numpy as np
import scikits.bvp1lg as bvp
import matplotlib.pyplot as plt

# Adaptation courtesy of B. Weber

q = 5

def fsub(x, Y):
    """The equations, in the form: Y' = f(x, Y)."""
    y, dy, lmd = Y
    ddy=-y*(lmd-2*q*np.cos(2*x))
    # lmd shall be a constant -> (d lmd)/dx = 0 for all x
    dlmd=0*x # short for np.zeros_like(x)
    return np.array([ddy, dlmd])

def gsub(Y):
    """The boundary conditions."""
    y, dy, lmd = Y
    y_at_a = 1
    dy_at_a = 0
    dy_at_b = 0
    # obviously there is no bc for lmd
    return np.array([y[0] - y_at_a,
                     dy[1] - dy_at_a,
                     dy[2] - dy_at_b
                     ])

def guess(x):
    y  = np.cos(4*x)
    dy = -4*np.sin(4*x)
    lmd = np.ones_like(x)*15
    Y = np.array([y, dy, lmd])
    dm = fsub(x, Y)
    return Y, dm

# 2 times zero as we have two bcs at x=0 :
boundary_points = [0, 0, np.pi]
tol = 1e-3 * np.ones_like(boundary_points)
degrees = [2, 1]

solution = bvp.colnew.solve(
    boundary_points, degrees, fsub, gsub,
    is_linear=False, tolerances=tol, initial_guess=guess,
    vectorized=True, maximum_mesh_size=5000)

# plot
x_orig = solution.mesh
y_orig = solution(x_orig)[:,0]
plt.plot(x_orig, y_orig, 'g.')

# higher resolution
x = np.linspace(0, np.pi, 101)
y, dy, lmd = solution(x).transpose()
plt.plot(x, y, 'b-')

print "4th eigenvalue: %5.3f" % lmd[0]

plt.gca().annotate('$\lambda_4=%4.3f$' %lmd[0],
                   xycoords="figure fraction",
                   xy=(.7, 0.2), size=16,
                   xytext=None, arrowprops=None)
plt.savefig('ex3.png', dpi=60)
plt.show()
../../_images/ex3.png

More examples

See the following links: