Download Bvp6c: High-Order Software for Solving ODE Boundary Value Problems and more Lab Reports Engineering in PDF only on Docsity!
Report no. 08/
A Sixth-Order Extension to the MATLAB Package
bvp4c of J. Kierzenka and L. Shampine
by NICHOLAS HALE Oxford University Computing Laboratory and DANIEL R. MOORE Dept. of Mathematics, Imperial College London
A new two-point boundary value problem algorithm based upon the MATLAB bvp4c package of Kierzenka and Shampine is described. The al- gorithm, implemented in a new package bvp6c, uses the residual control framework of bvp4c (suitably modified for a more accurate finite difference approximation) to maintain a user specified accuracy. The new package is demonstrated to be as robust as the existing software, but more efficient for most problems, requiring fewer internal mesh points and evaluations to achieve the required accuracy.
Oxford University Computing Laboratory Numerical Analysis Group Wolfson Building Parks Road Oxford, England OX1 3QD E-mail: nick.hale@comlab.oxford.ac.uk April, 2008
2 N. HALE AND D. R. MOORE
1 Introduction
In [KS01] Kierzenka and Shampine describe a software package, bvp4c, to solve a large class of boundary value problems (BVPs) for ordinary differential equations (ODEs) in MATLAB. Specifically, they consider first order ODE systems of the form;
y′^ = f (x, y; p), a ≤ x ≤ b (1.1a)
subject to two-point boundary conditions (BCs)
g(y(a), y(b); p) = 0. (1.1b)
They, as shall we, assume that f (x, y; p) and g(x, y; p) are as smooth as necessary, and in particular that f is continuous and satisfies a Lipschitz condition in y. The argument p is a vector of unknown parameters and is usually suppressed for clarity of exposition. The view of K&S is that in general, a user solving a BVP of the form (1.1) in MATLAB is interested in only a graphical representation of a solution. As such, a modest order solver such as the MIRK4 based Simpson method is appropriate for graphical accuracy. It is our opinion that whilst a fourth-order solver is reasonable, recent developments mean that a sixth-order solver will supply not only greater accuracy, but also perform more efficiently. We thus develop a new piece of software, bvp6c, intended to solve the same class of problems efficiently, whilst maintaining the accuracy and robustness of the original software. Importantly we maintain other valuable attributes, such as treating non-separated BCs and unknown parameters without requiring the user to reformulate these as higher order problems in the canonical form (1.1). Although MIRK (mono-implicit Runge-Kutta) formulae become more complex as order increases [CS82, DC01] - needing more function evaluations, more interior points and more arithmetic operations at each mesh point - we believe a higher order method will give comparable accuracy on a sufficiently coarser mesh (or alternatively, greater accuracy on a comparable mesh) and thus prove more efficient in the majority of cases.
2 Collocation Method
A number of existing methods for solving problems of the form (1.1) provide solutions only at mesh points, whereas others provide data everywhere but without uniform ac- curacy [CS82, CW91a, EM86, EM96]. Our method, like bvp4c, will provide a uniform prescribed accuracy throughout the computational interval. To achieve this we fit the MIRK6 data of [CS82] with the sixth-order interpolant described in [CM04], giving an accurate representation of the solution throughout the computational interval at little extra cost. If hi = xi+1 − xi, the MIRK6 formula given by [CS82] is
yi+1 = yi +
hi 90
[
7 f (xi, yi) + 32f (xi+1/ 4 , yi+1/ 4 ) + 12f (xi+1/ 2 , yi+1/ 2 ) + ...
32 f (xi+3/ 4 , yi+3/ 4 ) + 7f (xi+1, yi+1)
]
= yi +
hi 90
[
7 y i′ + 32y′ i+1/ 4 + 12y i′+1/ 2 + 32y′ i+3/ 4 + 7y i′+
]
4 N. HALE AND D. R. MOORE
3 Residual
As with bvp4c, the cornerstone of bvp6c is that both error estimation and mesh selection are based on the residual of S(x), defined by
r(x) = S′(x) − f (x, S(x)), (3.1) rBCs(x) = g(S(a), S(b))
for the ODEs and for the BCs respectively. Applying (2.7) to the ODE residual gives
r(x) = (S′(x) − f (x, S(x))) − (y′(x) − f (x, y(x))), (3.2) = S′(x) − y′(x) + O(h^6 i ),
for x in each [xi, xi+1]. We now seek to establish the asymptotic behaviour of this residual. Let q(x) be the interpolant (2.3,2.4), but of the true solution y(x). Subtracting from S and differentiating leads to
S′(x) − q′(x) = S′(xi + ωhi) − q′(xi + ωhi), (3.3)
=
hi
[A′ 66 (w) (yi+1 − y(xi+1)) − A′ 66 (1 − w) (yi − y(xi))] +
B′ 66 (w)
y i′+1 − y′(xi+1)
- B 66 ′ (1 − w) (y i′ − y′(xi)) + C 66 ′ (w)
{y′ i+3/ 4 − y′ i+1/ 4 } − {y′(xi+3/ 4 ) − y′(xi+1/ 4 )}
D 66 ′ (w)
y′ i+1/ 2 − y′(xi+1/ 2 )
The errors in y′ i − y′(xi), y′ i+1/ 2 − y′(xi+1/ 2 ) and y i′+1 − y′(xi+1) are O(h^6 i ) [CM04], and
B′ 66 and D 66 ′ are O(1), hence
S′(x) − q′(x) =
hi
[A′ 66 (w) (yi+1 − y(xi+1)) − A′ 66 (1 − w) (yi − y(xi))] (3.4)
C 66 ′ (w)
{y i′+3/ 4 − y i′+1/ 4 } − {y′(xi+3/ 4 ) − y′(xi+1/ 4 )}
The sixth-order MIRK6 formula (2.1) is satisfied by y(x) with an local truncation error O(h^7 i ), i.e. yi+1 − y(xi+1) = yi − y(xi) + O(h^7 i ) [CS82]. Moreover,
{y i′+ 3 4 − y′ i+ 1 4 } − {y′(xi+ 34 ) − y′(xi+ 14 )} =
3 h^5 i 1024
[
∂f ∂y
yv^ −
d dx
∂f ∂y
yiv
)]
=: errC 66 h^5 i + O(h^6 i ). (3.5)
This leads to
S′(x) − q′(x) =
hi
(A′ 66 (w) − A′ 66 (1 − w)) (yi − y(xi)) + C 66 ′ (w)errC 66 h^5 i + O(h^6 i ), (3.6)
but noting A′ 66 (w) − A′ 66 (1 − w) = 0 reduces this to
S′(x) − q′(x) = C 66 ′ (w)errC 66 h^5 i + O(h^6 i ). (3.7)
SIXTH-ORDER EXTENSION TO BVP4C 5
Substituting (3.7) to the expression for the residual in (3.1), we have
r(x) = S′(x) − f (x, S(x)), (3.8) = S′(x) − y′(x) + O(h^6 i ), = q′(x) − y′(x) + C 66 ′ (w)errC 66 h^5 i + O(h^6 i ),
and [CM04] gives the interpolation error of q(x) to y(x) in [xi, xi+1] as
h^6 i 720
w^2 (1 − w)^2
w^2 − w +
yvi
ξ∈[0,1].^ (3.9)
Thus, differentiating (3.9) with respect to x and substituting to (3.8) gives
r(x) =
h^5 3840
w(4w − 1)(2w − 1)(4w − 3)(w − 1)yvi(ξ) + h^5 C′ 66 (w)errC 66 + O(h^6 )
= h^5
[
w(4w − 1)(2w − 1)(4w − 3)(w − 1)yvi(ξ) − ...
16 3
w(2w − 1)(w − 1)errC 66
]
The behaviour of the residual then is local to each subinterval and of order h^6 i at the two end-points and the midpoint. To leading order the residual is a polynomial of degree 5, with dependence on yiv, yv^ and yvi. This extra dependence not present in the bvp4c residual is complex, and in particular leaves us unable to determine local extrema or an asymptotically correct estimate of the L∞ norm. However, in [KS01] the use of the L∞ norm is argued against, and we instead focus on the L 2 norm. On each subinterval we compute
‖r(x)‖i =
(∫ (^) xi+
xi
r^2 (x)dx
and use a suitable quadrature method to integrate the now leading order degree ten polynomial r^2 (x) accurately. An n-point Lobatto quadrature rule is exact for polynomi- als of degree 2n−3, so for an asymptotically correct estimate of our residual we need at least a 7-point procedure, requiring four extra function evaluations per subinterval, five including the evaluation to compute the more accurate y′ i+1/ 2.
4 Implementation
The new software bvp6c is intended to be a direct extension to bvp4c, and as such implementation is almost identical, with support routines changed only where necessary to maintain sixth-order accuracy. We also wish to retain other functionality of the original package, such as solving problems with unknown parameters, generalised two- point BCs, not requiring (although allowing) user entered partial derivatives and the vectorisation of f (x, y). Fortunately the change from fourth- to sixth-order has only a small impact on these areas, and little alteration is necessary. We discuss here some of the larger impacts the higher order extension has on implementation of the method.
SIXTH-ORDER EXTENSION TO BVP4C 7
∂Φi+ ∂p
h 90
[
7 (Ki + Ki+1) + 32
Ki+1/ 4 + Ki+3/ 4
h 2
Ji+1/ 2
5 (Ki+1 − Ki) + 16
Ki+1/ 4 − Ki+3/ 4
3 h 4
Ji+1/ 4 (3Ki − Ki+1) − Ji+3/ 4 (Ki − 3 Ki+1)
3 h 2
Ji+1/ 4 (3Ki − Ki+1) + Ji+3/ 4 (Ki − 3 Ki+1)
}]
where here Kl = ∂f (xl, yl)/∂p. In bvp4c the local internal Jacobian Ji+1/ 2 is approximated by an average when not changing rapidly. Here we have two additional internal Jacobians to calculate (Ji+1/ 4 and Ji+3/ 4 ), making such an approximation even more attractive. In [CS82] further simplifications for these Jacobians are set out that might be used to reduce overall computation time, but we have not applied these in current implementation. Our trials (see §5.3) suggest that although a user supplied Jacobian is not necessary, it can reduce computation speed compared with calculating these Jacobian by numerical differences. This result for the MATLAB package differs from Fortran implementations of the MIRK algorithm for two-point BVPs, where numerical difference calculation of the Jacobians is substantially faster than calculating the analytic matrix multiplication forms for all but the most complicated right hand sides [CGM02].
4.2 Residual Norm
Having found an approximation to the solution by solving the (4.1), we must now com- pute the norm of the residual (3.11) on each subinterval. As discussed in the previous chapter, we choose the 7-point Lobatto quadrature method, the abscissa and weights of which can be found in [Mic63]. This method will integrate a polynomial of degree 10 exactly, and the quadrature points include the end and the midpoints of the interval. Suppose we evaluate the interpolant S and its derivative at xi+1/ 2 ,
S(xi+1/ 2 ) =
(yi+1 + yi) −
hi 24
[
y′ i+1 − y′ i + 4(y′ i+3/ 4 − y′ i+1/ 4 )
]
= yi+1/ 2 , S′(xi+1/ 2 ) = y′ i+1/ 2. (4.6)
We see that this is exact for both y and y′^ (the improved approximations) at the mid- point, and as such
r(xi+1/ 2 ) = S′(xi+1/ 2 )−f (xi+1/ 2 , S(xi+1/ 2 )) (4.7) = y′ i+1/ 2 −y′ i+1/ 2 = 0
needs no calculation. Similarly the residual gives zero contributions at the two end points, but the values of the residual at the four remaining quadrature points must be computed, requiring four evaluations of the right hand side functional f.
8 N. HALE AND D. R. MOORE
In [KS01] a similar analysis is used to show how the residual at the midpoint in bvp4c gives an approximation to how well the collocation equations (2.1) are satisfied. Such a result using values at the quarter-points xi+1/ 4 , xi+3/ 4 is not available here because our interpolant S(x) uses the altered value yi+1/ 2.
4.3 Mesh Selection
The norm of the residual computed as above is then used to update the mesh. We maintain the style of mesh selection used in bvp4c, but note that since bvp6c is a higher order method, introducing/removing points will have a greater effect on the size of the residual. Points are removed if the residual on an subinterval (and an estimate of the new residual on the coarser mesh) is deemed small enough, and added if the residual is too large. The policy in bvp4c is to remove mesh points if the predicted residual is less than half of reltol, the residual tolerance. We find with bvp6c this can cause trouble, with the algorithm repeatedly removing and replacing the same point on alternate iterations. This is solved by only removing mesh points when the predicted residual is less than one tenth of reltol. We do however follow bvp4c in introducing no more than two new mesh points at a time, and the experiments set out below suggest that this slight variation of the mesh refinement strategy proposed by [KS01] is robust enough to ensure user mandated accuracy for almost all test problems.
5 Results & Examples
We demonstrate the accuracy and efficiency of our new method, in particular comparing the results with those computed using the existing software bvp4c. We reconsider some of those problems discussed in [KS01] as well as applying a new set of problems from as the BVP problem test suite [CW91a, CW91b]. For the first family of problems we consider explicit solutions are not available, and we resort to comparing the computed solutions between the two methods. We also compare other important factors, such as solution time, number of mesh points and number of function evaluations. By defining ‘reltol’ and ‘abstol’ to some tolerance, we expect the solution of both methods, and hence their difference, to be approximately of this order.
5.1 Measles
The first example [AMR88, 1.10] models the spread of measles via the differential equa- tions
y′ 1 = μ − β(t)y 1 y 3 (5.1a) y′ 2 = β(t)y 1 y 3 −
y 2 λ
(5.1b)
y′ 3 =
y 2 λ
y 3 ν
(5.1c)
10 N. HALE AND D. R. MOORE
We take R = 100 and obtain the results seen in Figure 1. In Table 5.2 we see that the discrepancy between the two methods behaves as expected, and for the most stringent tolerances the efficiency in the new method is even more pronounced.
bvp4c bvp6c tol time mesh maxres time mesh maxres discrepancy 1e-3 0.401 25 5.8e-4 0.448 24 2.6e-5 3.0e- 1e-6 1.284 152 9.9e-7 0.623 37 8.9e-7 3.0e- 1e-9 7.756 910 9.9e-10 1.449 113 8.2e-10 8.0e- 1e-12 99.37 6878 1.0e-13 3.876 358 9.9e-13 5.2e-
Table 2: As for Table 5.1, for (5.2), the fluid injection problem of [AMR88].
Figure 1: Solutions for (5.2), the fluid injection problem [AMR88, 1.4] solved for a tolerance of 10−^6. The bvp6c solutions for f ′^ and h have been shifted upwards to aid comparison of the final solution meshes for the two algorithms.
We include the graphs of f ′^ and h (note that one set of solutions has been shifted upwards) in Figure 1 to show that the bvp6c solution is, at least graphically, the same as that computed in [KS01], and to demonstrate the distribution of the mesh points for both methods. In this example too we see that bvp6c requires a far sparser mesh to compute the solution to the mandated accuracy. For six digits of accuracy bvp6c uses only 37 points, while bvp4c needs a mesh of 152 points. For nine digits this becomes 113 and 910 respectively, and for twelve digits of accuracy the new algorithm is a little over 25 times faster.
5.3 Shockbvp
The final BVP in this section [AMR88, 9.2] solves the singular problem (as ε → 0),
εy′′^ + xy′^ = −επ^2 cos(πx) − (πx)sin(πx) (5.3a)
with boundary conditions y(−1) = − 2 , y(1) = 0. (5.3b) We use this example to demonstrate a number of key points, particularly the effects of
SIXTH-ORDER EXTENSION TO BVP4C 11
tol=1e-3 tol=1e-6 tol=1e- Final ε 1e-3 1e-4 1e-5 1e-3 1e-4 1e-5 1e-3 1e-4 1e- bvp4c 0.46 0.98 2.79 2.50 4.40 7.45 13.33 24.42 41. bvp6c 0.93 1.57 4.75 0.98 1.73 3.67 2.36 5.41 22. vbvp4c 0.58 1.16 3.21 2.63 4.72 8.05 14.47 26.19 42. vbvp6c 0.49 1.17 3.98 0.96 1.69 3.34 2.13 4.73 19. ajbvp4c 0.23 0.47 1.27 1.14 2.00 3.39 5.84 10.69 17. ajbvp6c 0.24 0.56 1.38 0.54 0.96 2.83 1.44 3.09 13. vjbvp4c 0.18 0.33 0.81 0.76 1.29 2.11 3.59 6.57 10. vjbvp6c 0.16 0.37 0.87 0.32 0.54 1.51 0.77 1.67 7.
Table 3: Time (secs) to solve (5.3), the shock BVP [AMR88, 9.2] using continuation for various MATLAB implementations of the bvpc algorithms.
analytic Jacobians and vectorisation of the problem. Following the steps of the bvp4c analysis in [KS01], the shock problem is solved using continuation and the times shown are not the intermediate results, but the total time taken to reach a solution for the given value of ε. vbvp represents the results from a vectorised version of the problem, ajbvp from using an analytic Jacobian and vjbvp from combining the two. As a general observation on both solvers, it can be seen that vectorisation has a smaller effect on the computational time than that of analytic Jacobians, although it is when the two options are used in tandem that the time is most significantly reduced. For the weakest tolerance setting the new algorithms are the slightly slower of the two, but as the tolerance is made stricter, the improvement in bvp6c is significant. For the smallest value of ε and the strictest tolerance the gain in speed is perhaps not as large as we might expect, and we believe this to be caused by a cycle of adding and removing points unnecessarily on successive iterations as discussed in §4.3. This could perhaps be improved by tweaking the mesh adaption procedure further, but we make no attempt to do so here. This problem demonstrates that if analytic Jacobians are available then both bvp4c and bvp6c significantly improve their performance, which we again note is in contrast to Fortran implementations of MIRK algorithms.
5.4 BVP ODE Suite
The bvp6c package has been tested extensively using the test suite [CW91a, CW91b]; a collection of linear and nonlinear BVPs designed to test the performance and robustness of a numerical BVP solver. Since a number of the problems posed do not have explicit analytic solutions, numerical solutions for these are computed using a recently developed twelfth-order MIRK method [CM06a, CM06b] to a high degree of precision and assumed sufficiently accurate for error analysis. Solutions for each of the 32 problems are then computed using both bvp4c and bvp6c with a range of specified tolerances. For simplicity we choose an initial guess of zero (where possible) on a grid of 33 equally spaced points. We do not make any attempt to vectorise the right hand side, nor do we supply an analytic Jacobian.
SIXTH-ORDER EXTENSION TO BVP4C 13
(^00 8 16 24 )
CW Problem No.
Time (secs)
tol = 10−
(^00 8 16 24 )
1
2
3
CW Problem No.
Time (secs)
tol = 10−
(^00 8 16 24 )
5
10
15
20
25
CW Problem No.
Time (secs)
tol = 10−
(^00 8 16 24 )
5
10
15
20
25
30
35
40
45
50
CW Problem No.
Time (secs)
tol = 10−
Figure 3: Times on the 32 BVP test problems [CW91b] for bvp4c (white & dashed line) and bvp6c (black & solid line) for user tolerances 10−^3 , 10−^6 , 10−^9 & 10−^12. Squares indicate the accuracy is less than the prescribed tolerance for this problem, and triangles show bvp4c repeated with the maximum allowable mesh size doubled.
In those few problems in which bvp6c appears the the slower solver, returning to Figure 2d (and indicated by squares) one sees that these are problems for which bvp4c fails to obtain a solution of sufficient accuracy within the default maximum mesh size (squares). Further tests (triangles in Figure 3d) show that allowing a mesh large enough to reach a solution of the required accuracy with bvp4c takes significantly longer. To achieve the specified tolerance on problems #15 and #16, bvp4c requires 14,648 and 15,919 mesh points respectively, and a computational time of around two minutes.
Further experiments have shown that at strict tolerances there are few examples where the solvers supply a result of the required accuracy, but the residual is not re- duced enough to terminate the algorithm. This suggests the residual is sometimes too overestimated, and a significantly example of this is solving problem #32 with bvp6c. Using default settings, the residual is not reduced below the tolerance before the max- imum number of mesh points (500) is exceeded, which takes around 30 seconds. By reducing this maximum to only 200 points, a solution accurate to 13 digits is obtained in less than one second. This is the result given in Figures 2 and 3.
14 N. HALE AND D. R. MOORE
6 Conclusions
Extensive testing, particularly in §5.4, has shown that whilst not always providing the more accurate of the solutions, bvp6c consistently obtained the accuracy specified by the user. A meaningful comparison of accuracy between the two bvpc packages is difficult to quantify, since the design of the framework means that both solvers continue to improve their solutions until the residual is sufficiently reduced. We have shown that locally the bvp6c residual is asymptotically sixth- and the bvp4c fourth-order accurate, but this type of analysis is relevant only when considering solu- tions obtained on the same mesh. However, since the residual determines the adaptive mesh strategy, we expect and demonstrate that the higher order finite difference for- mula of [CS82] and interpolation formula of [CM04] obtain a comparable reduction in the global residual with fewer internal mesh points, and correspondingly takes less time to compute. We conclude that whilst both methods generally give solutions satisfying the required tolerance, the new method does so more efficiently. This can be further seen in a comparison of solution times for the bvp4c solver with a tolerance of 10−^6 with those from bvp6c at 10−^9. Figure 4 show these are almost identical; meaning that in the same time it takes bvp4c to solve to an accuracy of 10−^6 , bvp6c can solve the same problem to an accuracy of 10−^9. Such a result holds to an even greater extent for tolerances of 10−^9 and 10−^12.
(^00 8 16 24 )
1
2
3
CW Problem No.
Time (secs)
(^100 8 16 24 ) −
10 −
10 −
10 −
10 −
10 −
CW Problem No.
Lerror^2
Figure 4: Solution times and L 2 errors for bvp4c with a tolerance of 10−^6 (white & dashed line) and bvp6c with a tolerance of 10−^9 (black & solid line).
In view of the apparent success of the bvp6c extension to the existing MATLAB software, it is tempting to consider a further development using eighth [CS82], tenth and even twelfth-order [CM06b] difference methods. However, the eighth-order scheme requires a total of seven interior points per interval, and the form of the Jacobians in §4.1 becomes increasingly complex. It seems likely the improvement in local accuracy of an even higher order scheme may come at the expense of an unjustifiable increase in the time needed to find a global solution of the required accuracy. An area where we do see potential for the enhancement of the bvpc framework is in solving more efficiently problems of the form y′′(x) = f (x, y) and y′′^ = f (x, y, y′). Again
16 N. HALE AND D. R. MOORE
All results in this paper were computed using MATLAB 7.5 on a dual core 1.7Ghz laptop. The parameters used for the 32 problems in §5.4 are 0.001, 0.01, 0.05, 0.025, 0.01, 0.022, 0.025, 0.01, 0.055, 0.022, 0.001, 0.0025, 0.0025, 0.0025, 0.005, 1/19, 0.0005, 0.013, 0.03, 0.05, 0.0008, 0.025, 5, 0.03, 0.0025, 0.02, 0.02, 0.03, 0.015, 0.042, 0.025 and 100 respectively.
The bvp6c package and the examples of §5 are freely available online on both the MATLAB Central file exchange at http://www.mathworks.com, and on NH’s website http://www.comlab.ox.ac.uk/nick.hale/bvp6c.
Finally, we wish to thank Steven Capper for pointing out the more economical form for the Jacobians in the [CS82] MIRK6 scheme, and NH wishes to thank the Nuffield Foundation for an Undergraduate Research Bursary for 2005.
References
[AMR88] U. M. Ascher, R. M. R. Mattheij, and R. D. Russell. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Prentice–Hall, Englewood Cliffs, NJ, USA, 1988.
[CCM06] J. R. Cash, S. D. Capper, and D. R. Moore. Lobatto-Obrechkoff formulae for 2nd order two-point boundary value problems. Journal of Numerical Analysis, Industrial and Applied Mathematics, 1:13–25, 2006.
[CGM02] J. R. Cash, M. P. Garcia, and D. R. Moore. Mono-implicit runge-kutta for- mulae for the numerical solution of second order nonlinear two-point bound- ary value problems. Journal of Computational and Applied Mathematics, 143(2):275–289, 2002.
[CM04] J. R. Cash and D. R. Moore. High-order interpolants for solutions of two-point boundary value problems using MIRK methods. Computers and Mathematics with Applications, 48(10–11):1749–1763, 2004.
[CM06a] S. D. Capper and D. R. Moore. NewNRK a Fortran 90 code for solving two-point boundary value problems, 2006.
[CM06b] S. D. Capper and D. R. Moore. On high order MIRK schemes and Hermite- Birkhoff interpolants. Journal of Numerical Analysis, Industrial and Applied Mathematics, 1:27–47, 2006.
[CS82] J. R. Cash and A. Singhal. High order methods for the numerical solution of two-point boundary value problems. Behaviour & Information Technology, 22:184–199, 1982.
SIXTH-ORDER EXTENSION TO BVP4C 17
[CW91a] J. R. Cash and M. H. Wright. A deferred correction method for nonlinear two- point boundary value problems: implementation and numerical evaluation. SIAM Journal of Scientific and Statistical Computing, 12(4):971–989, 1991.
[CW91b] J. R. Cash and M. H. Wright. Thirty–two test problems for analysis of two- point boundary value problem solvers, 1991. http://www.ma.ic.ac.uk/˜jcash/BVP software/PROBLEMS.PDF.
[DC01] M. Van Daele and J. R. Cash. Superconvergent deferred correction methods for first order systems of nonlinear two–point boundary value problems. SIAM Journal of Scientific Computing, 22(5):1697–1716, 2001.
[EM86] W. H. Enright and P. H. Muir. Efficient classes of Runge–Kutta methods for two–point boundary value problems. Computing, 37:315–334, 1986.
[EM96] W. H. Enright and P. H. Muir. Runge–Kutta software with defect control for boundary value ODEs. SISSC, 2:479–497, 1996.
[KS01] J. Kierzenka and L. F. Shampine. A BVP solver based on residual control and the MATLAB pse. ACM Transactions on Mathematical Software, 27(3):299– 316, 2001.
[KS08] J. Kierzenka and L. F. Shampine. A BVP solver that controls residual and error. http://faculty.smu.edu/shampine/finalbvp5c.pdf, 2008.
[Mic63] H. H. Michels. Abscissas and weight coefficients for Lobatto quadrature. Math- ematics of Computation, 17(83):237–244, 1963.
