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:

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:

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

Results from running TWPBVP{LC}

Numerical results from running this code can be found here.

The Test Problems

The 35 test problems can be found in the following file.

Other Codes

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

References

  1. S. D. Capper, J. R. Cash and F. Mazzia. On the Development of Effective Algorithms for the Numerical Solution of Singularly Perturbed Two-Point Boundary Value Problems. Submitted.
  2. J. R. Cash and F. Mazzia. A new mesh selection algorithm, based on conditioning, for two-point boundary value codes. J. Comput. Appl. Math., 184(2):362--381, 2005.
  3. J. R. Cash, F. Mazzia. Hybrid Mesh Selection Algorithms Based on Conditioning for Two-Point Boundary Value Problems. Journal of Numerical Analysis, Industrial and Applied Mathematics Vol 1, No 1, pp 81-90, 2006.
  4. J. R. Cash, F. Mazzia, N. Sumarti, and D. Trigiante. The role of conditioning in mesh selection algorithms for first order systems of linear two point boundary value problems. J. Comput. Appl. Math. Vol 185 No 2, pp 212-224, 2006.
  5. D. F. Lawden. Elliptic Functions and Applications. Springer-Verlag Applied Mathematical Sciences 80.
  6. A. E. H. Love. A Treatise on the Mathematical Theory of Elasticity, Fourth Edition, Dover.

Contacting the authors

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

Valid XHTML 1.0 Strict
Valid CSS!
Get Firefox!

This page has been visited 4132 times.