|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.
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.
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:
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 - N.special.jv(nu, 1), v - 5**(nu+1) * N.special.jv(nu+1, 5), u - 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 - v0, u - 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 = 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], , 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 - y_at_a, dy - dy_at_a, dy - 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 plt.gca().annotate('$\lambda_4=%4.3f$' %lmd, xycoords="figure fraction", xy=(.7, 0.2), size=16, xytext=None, arrowprops=None) plt.savefig('ex3.png', dpi=60) plt.show()