# 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

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:

1. y(0) = 0,
2. 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:

f77 -o troeschs twpbvpc.f troeschs.f

the compiled code can be executed by typing:

./troeschs

and one should get output similar to:

```PROGRAM TROESCHS

Number of mesh points =  12
Dumping (x,y) mesh now...
0.0000000     0.0000000
0.20000000    0.10753422E-01
0.40000000    0.33200511E-01
0.53333333    0.65658365E-01
0.66666667    0.12915421
0.73333333    0.18187245
0.80000000    0.25821664
0.86666667    0.37332998
0.90000000    0.45506034
0.93333333    0.56435971
0.96666667    0.72289433
1.0000000     1.0000000
```

## 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:

1. x(0) = 0
2. y(0) = 0
3. κ(0) = 0
4. y(0.5) = 0
5. φ(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`.

The differential system coded up in a suitable form is as follows:

```SUBROUTINE FSUB(NCOMP,X,Z,F,RPAR,IPAR)
IMPLICIT NONE
INTEGER NCOMP, IPAR
DOUBLE PRECISION F, Z, RPAR, X
DIMENSION Z(*),F(*)
DIMENSION RPAR(*), IPAR(*)

F(1)=cos(Z(3))
F(2)=sin(Z(3))
F(3)=Z(4)
F(4)=Z(5)*cos(Z(3))
F(5)=0

RETURN
END
```

The analytic Jacobian for the F-function:

```      SUBROUTINE DFSUB(NCOMP,X,Z,DF,RPAR,IPAR)
IMPLICIT NONE
INTEGER NCOMP, IPAR, I, J
DOUBLE PRECISION X, Z, DF, RPAR
DIMENSION Z(*),DF(NCOMP,*)
DIMENSION RPAR(*), IPAR(*)

DO I=1,5
DO J=1,5
DF(I,J)=0.D0
END DO
END DO

C     dF1/dZ3
DF(1,3)=-sin(Z(3))
C     dF2/dZ3
DF(2,3)=cos(Z(3))
C     dF3/dZ4
DF(3,4)=1.0D0
C     dF4/dZ3
DF(4,3)=-Z(5)*sin(Z(3))
C     dF4/dZ4
DF(4,4)=1.0D0
C     dF4/dZ5
DF(4,5)=cos(Z(3))

RETURN
END
```

The boundary conditions can be coded up in Fortran 77:

```      SUBROUTINE GSUB(I,NCOMP,Z,G,RPAR,IPAR)
IMPLICIT NONE
INTEGER I, NCOMP, IPAR
DOUBLE PRECISION Z, RPAR, G
DIMENSION Z(*)
DIMENSION RPAR(*), IPAR(*)

C     I == the boundary condition "number".
C     The conditions at the left are enumerated first,
C     then the ones at the right.
C     The order of left or right conditions does not matter,
C     but we must be consistent when defining the jacobian
C     of the boundary conditions!

C     We have specified 3 left bcs, and 5 bcs total.
C     This means that:
C     BC(1) = x(0) = 0
C     BC(2) = y(0) = 0
C     BC(3) = kappa(0) = 0
C     BC(4) = y(0.5) = 0
C     BC(5) = phi(0.5) = -pi/2

IF (I.EQ.1) G=Z(1)
IF (I.EQ.2) G=Z(2)
IF (I.EQ.3) G=Z(4)
IF (I.EQ.4) G=Z(2)
IF (I.EQ.5) G=Z(3)+1.5707963267948966192313216916397514D0

RETURN
END
```

The analytic Jacobian for the G-function:

```      SUBROUTINE DGSUB(I,NCOMP,Z,DG,RPAR,IPAR)
IMPLICIT NONE
INTEGER I, NCOMP, IPAR
DOUBLE PRECISION Z, DG, RPAR
DIMENSION Z(*),DG(*)
DIMENSION RPAR(*), IPAR(*)

C     I == the boundary condition "number".
C     The conditions at the left are enumerated first,
C     then the ones at the right.
C     The order of left or right conditions does not matter,
C     but we must be consistent when defining the jacobian
C     of the boundary conditions!

C     We have specified 3 left bcs, and 5 bcs total.
C     This means that:
C     BC(1) = x(0) = 0
C     BC(2) = y(0) = 0
C     BC(3) = kappa(0) = 0
C     BC(4) = y(0.5) = 0
C     BC(5) = phi(0.5) = -pi/2

DG(1)=0.D0
DG(2)=0.D0
DG(3)=0.D0
DG(4)=0.D0
DG(5)=0.D0

C     dG1/dZ1
IF (I.EQ.1) DG(1)=1.D0
C     dG2/dZ2
IF (I.EQ.2) DG(2)=1.D0
C     dG3/dZ4
IF (I.EQ.3) DG(4)=1.D0
C     dG4/dZ2
IF (I.EQ.4) DG(2)=1.D0
C     dG5/dZ3
IF (I.EQ.5) DG(3)=1.D0

RETURN
END
```

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):

f77 -o elastica twpbvpc.f elastica.f

the code can then be executed:

./elastica

One should get output similar to:

```PROGRAM ELASTICA

Number of mesh points =  16
P =  -21.5490875
Theta0 =   0.71052198
Dumping (x,y) mesh now...
0.000000   0.00000
0.025333  0.216643E-01
0.051056  0.428629E-01
0.077538  0.631017E-01
0.105104  0.818315E-01
0.134003  0.984243E-01
0.164358  0.112158
0.196113  0.122216
0.228957  0.127712
0.262247  0.127755
0.294945  0.121568
0.325599  0.108652
0.352419  0.889974E-01
0.373463  0.632690E-01
0.386946  0.328976E-01
0.391594 -0.321401E-24
```

## Downloads

All the download links on this page can be found summarised in this section:

• Code for the paper On the Development of Effective Algorithms for the Numerical Solution of Singularly Perturbed Two-Point Boundary Value Problems .
• twpbvp.pdf, the manual to TWPBVP.
• twpbvp.tex, the source for the manual to TWPBVP.
• twpbvp.f, the source code for the original TWPBVP.
• twpbvpc.f, the source code for TWPBVPC.
• twpbvplc.f, the source code for TWPBVPLC.
• troeschs.f, Troesch's equation sample (requires TWPBVPC).
• elastica.f, the elastica sample (requires TWPBVPC).
• 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. The source code of TWPBVPC (TWPBVP with Conditioning). The source code of TWPBVPL (TWPBVP using a Lobatto scheme and Conditioning). Driver code for TWPBVPC on 35 test problems. 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:

make

### Compiling manually

The code can be built manually, using the following commands.

f77 -o dtwpbvpbc twpbvpc.f dtwpbvpbc.f

and

f77 -o dtwpbvpbl twpbvpl.f dtwpbvpbl.f

## Executing the Batch Test Code

The driver code for TWPBVPL can be tested by running:

./twpbvpbl

The driver code for TWPBVPC can be tested by running:

./twpbvpbc

## 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.
• COLNEW by Ascher and Bader, this is a modification to COLSYS.
• COLDAE by Ascher and Spiteri, another modification to COLSYS which deals with differential algebraic equations.
• ACDC by Cash and R.W.Wright, a code for solving extremely stiff problems. This code uses continuation applied to Lobatto schemes in a deferred correction framework.
• COLMOD by R.W.Wright and Cash, a modification of COLNEW, for solving very stiff problems.
• MIRKDC by Enright and Muir, a code which makes use of defect correction for error control.

## Contacting the authors

All comments (both good and bad) should be addressed to:

This page has been visited 3475 times.