Files

• scikits.bvp1lg-0.2.5.tar.gz: Support complex analytic ODEs. Repackage this as a Scikit.
• bvp-0.2.4.tar.gz: Fix for 64bit platforms.
• bvp-0.2.3.tar.gz: Support recursive calls into COLNEW (provided your Fortran compiler handles recursive routines properly).
• bvp-0.2.2.tar.gz: Try to fix more build issues.
• bvp-0.2.1.tar.gz: Fix bugs in setup script, speed up evaluation of the solution.
• bvp-0.2.tar.gz
• bvp-0.1.tar.gz
• Mercurial repository: Development version

Example 1

Consider the following problem:

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

Example 2

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.

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()

More examples

See the following links:

• Lorenzo Bolla's blog: http://lbolla.wordpress.com/2008/04/14/bvp/