4. Fast Implementation of Picard Iteration for Systems of Linear ODEs with Periodic Coefficients and the FastPicard Package

By constructing an alternative function, we will see that the main reason for the slowness of picard is that it uses the built-in function Integrate for integration. Integrate is of course very good for general purposes, but it is too slow in cases such as ours, in which we have to calculate a very large number of integrals of functions in a limited class. As shown in the previous section, we do not need all the generality of Integrate because all integrands belong to the family .

We created a new function, specialized for definite and indefinite integration of functions in , to be used instead of Integrate in the Picard iteration scheme. We called it NewIntegrate and gave it the same syntax as Integrate. It is implemented in our FastPicard package (see Additional Material) and exported by it.

First, NewIntegrate recognizes a linear combination of terms and uses the fact that the integral is a linear operator. This is a standard trick used, for example, in the implementation of Laplace transforms [9, 13]. It then readily performs indefinite integrals of constants, powers, sines, and cosines. The more difficult functions in are dealt with by using the recursive formulas of equations (8) and (9). As usual in these cases, we track all results, so that the same integral need not be calculated more than once. This is accomplished by the old trick of using a delayed and an immediate assignment, as explained in Section 2.4.9 of [14].

In calculating definite integrals, a trick for further speed is to use instead of as the expression for the indefinite integral of . In most cases, we will be interested in using in equation (2). This trick ensures that all definite integrals of functions in with a lower limit equal to zero will coincide with the corresponding indefinite integrals.

The reader is invited to test the accuracy and performance of NewIntegrate by producing a large random function in and integrating it with both functions. In our tests, we found that NewIntegrate is tens of times faster than the built-in function.

The package also exports the function FastPicard[matrix, tzero, xzero, n] that is similar to picard but uses NewIntegrate. For functions in to be recognized by NewIntegrate, it is necessary to transform products of sines and cosines into sums before integration. As we are always dealing with linear systems, the syntax is a bit simplified--instead of providing the function , the reader should provide matrix .

Of course, times are much shorter than the ones obtained by picard. This is shown in the following example, which also illustrates the syntax.

Besides using NewIntegrate, FastPicard incorporates other improvements we discovered after some experimenting. Here is a list of the differences.

The package also exports the functions FundamentalMatrix and FloquetMatrix.

A fundamental matrix solution to a linear system such as equation (4) is a matrix whose columns are solutions to the system with linear independent initial conditions. The function FundamentalMatrix[matrix, tzero, n] calculates Picard iterates approximating the particular fundamental matrix solution that takes the columns of the identity matrix as initial conditions at time tzero. The second argument tzero is optional with a default value of 0. For example, with 20 iterates and initial time 0, an approximate fundamental matrix solution to the system with

is obtained by

One advantage of having a fundamental matrix solution calculated at initial time is that with it we may recover solutions to the system with any desired initial conditions at time . This is even simpler when we use the particular fundamental matrix solution calculated by FundamentalMatrix. For example, this gives the approximate numerical value for the solution of this system with at time and initial condition .

These values can be compared with the corresponding exact ones. In fact, the exact fundamental matrix solution for the system is given in [1]. Using it, we obtain

The Floquet transition matrix for a linear system of differential equations with periodic coefficients of period is just the special fundamental matrix solution calculated by FundamentalMatrix with initial time 0 evaluated at time . In the next section we will see an interesting application for this matrix. An approximation to it using Picard iterates is calculated by FloquetMatrix[matrix, per, n], where the second argument, the period of the matrix, is optional with a default value of .

For example, here is an approximation to the Floquet transition matrix for the Mathieu equation with and .