Algorithms for the Solution of Two-Point Boundary Value Problems

The Fortran 77 code TWPBVP was originally developed by Jeff Cash and Margaret Wright
and is a global method to compute the numerical solution of two point boundary value problems (either linear or non-linear) with
separated boundary conditions. In the code TWPBVP, MIRK schemes of orders 4, 6 and 8 are solved in a deferred correction
framework in an attempt to give a solution accurate to a prescribed local tolerance at a discrete set of mesh points.
TWPBVP has been found to be especially effective at solving non-stiff and mildly-stiff ODEs efficiently.

For problems of a greater stiffness, fully implicit Runge Kutta schemes have been found to be more
suitable (although requiring more work per step) than MIRK schemes.
A deferred correction code TWPBVPL based on Lobatto IIIA schemes of orders 4, 6 and 8 has been developed. For stiff problems
it is recommended that this code should be tried first.

As a solution is sought only at a discrete set of points, this raises the question of how
to choose where these mesh points lie. The codes TWPBVP and TWPBVPL rely solely on the estimated error at each mesh point
in order to choose where the points lie
and they attempt to minimise the number of mesh points required to constitute a valid solution within the
user supplied tolerance. However recent work suggests that the
conditioning parameters of the problem should also be taken into account when selecting a mesh.

Conditioning parameters can be thought of as a measure of the effect that perturbations in both the boundary
conditions and derivative function will have on the final solution.
In this new generation of deferred correction codes the condition numbers are estimated both for the
true analytic solution and the discrete numerical solution. A modified code called TWPBVPC,
attempts to match both sets of condition numbers as part of a hybrid mesh selection strategy. Initial
results from this code are excellent, and the authors suggest that TWPBVPC should
be used instead of TWPBVP (conditioning can be enabled or disabled in TWPBVPC,
by just changing one line of code, so one can revert back to the old mesh selection code
in TWPBVP if required. The
line of code that must be changed to disable the conditioning is clearly marked in the code.).

A similar hybrid mesh selection algorithm has been added to TWPBVPL and this gives rise to a new code
TWPBVPLC. For non-stiff and mildly stiff problems the user is encouraged to try
TWPBVPC, for stiff problems the code TWPBVPLC is recommended and for very stiff problems
the codes COLMOD or ACDC are recommended. If you do not know whether or not your problem is stiff
then try TWPBVPC. (The calling conventions for both TWPBVPC and TWPBVPL are very similar
as can be evidenced from examining the driver code on this site.)

That is the background to the codes. We now give a straightforward example that can easily be modified by the
user to solve different second order problems.

Quick Start: Troesch's Equation.

A simple test problem we consider is Troesch's equation, given by

d²y/dx²

= μsinh(μy).

Increasing μ increases the stiffness of the ODE. For μ=10, the problem is relatively easy to solve,
while for values such as μ=40 the problem becomes very hard to solve. We formulate this problem as a system
of two first order ODES:

Y1 = y, Y2 = y',

where:

Y1' = Y2, Y2' = μsinh(μY1),
μ = 5.

with the boundary conditions:

y(0) = 0,

y(1) = 1.

A heavily commented example Fortran 77 file, troeschs.f, contains
all the code required to implement the above problem. To run this code, one must also download
the code twpbvpc.f. The example can be compiled using the following command:

The Full Documentation to TWPBVP

The full documentation to the original code TWPBVP can be found in the files
twpbvp.tex and twpbvp.pdf. This does not explain any of the
conditioning logic found in the newer codes, but it does go in to great
detail, and will be useful for anyone solving harder BVPs.

A more detailed example: a coincident endpoint elastica

We now wish to give a more complicated problem which is a system of 5 differential equations.
We consider an elastica in the (x, y) plane, with natural parameterisation s∈[0,0.5]. We
track an angle φ rather than separate derivatives for the x and y components as this gives
us control over the curvature κ. The curve can be defined as the solution the following 5 dimensional
first order system:

dx/ds

= cos(φ),

dy/ds

= sin(φ),

dφ/ds

= κ,

dκ/ds

= Fcos(φ),

dF/ds

= 0

With the following boundary conditions:

x(0) = 0

y(0) = 0

κ(0) = 0

y(0.5) = 0

φ(0.5) = -π/2

An expression for the solution to this problem can be formulated in terms of Jacobi elliptic functions and the
complete and incomplete elliptic integrals of the first and second kinds. The reader is referred to
Lawden for more details.

We make the following variable assignments:

