Q&A with Dr. Lawrence Shampine - Part II

Here is part II (see part I here) of the Q&A with Dr. Lawrence Shampine.

Q2:

In the absence of an analytical Jacobian, you have mentioned that accurate approximations of the Jacobians are difficult to compute using finite differences. Can you please elaborate? More specifically, is this a function of higher order equations or is it related to stiffness of an ODE?

A:  No, it is just hard to approximate derivatives accurately and reliably by finite differences.  The increment in the independent variable must be small enough to provide a reasonably accurate approximation to the derivative.  On the other hand, if it is too small, there is a loss of significance in forming the difference.  Indeed, if increment is too small, the function values are not different at all and the derivative is approximated as zero.  This is bad enough, but a Jacobian matrix is formed numerically a column at a time using a vector of function values formed in a call to the function defining the differential equations.  This means that the same increment is used in the approximation of all components of the column.  This increment might be too big to approximate accurately some partial derivatives and too small for others.  The stronger algorithms test the changes in the function values and adjust the increment if necessary to obtain a reasonable approximation to the larger derivatives.  Because the small partial derivatives in the column might be approximated poorly, finite differences should be used only when this is tolerable.
 

Q3:

In the paper “Implementation of Rosenbrock Methods”, you suggest that it is best to recondition the E matrix, and that it is easier to work with f_y – (2/h)*I as opposed to using E = I – hγf_y. In the implementation of a Rosenbrock type method, this is achieved by using (I/hγ-f_y) for the matrix inversion. What is the reason for this scaling?

A:  It is slightly less work to use the recommended form, but the real issue is one of scale.  When the problem is very stiff, there are eigenvalues of the Jacobian for which the product of the step size and eigenvalue is very large compared to one.  Accordingly, the term involving the Jacobian in the form E may be so big compared to I that when the subtraction is done in floating point arithmetic, the role of the identity matrix is lost.

Q&A with Dr. Lawrence Shampine - Part I

Prof. Lawrence Shampine of Southern Methodist University is a world renowned expert in Numerical Integration Methods, and has researched extensively on the Rosenbrock techniques. I recently had the privilege of corresponding with Dr. Shampine to get a few topics clarified and answered on the Rosenbrock techniques, specifically related to the definition of the Jacobians. Dr. Shampine graciously agreed to me posting these Q&A’s on my blog. Here is part I of the correspondence.

Q1:

How does the definition of the Jacobian differ for Rosenbrock and W-methods?

A:

Rosenbrock methods are one-step methods that use a fixed number of values of the Jacobian matrix at each step.  Many can be viewed as coming from the approximate evaluation of an implicit formula by a fixed number of iterations of Newton's method.  When evaluating an implicit formula, the Jacobian is used only to enhance convergence of the iteration, so it need not be evaluated exactly. This is important in practice because the default in many programs is to approximate it by finite differences and often some of the components are not at all accurate.  Moreover, many popular programs based on implicit formulas save approximate Jacobians and reuse them as long as the rate of convergence is acceptable.  A Rosenbrock method is much simpler to implement, but its order depends on having exact values of the Jacobian, which might be computed from an analytical expression or numerically with automatic differentiation or complex-step approximation.  Wolfbrandt devised methods, now called W-methods, that look just like Rosenbrock methods, but the order does not require the matrix to be exactly equal to the Jacobian.  Because it has more information at its disposal, a Rosenbrock method of the same number of evaluations will be more accurate, but if you like the convenience of finite difference Jacobians, you do have to worry about whether they are sufficiently accurate that a Rosenbrock method will behave as expected.  A W-method is not sensitive to the quality of the approximate Jacobian.

Reformulating the W-Method to avoid matrix multiplication

The W-method as presented by Steihaug & Wolfbrandt involves the matrix multiplication of the Jacobian with the ki’s at each step. Shampine (Implementation of Rosenbrock Algorithms, ACM Transactions on Mathematical Software, Vol. 8, No. 2, 1982, pp 93-113) also suggests “conditioning” the A matrix to avoid difficulties during the LU-decomposition. Shampine recommends working with the following format for A:

 

To avoid the matrix multiplication on the right hand side, the following re-arranged format is typically used for coding the 2-stage Rosenbrock method:

 

An intermediate step is required where the matrix g2 is calculated, from which k2 is later determined.