mus

Solve two-point boundary value problems for ODEs

• solve_linear: Solve linear problems
• solve_nonlinear: Solve non-linear problems

Description

This module uses the multiple shooting code MUS [MS] for solving boundary-value problems.

Warning

For some reason, this code can be slow for certain problems. This may be due to that

1. This wrapper uses MUS somehow incorrectly
2. Speed of MUS is sensitive to tweaking parameters
3. Python call overhead is significant for MUS as it is not vectorized over mesh points

I have not yet figured out what the problem is, so you may be better off using the colnew package.

References

[MS]	R. M. M. Mattheij, G. W. M. Staarink. EUT Report 92-WSK-01 (1992).

Module contents

scikits.bvp1lg.mus.solve_linear(f_homogenous, f_nonhomogenous, a, b, m_a, m_b, bcv, max_amplification=None, rtol=None, atol=None, output_points=None, verbosity=0)

Solve a linear two-point boundary value problem.

The problem is assumed to be:

u'(t) = L(t) u(t) + r(t),        a <= t <= b
M_A u(a) + M_B u(b) = BCV

where u is the solution vector.

Parameters:

• f_homogenous: Function f_h(t, u) that evaluates the homogenous part L(t) u. It should return an array of shape (n,).
• f_nonhomogenous: Function f(t, u) that evaluates L(t) u + r(t). For homogenous problems, you can set f_nonhomogenous to None. It should return an array of shape (n,).
• a: Left boundary point.
• b: Right boundary point.
• m_a: The (n, n) M_A matrix in the boundary condition.
• m_b: The (n, n) M_B matrix in the boundary condition.
• bcv: The (n,) BCV vector in the boundary condition.
• max_amplification: The allowed incremental factor of the homogeneous solutions between two succesive output points. If the increment of a homogeneous solution between two succesive output points becomes greater than 2*AMP, a new output point is inserted. Set to None to use a sensible default.
• rtol: A relative tolerance for solving the differential equation.
• atol: An absolute tolerance for solving the differential equation.
• output_points: If integer: desired number of output points, based on amplification If array: desired output points If None: let the algorithm decide the output points
• verbosity: 0 silent, 1 some output, 2 more output

Returns:

A tuple (t, y) where t is a (m,) array of mesh points, and y is (m, n) array of solution values at the mesh points.

Raises:

• ValueError – Invalid input
• NoConvergence – Numerical convergence problems
• SingularityError – Infinities occurred
• SystemError – Invalid output from user routines. (FIXME: these should be fixed)

scikits.bvp1lg.mus.solve_nonlinear(func, gsub, initial_guess, a, b, max_amplification=0, rtol=1e-05, atol=None, output_points=None, verbosity=0, iteration_limit=100)

Solve a non-linear two-point boundary value problem.

The problem is assumed to be:

u'(t) = f(t, u(t)),                   a <= t <= b
g_j(u(a), u(b)) = 0                   j = 0, 1, ..., n

where u is the solution vector with shape (n,).

Parameters:

• func: Function f = f(t, u) that evaluates derivatives. It should return an array of shape (n,).

• gsub: Function g, dga, dgb = g(u_a, u_b) that evaluates the boundary condition function and its Jacobian. The output should be:

  g[i] = g_i(u_a, u_b),                          i = 0, ..., n-1
  dga[i,j] = d g_i(u_a, u_b) / d u_a[j],         j = 0, ..., n-1
  dgb[i,j] = d g_i(u_a, u_b) / d u_b[j],
  

• initial_guess: Either callable u = initial_guess(t) that provides an initial guess for the solution vector, or a tuple (t, u) providing mesh and values for the guess — these are interpolated linearly to form the guess at all points. (Previous output from musn can be used as initial_guess.)

• a: Left boundary point.

• b: Right boundary point.

• max_amplification: The allowed incremental factor of the homogeneous solutions between two succesive output points. If the increment of a homogeneous solution between two succesive output points becomes greater than 2*AMP, a new output point is inserted. Set to None to use a sensible default.

• rtol: A relative tolerance for solving the differential equation.

• atol: An absolute tolerance for solving the differential equation.

• output_points: If integer: desired number of output points, based on amplification If array: desired output points If None: let the algorithm decide the output points

• verbosity: 0 silent, 1 some output, 2 more output

• iteration_limit: Maximum allowed number of Newton iterations

Returns:

A tuple (t, u) where t is a (m,) array of mesh points, and u is (m, n) array of solution values at the mesh points.

Raises:

• ValueError – Invalid input
• NoConvergence – Numerical convergence problems
• SingularityError – Infinities occurred
• SystemError – Invalid output from user routines. (FIXME: these should be fixed)