colnew

Solve multi-point boundary value problems for ODEs

Description

This module uses a modified version of COLNEW [CN], a mature solver for multi-point boundary value problems for ODEs.

COLNEW handles only problems with separated boundary conditions. Non-separated problems can be converted to separated form for example by adding dummy variables:

\[u''(x) = f(x, u)\]\[u(0) + u(1) = 1\]\[u(0) u(1) = 2\]

can be transformed to

\[u''(x) = f(x, u)\]\[v'(x) = 0\]\[v(0) - u(0) = 0\]\[v(1) + u(1) = 1\]\[v(1) u(1) = 2\]

Similarly, problems with constant parameters

\[u''(x) + (a + 2 \cos(2 x)) u(x) = 0\]\[u'(0) = 0, u(0) = 1, u'(\pi) = 0\]

can be transformed to

\[u''(x) + (a(x) + 2 \cos(2 x)) u(x) = 0\]\[a'(x) = 0\]\[u'(0) = 0, u(0) = 1, u'(\pi) = 0\]

This may make the problem non-linear.

References

[CN]U. Ascher and G. Bader (and J. Christiansen and R. D. Russell). SIAM J. Sci. Comput. 8, 483 (1987). http://www.netlib.org/ode/colnew.f

Module contents

scikits.bvp1lg.colnew.DEBUG = 2

Print debug output

scikits.bvp1lg.colnew.INFO = 1

Print selected output

scikits.bvp1lg.colnew.REGULAR = 0

The problem is regular

scikits.bvp1lg.colnew.SENSITIVE = 1

The problem is sensitive. The nonlinear iteration should not rely on past covergence.

scikits.bvp1lg.colnew.SILENT = 0

Print no messages

class scikits.bvp1lg.colnew.Solution(ispace, fspace)

A solution to a boundary value problem for an ODE

Methods

get_mesh() Get the mesh points on which the solution is specified
get_mesh_values() Get the solution at the mesh points
fspace = None

The FSPACE vector provided by COLNEW

get_mesh()

Get the mesh points on which the solution is specified

Returns:mesh : ndarray of float, shape (nmesh,)
get_mesh_values()

Get the solution at the mesh points

Returns:

values : ndarray

self(self.mesh)

ispace = None

The ISPACE vector provided by COLNEW

mesh

The mesh on which the solution is specified

scikits.bvp1lg.colnew.check_jacobians(boundary_points, degrees, fsub, gsub, dfsub, dgsub, vectorized=True, **kw)

Check that the Jacobian functions match numerically evaluated derivatives.

Parameters:

degrees, boundary_points, fsub, dfsub, gsub, dsub, vectorized :

As for solve.

kw :

Passed on to jacobian.check_jacobian

Raises:

ValueError :

If the jacobians seem to be invalid.

scikits.bvp1lg.colnew.solve(boundary_points, degrees, fsub, gsub, dfsub=None, dgsub=None, left=None, right=None, is_linear=False, initial_guess=None, coarsen_initial_guess_mesh=True, initial_mesh=None, tolerances=None, adaptive_mesh_selection=True, verbosity=0, collocation_points=None, extra_fixed_points=None, problem_regularity=0, maximum_mesh_size=100, vectorized=True, is_complex=False)

Solve a multi-point boundary value problem for a system of ODEs.

The mixed-order system is:

ncomp = len(degrees)
mstar = sum(degrees)

1 <= min(degrees) <= max(degrees) <= 4

u_i^{(m_i)}(x) = f_i(x, z(x))       i = 0, ..., ncomp-1
                                    left <= x <= right

g_j(zeta_j, z(zeta_j)) = 0          j = 0, ..., mstar-1

where u(x) is the solution vector at position x and zeta = boundary_points specifies the boundary points.

The solution vector is represented by z-vector:

z = [u_1, u_1', ..., u_1^{m_1-1}, u_2, u_2', ..., u_{mstar-1}]

It is of shape (mstar,) and contains derivatives of orders < m_i.

Note

Colnew has the hard-coded problem size limits::
ncomp <= 256 mstar <= 512
Parameters:

boundary_points :

Points where i:th boundary condition is given, as a (mstar,) array, with left <= boundary_points[i] <= right for all i. It must be sorted in increasing order.

degrees : list of integers

Degree of i:th equation is degree[i]. It is required that 1 <= degree[i] <= 4.

fsub : callable

Function f, given as def fsub(x, z): return f, where:

x[j]    = x_j                             (nx,)
z[i, k] = z_i(x[k])                       (mstar, nx)
f[i, k] = f_i(x[k], z[:,k])               (ncomp, nx)

If not vectorized, the last dimension is omitted for all variables. The function must be local: f[:,k] can only depend on z[:,k].

gsub : callable

Function g, given as def gsub(z): return g, where:

z[i, j] = z_i(u(zeta_j))                  (mstar, mstar)
g[i]    = g_i(zeta_i, z(u(zeta_i)))       (mstar,)

Boundary conditions must be separated: g[:, j] may depend only on z[:, j].

dfsub : callable, optional

Jacobian of f, given as def dfsub(x, z): return df, where:

x[j]        = x_j                         (nx,)
z[j, k]     = z_j(x_k)                    (mstar, nx)
df[i, j, k] = d f[i,k] / d z[j,k]         (ncomp, mstar, nx)

If not vectorized, the last dimension is omitted for all variables. If None, a simple difference approximation is used.

dgsub : callable, optional

Jacobian of g, given as def dgsub(z): return dg, where:

z[i, j]  = z_i(u(zeta_j))                 (mstar, mstar)
dg[i, j] = (d g_i / d z_j)(zeta_i, z)     (mstar, mstar)

If None, a simple difference approximation is used.

left : float, optional

The left boundary point. If None, left = min(boundary_points).

right : float, optional

The right boundary point. If None, right = max(boundary_points).

vectorized : bool, optional

Are the functions fsub, dfsub and initial_guess vectorized?

is_linear : bool, optional

Is the system of equations linear?

initial_guess : callable or Solution, optional

Initial guess for continuation. Can be

  1. Callable def guess(x): return z, dm, where:

    x[j]     = x_j                     (nx,)
    z[i, j]  = z_i(u(x_j))             (mstar, nx)
    dm[i, k] = u_i^{m_i}(x_j)          (ncomp, nx)
    

    If not vectorized, the last dimension is omitted for all variables.

  2. Previously obtained Solution

  3. None, indicating that a default initial guess is to be used.

tolerances : list of float, optional

Tolerances for components of the solution. tolerance[i] gives the tolerance for i:th component of z-vector. If tolerance[i] == 0, then no tolerance is imposed for that component.

adaptive_mesh_selection : bool, optional

Use adaptive mesh selection. If disabled, trivial mesh refinement is used – in this case the initial mesh to use must be given in initial_mesh.

verbosity : int, optional

Amount of messages to show. 0 means silent, 1 selected printout, and 2 diagnostic printout.

collocation_points : int, optional

Number of collocation points in each subinterval, or None for a sensible default. It is required that max(degrees) <= collocation_points <= 7.

extra_fixed_points : list of float, optional

Points to fix in the mesh, in addition to boundary_points. (E.g. known boundary layers etc.)

problem_regularity : int, optional

How regular the problem is. Can be REGULAR (0) or SENSITIVE (1). Usually, SENSITIVE should not be needed.

maximum_mesh_size : int, optional

Maximum number of points to allow in the mesh.

is_complex : bool, optional

Whether the problem is complex-valued. The equation must be analytical in the unknown variables.

Returns:

sol : Solution

Object representing the solution.

Raises:

ValueError :

Invalid input

scikits.bvp1lg.NoConvergence :

Numerical convergence problems

scikits.bvp1lg.TooManySubintervals :

maximum_mesh_size too small to satisfy tolerances

scikits.bvp1lg.SingularCollocationMatrix :

Singular collocation matrix (check your jacobians)

SystemError :

Invalid output from user routines. (FIXME: these should be fixed)

Table Of Contents

Previous topic

scikits.bvp1lg

This Page