"A Natural Approach to the Numerical Integration of Riccati Equations" - Programs Used

The aim of this page is to give a little more technical information about the programs used to generate the data in the Riccati equation paper written by Steve Shnider and myself.

A Riccati equation is specified by its coefficient matrix A(t). From this we need to compute Gamma(t,h), which is the matrix used to compute the numerical approximation to y(t+h) given y(t). (If you're looking at our paper you should have no trouble working out our notation here!) There are three approaches to programming this up:

  1. The program should contain a routine to compute A(t) for given t. It should then be able to find Gamma(t,h) when given h, and thus y(t+h) when given y(t).
    Advantage: Simple to apply; just write a routine to find A(t) and run!
    Disadvantage: Lots of matrix additions and multiplications are necessary to compute Gamma(t,h) from A(t), which weighs the program down a bit.
  2. A formula for Gamma(t,h) should be found using Maple, and translated into c using the Maple -> C translation facility. The program doesn't need a routine for A(t).
    Advantage: Run time goes down.
    Disadvantage: Preparation time goes up, a lot. The Maple -> C translator is currently limited, and so you need to do quite some playing around with the output to get your final program.
  3. If you're using Maple already, then you might as well get it to give you a general expression for y(t+h) starting from y(t),t and h.
    Advantage: Run time really down. And since y(t+h) is a smaller matrix than Gamma(t,h), in some ways it is easier to use the Maple -> C translator for this.
    Disadvantage: Trouble is, whereas Gamma(t,h) was a function of two scalars, t and h, y(t+h) depends on t,h, and y(t), the latter being a matrix. The Maple -> C translator doesn't seem to me to be able to write good code using pointers for such things. Maybe I could make it do this if I knew how.
Overall, as you increase the method number, you reduce run time but increase preparation time. Much of the preparation could in priciple be automated, so if you happen to want to solve Riccati equations for the next few years, you should probably aim for the last method.

In the examples in our paper, we used approach 1 in example 1, example 3 (for 2nd order methods) and example 4. For example 2, a scalar equation, we used approach 3 - the disadvantage here isn't any more. For example 3, using 4th order methods, we used approach 2. 4th order methods require way more multiplications and additions to get Gamma(t,h) from A(t), so we definitely saved by doing these "once and for all" via sybolic manipulator before running any programs.