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.
A simple test problem we consider is Troesch's equation, given by
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:
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 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.
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:
With the following boundary conditions:
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.
In order to use the software on this page, one must implement the following:
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
All the download links on this page can be found summarised in this section:
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. |
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.
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
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
The driver code for TWPBVPL can be tested by running:
./twpbvpbl
The driver code for TWPBVPC can be tested by running:
./twpbvpbc
Numerical results from running this code can be found here.
The 35 test problems can be found in the following file.
The codes on this website are by no means the only codes with which to solve BVPs.
All comments (both good and bad) should be addressed to:
This page has been visited 4284 times.