## Random, Linear, Homogeneous Difference Equations and Time Integration Methods for ODEs

Most methods (that I know of) of numerically approximating the solution of ODEs are "general linear methods". For this type of method, the so-called 'linear stability' is examined by applying the method to the linear, constant-coefficient, (complex) scalar ODE

$\dot{y} = \lambda y$.

Because the same procedure is used to generate a new approximation at each time step, the methods result in linear, homogeneous, constant-coefficient difference equations for the approximate solution values at each time step. As a result, linear stability analysis for these methods amounts to analyzing the solutions of constant coefficient, homogeneous linear difference equations.

If you vary the method randomly at each time step, the coefficients of the difference equation are no longer constant. For example, second order explicit linear multistep methods can be written as a one-parameter family. If you choose this parameter randomly at each time step, the difference equation looks like this

$y_{n+1} + F(a(n))y_n + G(a(n))y_{n-1} = 0$

where $a(n)$ is a random variable (the parameter in the family of methods) and the functions $F$ and $G$ are known. My main question is whether there is any theory giving conditions on $F$, $G$, and the distribution of $a$ such that solutions remain bounded in the limit of large $n$. It would be nice if the theory generalized to higher order difference equations too.

If you were to randomly decide to use second-order multistep versus second-order Runge-Kutta methods, for example, then $F$ and $G$ would also depend on $n$. A theory to handle that case would be welcome too.

Since the numerical analyst is free to choose the distribution of $a$, and to some extent also $F$ and $G$, I'm wondering if it might be possible to design 'random' methods that have better linear stability properties than the usual ones which repeat the same process over and over. I'm posting here because I know almost nothing about stochastic/random processes.

-

The usual way of analyzing three-term recurrences is to look at the asymptotic behavior of your $F(a(n))$ and $G(a(n))$ and then invoke the machinery of (PoincarĂ©-)Perron. What often happens is that there is one solution that is decaying (the "minimal" solution) and one that is growing (the "dominant" solution).

Usually people are interested in the bounded or minimal solution, and attempting to propagate the recurrence forward (increasing n) in finite precision is courting certain numerical disaster. Due to rounding, the initial values of your recurrence can be expressed as a linear combination of the minimal and dominant solution; thus, if the recurrence is executed forward, there will come a point where the growth of the dominant solution will swamp the behavior of the bounded solution. (In numerical work, people usually recurse backwards or use special techniques like the Miller algorithm when they are interested in minimal solutions)

Related to this, it has long been known that certain multistep methods for solving ODEs (e.g. Milne's method) can exhibit the phenomenon of "parasitic solutions", where no matter how small a step size $h$ you take, a growing solution to the recurrence that alternates in sign with every step can easily contaminate your bounded solution, rendering your numbers meaningless.

At this point, I would like to direct you to this excellent survey article by Walter Gautschi and this old but nice book by Jet Wimp.

-