3. Picard Iteration for Systems of Linear ODEs with Periodic Coefficients

In this section we show that for systems of linear differential equations with periodic coefficients we never encounter integrals that are impossible to calculate. But we will also see that evaluating Picard iterates, even if technically possible, may take too long if we use the built-in Integrate command.

More exactly, we will be concerned with systems of the form

where is a matrix satisfying

for all . As a simplification, we restrict ourselves a bit more by considering the case where is a finite sum of terms proportional to sines and cosines of period , that is,

In Section 5 of this article we will see some very interesting examples belonging to this class, but any reasonable periodic matrix with period can be cast approximately in the previous form by retaining a finite number of terms in its Fourier series.

An important example is the Mathieu equation

which can be rewritten in the form (4) by making and and taking

Although the Mathieu equation is exactly solvable with Mathematica, we will use it as our prototype example for its simplicity and importance.

We may also compare our approximate solutions with the exact ones. It should be noted that Mathieu functions, that is, the solutions of the Mathieu equation, are notably difficult to implement with a computer. For this reason they are the subject of recent research [10, 11, 12].

Let us evaluate some Picard iterates approximating the solution of the Mathieu equation, so that we can learn something from them. Taking an initial condition such as , here are the first, second, and third iterates.

The reader can see that the expanded form of all these iterates is a sum of terms that either are constants, positive integer powers of , sines or cosines, or assume one of the forms or , where is a positive integer. Let us define as the family of the finite linear combinations of all such functions and see whether the next iterate still belongs to the family .

If we calculate the next iterate by hand, we first have to multiply matrix by the last element in the previous list and then integrate.

Keeping in mind the possible appearance of difficult integrals, the main difficulty at this point is the appearance of products involving sines and cosines when we expand the last expression. But we may transform such products into sums of trigonometric functions by using TrigReduce, which again brings our expression back to the explicit form of a function in . We still need to integrate this expression. That may be done by using the recursive formulas

obtained by an easy integration by parts. We thus see that all integrals of functions in always belong to .

Although we are not interested in giving a formal proof here, all the reasoning with the third iterate can be turned into a proof by induction of the following result.

Theorem: Any Picard iterate of linear systems in the form of equation (4) with in the form of equation (6) is a function in the family .

Although Picard iteration is thus technically possible at any finite order, as the order grows the resulting expressions become increasingly complicated. As a consequence, the time needed to compute them increases quickly.

Although one might assume that eight Picard iterates are enough, we caution the reader that for the applications we have in mind it will be necessary to evaluate the solution to the Mathieu equation at time , that is, at the end of the period of the matrix of the system. To see that eight iterations are too few, we recommend comparing the graphs of the approximate solution obtained by Picard iteration (substituting and by typical values 0.9 and 0.6, respectively) and the exact solution given by DSolve.