Y1 = x, Y2 = y, Y3 = φ,
Y4 = κ, Y5 = F,

and arrive at the following differential system:

Y1' = cos(Y3), Y2' = sin(Y3), Y3' = Y4,
Y4' = Y5cos(Y3), Y5' = 0.

Coding up the Elastica

In order to use the software on this page, one must implement the following:

Code to compute the derivative (a F-function).

An analytic Jacobian for the F-function.

An expression for the boundary conditions (a G-function).

The analytic Jacobian for the G-function.

A driver function to allocate the necessary memory and call TWPBVPC
or TWPBVPL.

Here is an example driver function which allocates memory and calls TWPBVPC:

SUBROUTINE RUNPROB(ETOL)
IMPLICIT NONE
INTEGER LWRKFL, LWRKIN, NUDIM, NXXDIM, LISERIES, IPRINT, ind,
+ IWRK, NTOL, NLBC, NMSH, NFXPNT, NCOMP, LTOL, NMAX, IPAR,
+ IFLBVP, NMINIT, IDUM
DOUBLE PRECISION TOL, ETOL, WRK, UVAL0, ALEFT, ARIGHT,
+ FIXPNT, XX, U, CKAPPA1, GAMMA1, CKAPPA, RPAR
PARAMETER(LWRKFL=250000)
PARAMETER(LWRKIN=50000)
DIMENSION WRK(LWRKFL),IWRK(LWRKIN)
PARAMETER(NUDIM=5,NXXDIM=100000)
DIMENSION U(NUDIM,NXXDIM), XX(NXXDIM)
DIMENSION TOL(NUDIM), LTOL(NUDIM), FIXPNT(1)
DIMENSION RPAR(1), IPAR(1)
PARAMETER (LISERIES=10)
LOGICAL LINEAR, GIVEU, GIVMSH
EXTERNAL TWPBVPC, FSUB, DFSUB, GSUB, DGSUB
LOGICAL PDEBUG, use_c, comp_c
c this common need to be defined in order to run TWPBVPC successfully
c give information about some parameters
COMMON /ALGPRS/ NMINIT,PDEBUG,IPRINT,IDUM,UVAL0,use_c,comp_c
C RPAR (real) and IPAR (integer) are arrays that are passed to
C the F function, G function, and their Jacobians.
C RPAR is typically used as a parameter in the problem formulation, and
C IPAR is normally used to select a problem from a set.
C
C We have no parameters in our problem, so IPAR and RPAR are left alone.
C If you do not like to use the conditioning in the mesh selection
C USE_C = .false.
C If you do not like to compute the conditioning parameters
C COMP_C = .false.
USE_C = .true.
COMP_C = .true.
C WRK is the floating point workspace and IWRK
C is the integer work space. We set these to zero here,
C to ensure consistent behaviour with different compilers.
DO ind=1,LWRKFL
WRK(ind) = 0d0
ENDDO
DO ind=1,LWRKIN
IWRK(ind) = 0
ENDDO
C ALEFT and ARIGHT are the values of x at the left
C and right boundaries.
ALEFT = 0.D0
ARIGHT = 5.D-1
C ETOL is the required error tolerance of the solution.
C Decreasing ETOL will give a more accurate solution, but
C more mesh points are likely to be required and the code
C will take longer to finish.
C ETOL is passed to this subroutine as an argument.
C NTOL is the number of components that have to satisfy
C the prescribed tolerance.
C LTOL is the component that needs to be tested.
C Most of the time one will set NTOL to the number of system components,
C LTOL(i) to component i, and set the tolerance for each component TOL to be equal.
NTOL = 5
TOL(1) = ETOL
DO ind=1,ntol
LTOL(ind)=ind
TOL(ind) =TOL(1)
ENDDO
C IPRINT controls whether or not the solver prints out
C verbose output. A value of -1 disables these diagnostics
IPRINT = -1
C The number of components in the system
NCOMP = 5
C The number of boundary conditions at the left, the
C number of the right being equal to NCOMP-NLBC
NLBC=3
C NMSH is the number of initial points we supply to
C the solver in the form of a guess, we set this
C to zero to indicate that we have no initial guess
NMSH = 0
C fixed points are x values which are required by the
C user to be part of the returned mesh.
C NFXPNT is the number of fixed points, which we set to be zero
NFXPNT= 0
c The initial approximation is equal to UVAL0
UVAL0 = 0.D0
C the problem is nonlinear so we specify .false. here
LINEAR = .FALSE.
C we do not supply any initial guesses, so again we
C choose .false.
GIVEU = .FALSE.
C No initial mesh (a set of x values) are given, so
C this is .false. too
GIVMSH = .FALSE.
C PDEBUG controls the low level solver diagnostics,
C most of the time this should be set to .false.
PDEBUG = .FALSE.
CALL TWPBVPC(NCOMP,NLBC,ALEFT,ARIGHT,NFXPNT,FIXPNT,NTOL,
+ LTOL,TOL,LINEAR,GIVMSH,GIVEU,NMSH,NXXDIM,
+ XX,NUDIM,U,NMAX,LWRKFL,WRK,LWRKIN,IWRK,
+ FSUB,DFSUB,GSUB,DGSUB,
+ ckappa1,gamma1,ckappa,rpar,ipar,IFLBVP)
C When returning from TWPBVPC, one should immediately
C check IFLBVP, to see whether or not the solver was
C succesful
IF (IFLBVP .LT. 0) THEN
WRITE(6,*) 'The code failed to converge!'
RETURN
END IF
C the solution x values are stored in XX, the Y
C values are stored in U.
C U(i,j) refers to component i of point j in the mesh.
WRITE(6,*) 'Number of mesh points = ', NMSH
WRITE(6,*) 'P = ', U(5,1)
WRITE(6,*) 'Theta0 = ', U(3,1)
WRITE(6,*) 'Dumping (x,y) mesh now...'
DO IND=1,NMSH
WRITE(6,*) U(1,IND), U(2,IND)
END DO
RETURN
END

