Re: Working on bvp4c

[Bill Greene]

Here are a few thoughts on numjac.

numjac implements an algorithm for finite difference (FD) step size selection based on

this paper by Salane:

D. E. Salane, Adaptive routines for forming Jacobians numerically, Tech. Report SAND86

1319, Sandia National Laboratories, Albuquerque, NM, 1986.

I have done some experimentation with this algorithm and numjac and am much less

positive on this approach than either Salane or the numjac authors; I have seen

cases where the algorithm prescribes very tiny step sizes that lead to large

roundoff errors.

A second part of numjac is that, for sparse matrices (which I assume applies

to the bvp4c algorithm) it is usually necessary to make far fewer than n 

(number of equations) calls to the residual function when evaluating the

jacobian by FD. This substantially improves performance.

A discussion of this can be found in this paper:

T. F. Coleman, B. S. Garbow, and J. J. More, Software for estimating sparse Jacobian

matrices, ACM Transactions on Mathematical Software, 11 (1984), pp. 329-345.

Unfortunately, implementing this is not trivial.

A possible alternative to both of these issues is the approach used by the algorithm

in Sundials (http://computation.llnl.gov/projects/sundials-suite-nonlinear-differential-algebraic-equation-solvers)

for calculation of jacobians by FD when the matrix is banded (I'm

assuming this applies to bvp4c). The routine is idaDlsBandDQJac and it is discussed

somewhat in the Sundials documentation. A downside of this (other than the

implementation being in C) is that, as far as I am aware, Octave (and MATLAB) has

no specific support for banded matrices.

On Sat, Aug 6, 2016 at 9:38 AM, lakerluke <address@hidden> wrote:

Thanks for the pointer! That really helps. So don't worry about answering my
questions from the previous post. In light of your comment I think I
understand the approach intended by the Shampine paper. We use Newton's
method to look for roots for the system (4.1) which gives us our solution at
mesh points. As pointed out in the paper in order to apply this method we
must compute the jacobian of the terms in the system (4.1). This in turn
requires us to compute the partial derivatives of the f and g functions. I
think the paper then goes on to explain that in order to form these partial
derivative terms we consider the Jacobians of f (I will refer to these
Jacobians as [1]):

J_{i} = (partial f_{i}) / (partial y)

J_{i - 1/2} = (partial f_{i - 1/2}) / (partial y)

We then compute these using finite difference to then use in the global
Jacobian (printed in the middle of page 7):

(partial phi_{i}) / (partial y_{i})

Now in order to form the Jacobian terms [1] the paper says they use the
function numjac to calculate the partial derivatives. I believe Octave does
not have such a function...

Does anyone know of any equivalent method in octave to numjac? Do you think
this is something I am going to have to implement first?

Luke



--
View this message in context: http://octave.1599824.n4.nabble.com/Working-on-bvp4c-tp4677540p4679048.html

Sent from the Octave - Maintainers mailing list archive at Nabble.com.