`bvp4cAD`- An Automatic Differentiation Enabled Boundary Value Solver

Shaun Forth
Applied Mathematics & Scientific Computation
Cranfield University (Defence Academy)
Shrivenham, Swindon SN6 8LA, UK
S.A.Forth@cranfield.ac.uk

and

Lawrence F. Shampine
Mathematics Department
Southern Methodist University, Dallas, TX75275, USA

Last Updated 30 August 2005

If you have any comments, suggestions then please

email me

Introduction

In a recent paper [SKF05] we investigated use of the forward mode of automatic differentiation (AD), as coded in the MAD package [For05], to improve the reliability and robustness of the quasi-Newton solve within MATLAB's boundary value solver bvp4c [SKR,KS01]. The resulting AD-enabled boundary value solver was called bvp4cAD and this web-page/document is devoted to its use.

The packages bvp4c/bvp4cAD are designed to solver boundary value problems of the form,

$\begin{displaymath} \frac{d {\bf y}}{dx}={\bf f}(x,{\bf y},{\bf p}) \mbox{ for }... ...b \mbox{ such that } 0={\bf g}({\bf y}(a),{\bf y}(b),{\bf p}), \end{displaymath}$

(1)

where the optional ${\bf p}$ is a vector of unknown parameters. The functions ${\bf f}(x,{\bf y},{\bf p})$ and ${\bf g}({\bf y}(a),{\bf y}(b),{\bf p})$ prescribe the differential equation and its boundary conditions respectively. The user is referred to the MATLAB documentation and [SGT03, Chap. 3],[SKR,KS01] for more details.

Getting and Installing the Software

To use bvp4cAD you need:

MATLAB

: We wrote bvp4cAD using MATLAB 6.5 and have used it under versions 6.1 and 7.0.

MAD

: MAD [For05] is a licensed software package written by Shaun Forth. Currently it is distributed, with a 21 free days trial, from the TOMLAB company www.tomlab.biz with a user guide [FE04] also available. After the 21 days a license fee is payable (please contact Shaun Forth if you wish to continue using MAD and this causes difficulty for you).

bvp4cAD

: The bvp4cAD package (a modified version of the MathWorks copyright bvp4c) is supplied, by kind permission of The Mathworks, together with example programs used in the paper [SKF05] (themselves adapted from those of [SGT03]) as a tarred gzipped file bvp4cAD.tar.gz which may be opened using

    gunzip bvp4cAD.tar.gz
    tar xvf bvp4cAD.tar

on a UNIX/LINUX machine or an application such as WinZip on a PC.

For the rest of this document we assume you are working in the bvp4cAD folder/directory of the distribution, or that you have added this folder/directory to you MATLAB path.

Using `bvp4cAD`

The main difference between bvp4cAD and the standard MATLAB bvp4c is that, in the absence of appropriate user defined functions, the Jacobians ${\bf f}$ and ${\bf g}$ of (1) are calculated by the forward mode of automatic differentiation [Gri00] via the fmad class of the MAD package instead of being approximated by finite-differences (FD). A subtlety of this is that when AD is used Jacobians of ${\bf f}$ with respect to ${\bf y}$ and ${\bf p}$ at both mesh, and inter-mesh points are calculated exactly, whereas when FD is used Jacobians at inter-mesh points may be approximate by averaging those at neighbouring mesh points. Otherwise everything else is identical, with the obvious exception of the function name.

Example 3.5.1

Consider Example 3.5.1 of [SGT03], the Bratu problem involving solution of

$\begin{displaymath}y''+e^y=0,\mbox{ such that } y(0)=y(1)=0.\end{displaymath}$

The code for standard solution of this using bvp4cAD is given as ch3ex1.m in the bvp4cAD distribution and is shown in Figure 1.

**Figure 1:** `cb3ex1.m` : standard (finite-difference) solution of the Bratu problem
$\begin{figure}\begin{verbatim}function sol = ch3ex1(vect) % FUNCTION: ch3ex1: ... ... ];function res = bcs(ya,yb) res = [ ya(1); yb(1) ];\end{verbatim}\end{figure}$

The problem is solved by running as,

ch3ex1('on')

to take advantage of vectorization of the coding of ${\bf f}$ , or,

ch3ex1('off')

to disable vectorization.

To use AD we simply call bvp4cAD as coded in ch3ex1AD.m and shown in Figure 2.

**Figure 2:** `cb3ex1AD.m` : AD solution of the Bratu problem
$\begin{figure}\begin{verbatim}function sol = ch3ex1AD(vect) % FUNCTION: ch3ex1... ... ];function res = bcs(ya,yb) res = [ ya(1); yb(1) ];\end{verbatim}\end{figure}$