The complete elastica code elastica.f, is available for download. The code can be
compiled with the following command (one also needs to download twpbvpc.f):

twpbvplc041006.tar.gz, a gzip-compressed tarball, containing
TWPBVPC, TWPBVPL and the driver routines for carrying out tests on a batch of problems.

To unpack the code, use the command:

tar -xvzf twpbvplc041006.tar.gz

The following files are contained in the archive:

Makefile

UNIX Makefile used to build the code.

twpbvpc.f

The source code of TWPBVPC (TWPBVP with Conditioning).

twpbvpl.f

The source code of TWPBVPL (TWPBVP using a Lobatto scheme and Conditioning).

dtwpbvpbc.f

Driver code for TWPBVPC on 35 test problems.

twpbvpbl.f

Driver code for TWPBVPL on 35 test problems.

Comparisons between Codes

In the previous examples we have shown how to run standard two point boundary value problems using our codes.
There are now two extensions we wish to make. The first is to allow you to test your code against ours and
the second is to show for what parameter ranges our codes are appropiate. In order to allow you to make a comparison with
our codes we define a set of challenging problems in the following way.

Singularly perturbed boundary value problems are characterised by the presence
of a small, positive parameter which multiplies the highest derivative
term in the (system of) differential equations. The essential difficulties
with such problems arise because their solutions exhibit narrow regions
of very fast variation (for example boundary layers) as the problem parameter
becomes progressively small. For numerical methods, this means that it
is often very difficult to determine both a mesh that will yield a sufficiently
accurate numerical solution, and an initial guess to the solution which
will lead to the convergence of Newton's method.

We provide a batch of test problems which we recommend as a basis for comparison of our codes with other
two point boundary value problem solvers. All the problems are singular perturbation problems with
parameters that can be varied to control the problem stiffness.

Using a Makefile

A Makefile is provided, that will compile both the BVP solving code and the driver routines. A
Fortran 77 compiler named f77 is assumed to exist in a valid path.
For other compilers modifications will need to be made to the FORTRAN variable
in the Makefile. All the code can be compiled by issueing the command:

Results from running TWPBVP{LC}

The Test Problems

Other Codes

The codes on this website are by no means the only codes with which to solve BVPs.

The Matlab code TOM can be found on Francesca Mazzia's
homepage. This code makes use
of the Top Order Method of order 2 and 6. Conditioning is used for the mesh selection.

The Fortran 95 code NewNRK, is currently being developed with the goal of being a replacement
to the Fortran 77 codes. It is not yet finished, but at the moment non-separated BCs, parameters,
solution interpolation, integral constraints, finite difference schemes of orders up to 12 and
specialised finite difference schemes for the solution to second order systems are supported.

A number of high quality codes can be found on netlib.
In particular, one may find the following codes useful:

COLSYS by Ascher, Christiansen, and Russell. A high order
collocation code, which computes the continuous solution to a BVP.