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.

- 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

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

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

See the following links:

- Lorenzo Bolla’s blog: http://lbolla.wordpress.com/2008/04/14/bvp/