Results

The function timeall runs all eight problems of Examples 3.5.1-3.5.8 of [SGT03], and produces a plot and table of the CPU times. CPU times running under MATLAB 6.5 are given in Table 1 and Figure 3.

	unvectorized			vectorized
Problem	`bvp4c`	`bvp4cAD`	Ratio	`bvp4c`	`bvp4cAD`	Ratio
1	0.17	0.14	0.81	0.12	0.12	1.03
2	0.27	0.34	1.27	0.25	0.12	0.50
3	0.47	1.32	2.80	0.42	0.27	0.64
4	1.47	2.13	1.45	1.29	1.01	0.78
5	0.97	3.17	3.27	0.71	0.56	0.79
6	2.26	6.39	2.83	1.64	0.93	0.57
7	1.11	7.36	6.61	0.70	0.73	1.05
8	0.21	0.41	1.98	0.18	0.24	1.31

**Figure 3:** CPU times (s) for Examples 3.5.1-3.5.8 of [SGT03] running under MATLAB 6.5
$\resizebox{10cm}{!}{\includegraphics{bvp65.eps}}$

Also in Table 1 the ratio of AD to FD times is given: a ratio less than one indicates that a speedup has been obtained by using AD. Without vectorization it is only for problem 1 that AD produces a speedup compared to FD and in some cases, particularly problem 7, AD significantly slows performance. With vectorization on AD produces a speedup in five of the 8 cases, equivalent times for two, and a 30% reduction in efficiency just for problem 8. We note that the improvement arising from vectorization is much greater when using AD than when using FD. To understand why vectorisation is so useful for AD we must appreciate that for each arithmetic operation coded in the ODE function only one AD operation, vectorized across all grid points, is performed. In contrast, for the unvectorized case the AD operation is repeated for every grid point. Consequently the overheads of using AD are reduced by a factor approximately equal to the number of grid points in the vectorized as opposed to the unvectorized case.

CPU times running under MATLAB 7.0 are given in Table 2 and Figure 4.

	unvectorized			vectorized
Problem	`bvp4c`	`bvp4cAD`	Ratio	`bvp4c`	`bvp4cAD`	Ratio
1	0.09	0.15	1.77	0.07	0.17	2.54
2	0.34	0.79	2.31	0.55	0.63	1.14
3	0.71	2.63	3.72	1.00	1.36	1.37
4	2.15	4.99	2.32	2.42	2.97	1.23
5	1.60	6.00	3.74	1.74	2.85	1.64
6	3.35	9.65	2.88	3.72	2.68	0.72
7	1.45	11.65	8.03	1.43	2.83	1.98
8	0.44	1.24	2.81	0.42	1.04	2.51

**Figure 4:** CPU times (s) for Examples 3.5.1-3.5.8 of [SGT03] running under MATLAB 7.0
$\resizebox{10cm}{!}{\includegraphics{bvp7.eps}}$

Here use of MATLAB 7 has increased run times on all problems but one! The performance of bvp4cAD has also been slowed by the switch to MATLAB 7. It is also interesting to note that, unlike when using MATLAB 6.5, vectorization no longer improves performance of the FD approach of bvp4c. Regarding use of AD, under MATLAB 7 bvp4cAD only outperforms bvp4c on one vectorized test problem. MATLAB 7 features a more aggressive just-in-time (JIT) compiler compared to MATLAB 6.5. However the JIT compiler does not work in the presence of overloaded arithmetic, as used by MAD, and we attribute the poorer relative performance of the AD-enabled solver to this.

The Fluid Injection Problem

In [SKF05] we considered in detail Problem 5, the Fluid Injection Problem, with solutions calculated with increasing Reynolds number by continuation. The file ch3ex5continuation.m can be used to compare solution robustness and efficiency with and without AD. If run as,

ch3ex5continuation('no')

then a solution, without continuation, is attempted at Reynolds numbers ranging from 100 to 600 in steps of 50 using both bvp4c and bvp4cAD. The calculation shows that bvp4c first fails at Reynolds number 250, whereas bvp4cAD is successful up to a Reynolds number of 350.

By running ch3ex5continuation as,

ch3ex5continuation('yes')

continuation is used for the fluid injection problem and timings, averaged over 10 repeats, of vectorized and unvectorized performance of both bvp4c and bvp4cAD is performed. The number of repetitions can be varied, e.g.,

ch3ex5continuation('yes',N)

performs N repeats.

In Figures 5 and 6

**Figure 5:** The Fluid Injection Problem - unvectorized and running under MATLAB 6.5
$\resizebox{10cm}{!}{\includegraphics{fipnovect65.eps}}$

