"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:
- 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.
- 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.
- 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.