**Figure 6:** The Fluid Injection Problem - vectorized and running under MATLAB 6.5
$\resizebox{10cm}{!}{\includegraphics{fipvect65.eps}}$

we show CPU times for increasing Reynolds number required to solve this problem under MATLAB 6.5. In the unvectorized case we see that bvp4c clearly outperforms bvp4cAD, whereas in the vectorized case the converse is true with the run times obtained with AD even approaching those from analytic partial derivatives (hand-coding).

In Figures 7 and 8 we repeat our solution of the fluid injection problem but now using MATLAB 7.0.

**Figure 7:** The Fluid Injection Problem - unvectorized and running under MATLAB 7
$\resizebox{10cm}{!}{\includegraphics{fipnovect7.eps}}$

**Figure 8:** The Fluid Injection Problem - vectorized and running under MATLAB 7.0
$\resizebox{10cm}{!}{\includegraphics{fipvect7.eps}}$

The chief difference is that we now need to get to high Reynolds number, and an associated high number of mesh points, before the vectorized bvp4cAD outperforms bvp4c.

Conclusions

This short report shows, with examples, how BVPs can be solved using the AD-enabled solver bvp4cAD. When running several vectorized test problems under MATLAB 6.5 performance frequently exceeds that by using the default FD of bvp4c. When running under MATLAB 7.0 performance with AD or FD is degraded compared to that obtained with MATLAB 6.5. Under MATLAB 7.0 FD was more efficient for these problems. We attribute the relative success of FD over AD under MATLAB 7.0 as being due to the success of the JIT compiler in improving FD performance. JIT is unable to improve performance of the overloaded arithmetic used by MAD.

We are currently working on a source transformation implementation of AD in MATLAB called MATLAB Source transformation AD (MSAD) [Kha04,Kha05]. Source transformation removes the overhead of overloading and should permit the JIT compiler to optimise the derivative code.

Bibliography

FE04: Shaun A Forth and Marcus M. Edvall.
User Guide for MAD - MATLAB Automatic Differentiation Toolbox TOMLAB/MAD, Version 1.1 The Forward Mode.
TOMLAB Optimisation Inc., 855 Beech St 12, San Diego, CA 92101, USA, Jan 2004.
See http://tomlab.biz/products/mad.
For05: Shaun A. Forth.
An efficient overloaded implementation of forward mode automatic differentiation in MATLAB.
Accepted ACM TOMS, 2005.
Gri00: Andreas Griewank.
Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation.
Number 19 in Frontiers in Appl. Math. SIAM, Philadelphia, Penn., 2000.
Kha04: Rahul V. Kharche.
Source transformation for automatic differentiation in MATLAB.
Master's thesis, Cranfield University (Shrivenham Campus), Applied Mathematics & Operational Research Group, Engineering Systems Department, RMCS Shrivenham, Swindon SN6 8LA, UK, Aug. 2004.
Kha05: Rahul V. Kharche.
Progress with MATLAB Source Transformation AD - msad.
Presented at the 1st European Workshop on Automatic Differentiation, Maison du Seminaire, Nice, France, 14th-15th April 2005.
http://www-sop.inria.fr/tropics/adNice2005/slides/kharcheWkNi05.pdf.
KS01: Jacek Kierzenka and Lawrence F. Shampine.
A BVP solver based on residual control and the MATLAB PSE.
ACM Trans. Math. Softw., 27(3):299-316, 2001.
http://doi.acm.org/10.1145/502800.502801.
SGT03: L.F. Shampine, I. Gladwell, and S. Thompson.
Solving ODEs with MATLAB.
Cambridge University Press, 2003.
SKF05: L. F. Shampine, Robert Ketzscher, and Shaun A. Forth.
Using AD to solve BVPs in MATLAB.
ACM Trans. Math. Softw., 31(1):79-94, March 2005.
http://doi.acm.org/10.1145/1055531.1055535.
SKR: L. Shampine, J. Kierzenka, and M. Reichelt.
Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c.
http://www.mathworks.com/bvp_tutorial.

About this document ...

bvp4cAD- An Automatic Differentiation Enabled Boundary Value Solver

This document was generated using the LaTeX2HTML translator Version 2002 (1.62)

The command line arguments were:
latex2html bvp4cAD -dir /home/forth/www/AD/ADODE/bvp4cAD -split 0 -local_icons

The translation was initiated by Sean Forth AMORG on 2005-08-30

Sean Forth AMORG 2005-08-30

bvp4cAD- An Automatic Differentiation Enabled Boundary Value Solver