A Filtered Milne-Simpson ODE Method with Enhanced Stability Properties

Ari Aluthge, Scott A. Sarra

Marshall University, Huntington, WV, USA.

DOI: 10.4236/jamp.2023.111013                

Abstract

The Milne-Simpson method is a two-step implicit linear multistep method for the numerical solution of ODEs that obtains the theoretically highest order of convergence for such a method. The stability region of the method is only an interval on the imaginary axis and the method is classified as weakly stable which causes non-physical oscillations to appear in numerical solutions. For this reason, the method is seldom used in applications. This work examines filtering techniques that improve the stability properties of the Milne-Simpson method while retaining its fourth-order convergence rate. The resulting filtered Milne-Simpson method is attractive as a method of lines integrator of linear time-dependent partial differential equations.

Keywords

Ordinary Differential Equations, Numerical Methods, Weak Instability

1. Introduction

The first Dahlquist barrier [1] states that a zero-stable and linear s-step multistep method cannot attain an order of convergence greater than $s+1$ if s is odd and greater than $s+2$ if s is even. If the method is also explicit, then it cannot attain an order greater than s. The Milne-Simpson (MS) method is an implicit 2-step Linear Multistep Method (LMM) that realizes the optimal fourth-order convergence rate limit given in the first Dahlquist barrier. However, the MS method is only weakly stable which severely limits its use in applications. In this work, filters are examined that improve the stability properties of the MS method. The filters were first introduced in [2] and [3] where they were arbitrarily applied every ten-time steps. However, in these works, there was no examination of how the application of the filters affected the stability region of the method. In this work, we show that the frequency of the application of the filters can be used to modify the location of the stability region of the method which allows the technique to be tailored to the type of ODE system that is being solved.

In previous work [4], similar filters were used to improve the stability properties of the weakly stable explicit midpoint rule (also called the leapfrog method) while maintaining the second-order accuracy of the method. The filters used with the midpoint rule were of lengths three and five. Using filters to improve the stability properties of the MS method while maintaining fourth-order convergence is more difficult with the MS method than with the midpoint rule due to the MS method being implicit and also due to the necessity that the filter must be of length seven.

In the next section, LMMs and their stability properties are summarized. Then, a seven-point filter is developed that allows the MS method to retain its fourth-order convergence rate while improving the stability properties of the method and insight is given as to how to most effectively apply the filter. Finally, numerical examples are carried out that demonstrate the effectiveness of the filter.

2. Linear Multistep Methods

An s-step linear multistep method for the numerical solution of the Ordinary Differential Equation (ODE) Initial Value Problem (IVP)

$y^{\prime }=f\left(t,y\left(t\right)\right),y\left(0\right)=y^0$ (1)

is of the form

${\displaystyle \underset{m=0}{\overset{s}{\sum }}}{\alpha }_m{y}^{n+m}=k{\displaystyle \underset{m=0}{\overset{s}{\sum }}}{\beta }_{m}f\left(t^{n+m},y^{n+m}\right),n=0,1,\cdots$ (2)

where ${\alpha }_m$ and ${\beta }_m$ are given constants and k is the size of the time step. Multistep methods use information from the previous s steps to calculate the next value. It is conventional to normalize (2) by setting ${\alpha }_s=1$. When ${\beta }_s=0$ the method is explicit. Otherwise it is implicit. In order to start the method, the first $s-1$ time levels need to be calculated by a one-step method.

The stability properties of LMMs are examined by applying the method to the ODE

$y^{\prime }=\lambda y,y\left(0\right)=y^0$ (3)

where $\lambda$ is a complex number with $\Re \left(\lambda \right)\le 0$. Equation (3) is called the Dahlquist test equation. The first concepts of stability result from considering the simple case of $\lambda =0$ in (3) for which the solution is the constant $y^0$. In this case, the application of the LMM (2) to the test equation results in a difference equation with a general solution of the form

$y^n=c_1{\zeta }_1^n+\cdots +c_s{\zeta }_s^n$ (4)

where ${\zeta }_i$ are the zeros of the first characteristic polynomial

$\rho \left(\zeta \right)={\displaystyle \underset{m=0}{\overset{s}{\sum }}}{\alpha }_m{\zeta }^m$ (5)

of the LMM (2). For a consistent method, one of the terms, say ${\zeta }_i^n$, approaches the solution of (3) as $n\to \infty$ while the other roots correspond to extraneous solutions. The following stability definitions quantify to what extent the computed solutions of (4) remain bounded. A method is stable (also called zero stable or strongly stable) if $\rho$ has a simple root equal to one and if all the other roots of $\rho$ are less than one in absolute value. The criteria that the roots must satisfy for the method to be stable are called a root condition. A weakly stable method has all roots of $\rho$ satisfying $\left|\zeta \right|\le 1$ and has at least one root of absolute value one other than the simple root $\zeta =1$. The solutions may remain bounded for a fixed interval of time $\left[\mathrm{0,}T\right]$ but will eventually become unbounded as $t\to \infty$. A method for which a root of $\rho$ is greater than one in magnitude is unstable. In this case the solutions will rapidly become unbounded.

A stronger definition of stability can be made by considering non-zero $\lambda \in \mathcal{C}$ with negative real part in Equation (3). In this case the exact solution of (3) is $y\left(t\right)=y^0{\text{e}}^{\lambda t}$ which approaches zero as $t\to \infty$ since the real part of $\lambda$ is negative. A method that produces numerical solutions with the same asymptotic behavior, that is $y^n\to 0$ as $t^n\to \infty$, for a fixed value of k is said to be absolutely stable. This type of stability has also been called asymptotic stability, or eigenvalue stability. Absolute stability may also be characterized in terms of a stability function or amplification factor $R\left(z\right)$ which satisfies $y^{n+1}=R\left(z\right)y^n$ where $z=k\lambda$. Then the absolute stability region is $\mathcal{S}=\left\{z\in \mathcal{C}\mathrm{<mo>\; :\; </mo>}\left|R\left(z\right)\right|\le 1\right\}$. For z in the stability region, the numerical solutions exhibit the same asymptotic behavior as does the exact solution of (3).

For a linear system of ODEs, a necessary condition for absolute stability is that all of the scaled (by the time step size k) eigenvalues $z=k\lambda$ of the coefficient matrix of the system must lie in the stability region. For nonlinear systems, the system must be linearized and the scaled eigenvalues of the Jacobian matrix of the system must be analyzed at each $t^n$. Standard references on numerical ODE methods [5] - [10] can be consulted for more details.

3. The Milne-Simpson Method

The Milne-Simpson method for the system (1) is

$y^{n+1}=y^{n-1}+\frac{k}{3}\left[f^{n+1}+4f^n+f^{n-1}\right].$ (6)

The method is constructed by integrating the ODE (1) on the interval $\left[t^n\mathrm{,}{t}^{n+1}\right]$ to get

$y\left(t^{n+1}\right)=y\left(t^n\right)+{\displaystyle {\int }_{t^n}^{t^{n+1}}}f\left(t,y\right)\text{d}t$

and then approximating the integral by the well-known Simpson’s rule that is usually introduced in the Calculus sequence. The resulting implicit method has fourth-order algebraic convergence. The Milne-Simpson method is of interest since it is the only 2-step linear multistep method that attains its optimal order [8]. Any optimal linear multistep method has no region of absolute stability [8]. The MS method needs two time levels to start, $y^0$ which is provided by the initial condition and $y^1$. Typically, the time level $y^1$ is obtained by an explicit fourth-order Runge-Kutta (RK) method.

The first characteristic polynomial of the Milne-Simpson method is $\rho \left(z\right)=z^2-1$ with roots $\zeta =\pm 1$ and the method is only weakly stable. The stability region only consists of the interval $\mathcal{S}=\left[-i\sqrt{3}\mathrm{,}i\sqrt{3}\right]$ on the imaginary axis.

After discretizing a system of nonlinear ODEs, the resulting system of nonlinear algebraic equations at each time step is solved using Newton’s method with the previous time level used as an initial guess. In the numerical examples, one or two iterations with Newton’s method was typically required at each time step. Each iteration requires one evaluation of $f\left(t\mathrm{,}y\right)$ from the right side of (1) and the solution of a linear algebraic system using LU factorization.

Linear systems of ODEs can be more efficiently evaluated using the MS method. For a constant coefficient linear system of ODEs, $y^{\prime }=Ay$, where A is a $N\times N$ coefficient matrix the MS method can be advanced in time by solving the algebraic linear system

$L{y}^{n+1}=b$ (7)

where

$L=I-\frac{k}{3}A$ (8)

I is the $N\times N$ identity matrix, and

$b=y^{n-1}+\frac{k}{3}\left[4A{y}^n+A{y}^{n-1}\right].$

The matrix L must be factorized once at the beginning of the computation with a $\mathcal{O}\left(N^3\right)$ LU factorization. Then, at each time step one matrix-vector multiplications and a forward and back substitution to solve (7) are required for a total flop count of $4N^2$. In comparison, an explicit fourth-order RK method requires four matrix-vector multiplications that require $8N^2$ flops. The process for linear variable coefficient systems of ODEs is similar except that the matrix L must be factorized at each time step.

4. Filter Construction

It is shown in [3] that the numerical solution of the IVP (1) obtained by the Milne-Simpson method using starting values only assumed to have an asymptotic expansion of the very general form

${\eta }_n\left(k\right)=y^n+{\displaystyle \underset{m=4}{\overset{5}{\sum }}}{\eta }_{mn}{k}^m+\mathcal{O}\left(k^6\right),n=0,1$ (9)

has an asymptotic expansion of the form

$y^n=\text{physical}\text{mode}\left(y_p\right)+\text{computational}\left(y_c\right)+\mathcal{O}\left(k^6\right).$ (10)

The physical mode $y_p$

$\begin{array}{c}{y}_p=y^0{\text{e}}^{\lambda t}+\left[\left(\frac{{\lambda }^5}{180}{y}^{0}t+\frac{1}{2}\left({\eta }_{40}+{\eta }_{41}\right)\right){\text{e}}^{\lambda t}\right]k^4\\ +\left[\frac{1}{2}\left(-\frac{{\lambda }^5}{180}{y}^0-\frac{2}{3}{\eta }_{40}-\frac{1}{3}{\eta }_{41}+{\eta }_{50}+{\eta }_{51}\right){\text{e}}^{\lambda t}\right]k^5\end{array}$ (11)

consists of the exact solution $y^0{\text{e}}^{\lambda t}$ and a $\mathcal{O}\left(k^4\right)$ truncation error that goes to zero as $k\to 0$. The computational mode

$\begin{array}{c}{y}_c=\left[\frac{{\left(-1\right)}^n}{2}\left({\eta }_{40}-{\eta }_{41}\right){\text{e}}^{-\lambda t/3}\right]k^4\\ +\left[\frac{{\left(-1\right)}^n}{2}\left(\frac{{\lambda }^5}{180}{y}^0+\frac{2}{3}\lambda {\eta }_{40}+\frac{1}{3}{\eta }_{41}+{\eta }_{50}-{\eta }_{51}\right){\text{e}}^{-\lambda t/3}\right]k^5\end{array}$ (12)

contributes a non-physical oscillatory factor to the solution which grows in size with increasing t. References [3] and [11] can be consulted for the details of the derivation of the expressions for $y_p$ and $y_c$.

In order to damp the computational mode while retaining the fourth-order accurate approximation of the physical mode a filter of the form

$P\left(E\right)={\displaystyle \underset{j=-r_1}{\overset{r_2}{\sum }}}{a}_j{E}^j$ (13)

is considered where $a_j\in ℝ$, $r_1\mathrm{,}{r}_2\in {\mathbb{Z}}_{+}$, and E is the forward shift operator defined by $Ey\left(t\right)=y\left(t+k\right)$ in continuous notation and $E{y}^n=y^{n+1}$ in discrete notation. Let $s=\lambda k$ and then elementary calculus gives

$P\left(E\right){\text{e}}^{\lambda t}={\text{e}}^{\lambda t}P\left({\text{e}}^s\right)\mathrm{,}$ (14)

$P\left(E\right)t{\text{e}}^{\lambda t}={\text{e}}^{\lambda t}\left[tP\left({\text{e}}^s\right)+k\frac{\text{d}}{\text{d}s}P\left({\text{e}}^s\right)\right],$ (15)

and for any $\mu \in \mathcal{R}$

$P\left(E\right){\left(-1\right)}^n{\text{e}}^{-\mu \lambda t}={\left(-1\right)}^n{\text{e}}^{-\mu \lambda t}P\left(-{\text{e}}^{-\mu s}\right)\mathrm{.}$ (16)

Applying the filter (13) to the physical mode (11) results in

$\begin{array}{c}{\hat{y}}_p=P\left(E\right)y_p\\ =y^0{\text{e}}^{\lambda t}P\left({\text{e}}^s\right)+\left[\left(\frac{{\lambda }^5}{180}{y}^0\left(tP\left({\text{e}}^s\right)+k\frac{\text{d}}{\text{d}s}P\left({\text{e}}^s\right)\right)+\frac{1}{2}\left({\eta }_{40}+{\eta }_{41}\right)P\left({\text{e}}^s\right)\right){\text{e}}^{\lambda t}\right]k^4\\ +\left[\frac{1}{2}\left(-\frac{{\lambda }^5}{180}{y}^0-\frac{2}{3}{\eta }_{40}-\frac{1}{3}{\eta }_{41}+{\eta }_{50}+{\eta }_{51}\right){\text{e}}^{\lambda t}P\left({\text{e}}^s\right)\right]k^5\end{array}$ (17)

and applying the filter to the computational mode (12) gives

$\begin{array}{c}{\hat{y}}_c=P\left(E\right)y_c\\ =\left[\frac{{\left(-1\right)}^n}{2}\left({\eta }_{40}-{\eta }_{41}\right){\text{e}}^{-\lambda t/3}P\left(-{\text{e}}^{-s/3}\right)\right]k^4\\ +\left[\frac{{\left(-1\right)}^n}{2}\left(\frac{{\lambda }^5}{180}{y}^0+\frac{2}{3}\lambda {\eta }_{40}+\frac{1}{3}{\eta }_{41}+{\eta }_{50}-{\eta }_{51}\right){\text{e}}^{-\lambda t/3}P\left(-{\text{e}}^{-s/3}\right)\right]k^5.\end{array}$ (18)

Equation (17) suggests that in order for the physical mode to be kept unchanged up to the order $k^N$ terms that P must satisfy

$P\left({\text{e}}^s\right)=1+\mathcal{O}\left(s^{N+1}\right)\mathrm{.}$ (19)

In a similar manner, Equation (18) reveals that in order for the computational mode to be reduced by a factor of $k^M$ that P must satisfy

$P\left(-{\text{e}}^{-s}\right)=\mathcal{O}\left(s^M\right)\mathrm{.}$ (20)

Equation (19) implies that

$\frac{\text{d}}{\text{d}s}P\left({\text{e}}^s\right)=\mathcal{O}\left(s^N\right).$ (21)

Substituting Equations (19), (20), and (21) into Equations (17) and (18), results in

$\begin{array}{c}{\hat{y}}_p=y^0{\text{e}}^{\lambda t}\left(1+\mathcal{O}\left(s^{N+1}\right)\right)\\ +\left[\left(\frac{{\lambda }^5}{180}{y}^0\left(t\left(1+\mathcal{O}\left(s^{N+1}\right)\right)+k\mathcal{O}\left(s^N\right)\right)+\frac{1}{2}\left({\eta }_{40}+{\eta }_{41}\right)\left(1+\mathcal{O}\left(s^{N+1}\right)\right)\right){\text{e}}^{\lambda t}\right]k^4\\ +\left[\frac{1}{2}\left(-\frac{{\lambda }^5}{180}{y}^0-\frac{2}{3}{\eta }_{40}-\frac{1}{3}{\eta }_{41}+{\eta }_{50}+{\eta }_{51}\right){\text{e}}^{\lambda t}\left(1+\mathcal{O}\left(s^{N+1}\right)\right)\right]k^5\end{array}$ (22)

and

$\begin{array}{c}{\hat{y}}_c=\left[\frac{{\left(-1\right)}^n}{2}\left({\eta }_{40}-{\eta }_{41}\right){\text{e}}^{-\lambda t/3}\mathcal{O}\left(s^M\right)\right]k^4\\ +\left[\frac{{\left(-1\right)}^n}{2}\left(\frac{{\lambda }^5}{180}{y}^0+\frac{2}{3}\lambda {\eta }_{40}+\frac{1}{3}{\eta }_{41}+{\eta }_{50}-{\eta }_{51}\right){\text{e}}^{-\lambda t/3}\mathcal{O}\left(s^M\right)\right]k^5\end{array}$ (23)

respectively.

The term $\frac{{\lambda }^5}{180}{y}^{0}t\left(1+\mathcal{O}\left(s^{N+1}\right)\right)$ in (22) shows that the filter P is unable to improve the order of accuracy of the method. It suggests taking $N=4$. Additionally, so that extrapolation for improving the accuracy of the physical mode is applicable after applying the filter P, it is necessary to reduce the computational mode up to $\mathcal{O}\left(s^6\right)$ [8]. This suggest taking $M=2$. Thus,

$P\left({\text{e}}^s\right)=1+\mathcal{O}\left(s^5\right)$ (24)

and

$P\left(-{\text{e}}^{-s/3}\right)=\mathcal{O}\left(s^2\right)\mathrm{.}$ (25)

Equations (24) and (25) consist of seven linear equations which can be solved for seven unknowns (i.e. for any seven $a_j$ ’s). This set of unknowns may be written in the form

$\left\{a_{l-3}^{\left(l\right)},a_{l-2}^{\left(l\right)},a_{l-1}^{\left(l\right)},a_l^{\left(l\right)},a_{l+1}^{\left(l\right)},a_{l+2}^{\left(l\right)},a_{l+3}^{\left(l\right)}\right\},l=-3,-2,-1,0,1,2,3$

where, for example, the case $l=3$ is the case $r_1=6$, $r_2=0$ and the case $l=2$ is the case $r_1=5$, $r_2=1$ and so on. Elementary linear algebra shows that, in all of the above cases, the systems of linear equations are linearly independent and can be solved uniquely for seven unknowns. In each case, the result yields a unique filter involving seven terms.

A general filter may be written in the form

$P_l\left(E\right)={\displaystyle \underset{j=-3}{\overset{3}{\sum }}}{a}_{l+j}^{\left(l\right)}{E}^j,l=-3,-2,-1,0,1,2,3.$ (26)

Equations (24), (25), and (26) lead to the equations

${\displaystyle \underset{j=-3}{\overset{3}{\sum }}}{a}_{l+j}^{\left(l\right)}{\text{e}}^{sj}=1+\mathcal{O}\left(s^5\right),l=-3,-2,-1,0,1,2,3$ (27)

and

${\displaystyle \underset{j=-3}{\overset{3}{\sum }}}{a}_{l+j}^{\left(l\right)}{\left(-{\text{e}}^{-s/3}\right)}^j=\mathcal{O}\left(s^2\right),l=-3,-2,-1,0,1,2,3.$ (28)

For any given filter (any fixed value of $l$ ), Equation (27) contains five linear equations and Equation (28) contains two linear equations in s. Each equation contains seven unknowns. For example, for the filter $P_{-3}$ ( $l=-3$ ), the five equations resulting from (27) are:

$2a_{-6}^{\left(-3\right)}+a_{-5}^{\left(-3\right)}+a_{-4}^{\left(-3\right)}+a_{-3}^{\left(-3\right)}+a_{-2}^{\left(-3\right)}+a_{-1}^{\left(-3\right)}+a_0^{\left(-3\right)}=1$

$6a_{-6}^{\left(-3\right)}+5a_{-5}^{\left(-3\right)}+4a_{-4}^{\left(-3\right)}+3a_{-3}^{\left(-3\right)}+2a_{-2}^{\left(-3\right)}+a_{-1}^{\left(-3\right)}=0$

$36a_{-6}^{\left(-3\right)}+25a_{-5}^{\left(-3\right)}+16a_{-4}^{\left(-3\right)}+9a_{-3}^{\left(-3\right)}+4a_{-2}^{\left(-3\right)}+a_{-1}^{\left(-3\right)}=0$

$216a_{-6}^{\left(-3\right)}+125a_{-5}^{\left(-3\right)}+64a_{-4}^{\left(-3\right)}+27a_{-3}^{\left(-3\right)}+8a_{-2}^{\left(-3\right)}+a_{-1}^{\left(-3\right)}=0$

$1296a_{-6}^{\left(-3\right)}+625a_{-5}^{\left(-3\right)}+256a_{-4}^{\left(-3\right)}+81a_{-3}^{\left(-3\right)}+16a_{-2}^{\left(-3\right)}+a_{-1}^{\left(-3\right)}=0$

and the two equations resulting from (28) are:

$2a_{-6}^{\left(-3\right)}-a_{-5}^{\left(-3\right)}+a_{-4}^{\left(-3\right)}-a_{-3}^{\left(-3\right)}+a_{-2}^{\left(-3\right)}-a_{-1}^{\left(-3\right)}+a_0^{\left(-3\right)}=0$

$6a_{-6}^{\left(-3\right)}-5a_{-5}^{\left(-3\right)}+4a_{-4}^{\left(-3\right)}-3a_{-3}^{\left(-3\right)}+2a_{-2}^{\left(-3\right)}-a_{-1}^{\left(-3\right)}=0.$

The solution of the $7\times 7$ linear systems is $a_{-6}^{\left(-3\right)}=\frac{5}{64}$, $a_{-5}^{\left(-3\right)}=\frac{-18}{64}$, $a_{-4}^{\left(-3\right)}=\frac{15}{64}$, $a_{-3}^{\left(-3\right)}=\frac{20}{64}$, $a_{-2}^{\left(-3\right)}=\frac{-45}{64}$, $a_{-1}^{\left(-3\right)}=\frac{30}{64}$, and $a_0^{\left(-3\right)}=\frac{57}{64}$. This filter can be written as

$P_{-3}\left(E\right)=\frac{1}{64}\left(5E^{-6}-18E^{-5}+15E^{-4}+20E^{-3}-45E^{-2}+30E^{-1}+57E^0\right)$

or

${\hat{y}}^n=P_{-3}\left(E\right)y^n=\frac{1}{64}\left(5y^{n-6}-18y^{n-5}+15E^{n-4}+20y^{n-3}-45y^{n-2}+30y^{n-1}+57y^n\right)$

where $y^n$ is the unfiltered solutions and ${\hat{y}}^n$ is the filtered solution at $t^n=nk$. The coefficients of all seven filters are listed in Table 1.

Table 1. Seven-point filter coefficients.

Filter $P_7^{\left(0\right)}$ is called a centered filter as it requires three time levels greater than n and three less than n to filter at time level n. Filters $P_7^{\left(1\right)}$, $P_7^{\left(2\right)}$, and $P_7^{\left(3\right)}$ are forward biased filters as they require more time levels greater than n than less than n to filter at time level n. Filters $P_7^{\left(-1\right)}$, $P_7^{\left(-2\right)}$, and $P_7^{\left(-3\right)}$ are backward biased filters as they require fewer time levels greater than n and more less than n to filter at time level n. Table 2 lists the number of additional time levels that must be calculated and disregarded each time a filter is applied.

5. Filter Application

This section gives insight on how to use the filters to manipulate both the size and location of the absolute stability region of the MS method which without filtering is only an interval on the imaginary axis.

In applications, it is not necessary to filter every step. Rather, a filter is applied every $N_0$ time steps. The choice of $N_0$ is problem dependent as it determines both the size and the location of the absolute stability region of the method. There is no exact $N_0$, but rather a range of $N_0$ that is appropriate for a specific problem.

The method starts by setting $y^0$ equal to the initial condition and calculating $y^1$ via an explicit fourth-order Runge-Kutta method [5]. Then, $N_0-1$ Milne-Simpson steps are taken as well as the additional time levels required by the filter (Table 2) and then the level $N_0$ step is filtered.

The method is restarted and continues with the filtered $y^{N_0}$ as the initial condition and the following level calculated with RK4 and then additional levels are calculated with the MS method. This approach is needed in order to calculate the stability region of what is a composition of three methods: RK4, MS, and a filter. It can be done over a finite time interval, which is why such a restart is needed. Another approach is to restart with level $y^{N_0-1}$ as the first level and the filtered $y^{N_0}$ as the second level in the MS method which omits the RK4 step in the restart after filtering. With this approach, an exact stability region can not be determined for what is now a composition of methods on an infinite time interval. However, if this simpler approach is taken the stability regions from the method with RK4 in the startup still provide a good sense of the stability region. The approach without the RK4 step in the restart has been taken in all the numerical examples that follow.

Table 2. Extra time levels to be calculated and then disregarded for each filter.

The absolute stability regions of the filtered MS methods are found in the following manner. For example, the stability region that results from applying filter $P_7^{\left(-3\right)}$ with $N_0=6$ is found by applying the method to the test problem (3). Then the amplification factor $R\left(z\right)$ of the method where $z=k\lambda$ is found by looking at

$2y^1=\left(1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}\right)y^0(\; RK4\; )$

$y^2=\frac{\left(1+\frac{z}{3}\right)y^0+\frac{4}{3}z{y}^1}{\left(1-\frac{z}{3}\right)}=-\frac{y_0\left(z^5+4z^4+12z^3+24z^2+30z+18\right)}{6z-18}(\; MS\; )$

$⋮$

$\begin{array}{l}{y}^6=\frac{\left(1+\frac{z}{3}\right)y^4+\frac{4}{3}z{y}^5}{\left(1-\frac{z}{3}\right)}\\ =-\frac{y_0\left(65z^9+260z^8+954z^7+2256z^6+4147z^5+5754z^4+5976z^3+4428z^2+2106z+486\right)}{2z^5-30z^4+180z^3-540z^2+810z-486}\end{array}$

and

$\begin{array}{c}{y}^6=\frac{1}{64}\left(5y^0-18y^1+15y^2+20y^3-45y^4+30y^5+57y^6\right)\left(\text{apply}\text{filter}\right)\\ =\frac{K{y}_0}{192\left(z^6-18z^5+135z^4-540z^3+1215z^2-1458z+729\right)}\\ =R\left(z\right)y^0\mathrm{,}\end{array}$

where

$\begin{array}{c}K=16538z^{10}+82757z^9+312567z^8+823791z^7+1652835z^6+2550312z^5\\ +3015144z^4+2672352z^3+1702944z^2+699840z+139968.\end{array}$

The absolute stability region contains all $z\in \mathcal{C}$ for which $\left|R\left(z\right)\right|\le 1$.

Figure 1 shows the absolute stability regions of filters $P_7^{\left(0\right)}$, $P_7^{\left(-1\right)}$, $P_7^{\left(-2\right)}$, and $P_7^{\left(-3\right)}$ applied every $N_0$ steps for $N_0=6$, $N_0=10$, $N_0=15$, and $N_0=20$. The regions of the three forward biased filters are not shown but the stability regions of $P_7^{\left(3\right)}$ are very similar to that of $P_7^{\left(-3\right)}$, the regions for $P_7^{\left(2\right)}$ are similar to $P_7^{\left(-2\right)}$, and the regions of $P_7^{\mathrm{<mo\; stretchy="false">\; (\; </mo>1\; <mo\; stretchy="false">\; )\; </mo>}}$ are similar to $P_7^{\left(-1\right)}$. The trend that is observed is that increasing the frequency of filtering causes more of the left half-plane to be included and some of the imaginary axis coverage to be lost. This indicates that if the eigenvalues of the coefficient matrix (or of the Jacobian matrix if the system is nonlinear) of the ODE system are purely imaginary or have small negative real parts that only very infrequent filtering will be required for stability. However, if the eigenvalues have significant negative real parts more frequent filtering will be necessary.

In general the centered $P_7^{\left(0\right)}$ and marginally forward $P_7^{\left(1\right)}$ and backward $P_7^{\left(-1\right)}$ biased filters have the largest stability regions while the more extreme forward and backward biased filters such as $P_7^{\left(-3\right)}$ and $P_7^{\left(3\right)}$ have smaller stability regions. Filter $P_7^{\left(-3\right)}$ has a smaller stability region, but is attractive since it does not require any additional time levels to be computed. Figure 2 shows the stability regions of the four filters in Figure 1 for filtering every 6 and 20 steps in the same image.

Figure 1. Absolute stability regions in the complex plane of filtering every $N_0^{th}$ step from outer to inner: $N_0=6$, $N_0=10$, $N_0=15$, and $N_0=20$. Top left: $P_7^{\left(0\right)}$. Top right: $P_7^{\left(-1\right)}$. Bottom left: $P_7^{\left(-2\right)}$. Bottom right: $P_7^{\left(-3\right)}$.

Figure 2. Absolute stability regions in the complex plane of filtering every 6 steps (left) and 20 steps (right) listed from most to least negative real axis coverage: $P_7^{\left(-1\right)}$, $P_7^{\left(-2\right)}$, $P_7^{\left(0\right)}$, and $P_7^{\left(-3\right)}$.

6. Numerical Examples

In the numerical examples, first a nonlinear ODE is considered, and then linear ODE systems with coefficient matrices with real and purely imaginary eigenvalues are considered. Next, a linear variable coefficient ODE is solved. Finally, the filtered MS method is used in a method of lines approach to approximate the time derivative in linear PDEs which is where the filtered MS method is very competitive with standard methods for this type of problem.

6.1. Nonlinear ODE

The nonlinear initial value problem

$y^{\prime }=1-y^2,y\left(0\right)=0$ (29)

has the exact solution $y\left(t\right)=\tanh \left(t\right)$. For this equation $\frac{\partial f}{\partial y}=-2y<0$ and for any $y>0$ this value is out of the stability interval of the MS method and instability should be expected as is illustrated in Figure 3. Since the stability region needs substantial negative real axis coverage, the filters will have to be applied frequently. Table 3 lists the accuracy of the seven filtered MS methods over several time intervals with two different filtering frequencies, $N=5$ (except for $N=6$ with $P_7^{\left(-3\right)}$ ) and $N=15$. All seven methods produce accurate solutions with $N=5$ at the final time of $t=100$. However, $N=15$ is not frequent enough and the errors are larger as some instability begins to set in over long time intervals. If N is increased to 25, instability sets in and the solutions become infinite by all seven filtered MS methods.

6.2. Linear Constant Coefficient ODE Systems

The linear constant coefficient ODE system

$y^{\prime }=Ay$ (30)

is considered with two coefficient matrices, A. One has purely imaginary eigenvalues and one has complex eigenvalues with negative real parts.

Table 3. Errors from problem (29) with $k=0.125$ and $N_0=5$ (Columns 3 to 5) and $N_0=15$ (Columns 6 to 8).

Figure 3. Exact (dashed) and MS (solid) solution of problem (29) with $k=0.125$ up to time $t=30$ (left) and $t=100$ (right).

The matrix

$A=\left[\begin{array}{cc}0& 2\\ -2& 0\end{array}\right]$ (31)

has eigenvalues $\pm 2i$. The initial condition for the problem is $y\left(0\right)={\left(\mathrm{1,2}\right)}^{\text{T}}$ which leads to the exact solution

$y_1=\cos 2t+2\sin 2t$

$y_2=-\sin 2t+2\cos 2t.$

The eigenvalues can be scaled byk so that they lie in the stability region of the MS method and filtering is not necessary for stability. However, the seven filters are applied to verify that a fourth-order convergence rate is preserved when the filters are applied. In the left image of Figure 4, a convergence plot is shown for the number of time steps varying from 40 to 5120 on the interval $\left[\mathrm{0,8}\right]$ and the filters are applied every 25 steps. Once the convergence trend sets in, all lines have a slope of approximately four. While all seven filtered methods converge at the same rate, they do not reach the same solution at $t=100$. At this time, the unfiltered MS method has error 2.2e-10, the $P_7^{\left(-3\right)}$ method 3.2e-11, and the $P_7^{\left(3\right)}$ method 4.2e-10. This is typical when the MS method is applicable and can be compared to the filtered methods. A particular filter may give a solution slightly more or less accurate than the MS method but it is problem dependent as to which filter gives the most accurate solution.

The matrix

$A=\left[\begin{array}{cccc}-4& 2& 0& 2\\ -4& 1& 1& 2\\ -2& 1& -1& 2\\ -1& 1& -1& 0\end{array}\right]$ (32)

has eigenvalues $-1+1i$ and $-1-1i$ which are each of multiplicity two. The initial condition for the problem is $y\left(0\right)={\left(1,0,1,0\right)}^{\text{T}}$. The exact solution is

$y_1={\text{e}}^{-t}\cos t-3{\text{e}}^{-t}\sin t$

$y_2=-3{\text{e}}^{-t}\sin t$

Figure 4. Convergence plot of the filtered MS methods $P_7^{\left(-3\right)}$ (cyan), $P_7^{\left(-2\right)}$ (green), $P_7^{\left(-1\right)}$ (blue), $P_7^{\left(0\right)}$ (black), $P_7^{\left(1\right)}$ (magenta), $P_7^{\left(2\right)}$ (yellow), and $P_7^{\left(3\right)}$ (red). Left: purely imaginary eigenvalues (31) filtered every 25 steps. Right: complex eigenvalues with negative real part (32) filtered every 10 steps.

$y_3={\text{e}}^{-t}\cos t-2{\text{e}}^{-t}\sin t$

$y_4=-2{\text{e}}^{-t}\sin t.$

Since the eigenvalues have negative real parts a more frequent application of the filters is necessary to obtain enough left half-plane coverage of the stability region. In the right image of Figure 4 a convergence plot is shown for the number of time steps varying from 40 to 2560 on the interval $\left[\mathrm{0,8}\right]$ and the filters are applied every 10 steps. Once the convergence trend sets in, all lines have a slope of approximately four.

6.3. Linear Variable Coefficient Linear System

This example considers the second-order ODE IVP

$y^{″}+t{y}^{\prime }+y=0=0,y\left(0\right)=0,y^{\prime }\left(0\right)=1.$ (33)

The exact solution is $y\left(t\right)=\exp \left(-0.5t^2\right){\displaystyle {\int }_0^t}\exp \left(0.5x^2\right)\text{d}x$. For computational purposes, the second-order ODE is converted to an equivalent first-order system of the form (30) with coefficient matrix

$A=\left[\begin{array}{cc}0& 1\\ -1& -t\end{array}\right]$ (34)

and initial condition $y\left(0\right)={\left(\mathrm{0,1}\right)}^{\text{T}}$. At $t=0$ the coefficient matrix has purely imaginary eigenvalues but as the system is advanced in time the eigenvalues transition to being real and negative and grow in magnitude as time is advanced. Thus, overall this problem will need frequent filtering. Filters $P_7^{\left(-1\right)}$ and $P_7^{\left(0\right)}$ have stability regions with good negative real axis coverage so they are applied ever 5 steps. With a relatively large time step size of $k=0.1$ the solution is advanced to time $t=20$. At this time, the error from the $P_7^{\left(-1\right)}$ solution is 1.87e-4 and the error from the $P_7^{\left(0\right)}$ solution is 1.36e-4.

6.4. Linear PDEs, Method of Lines, and Advection

The advection equation

$\frac{\partial u}{\partial t}-\frac{\partial u}{\partial x}=0.$ (35)

is considered on the interval $\left[\mathrm{0,1}\right]$ with periodic boundary conditions. The initial condition is taken from the exact solution $u\left(x\mathrm{,}t\right)={\sin }^{40}\left(\pi \left(x+t\right)\right)$. The space derivative is discretized with the Fourier pseudospectral method [12] using $M=80$ points. The eigenvalues of the first-order Fourier psuedospectral differentiation matrix are purely imaginary and are given by $iw$ where $w=-\frac{M}{2}+1,\cdots ,\frac{M}{2}+1$. The eigenvalues lie in the interval $\left[-\pi i/\Delta x\mathrm{,}\pi i/\Delta x\right]$ on the imaginary axis. The MS method is applicable without filtering. This is another example of where infrequent filtering may slightly improve the accuracy of the solution.

Table 4 gives the accuracy of the centered and backward biased filters up to time $t=1000$ with solutions filtered every 10 steps and every 100 steps. The less frequently filtered solutions are slightly more accurate but the heavier filtered solutions retain good accuracy as well. The more computationally expensive forward biased filters were not considered in this example as they have not shown to have any favorable properties when compared to the less expensive centered and backward biased filters. For comparison, the more computationally expensive explicit RK4 method which is commonly used on this type of problem is compared as well. The RK4 method was less accurate than the MS method or the three filtered MS methods.

The heat equation

$u_t=u_{xx}$ (36)

is considered on the domain $\left[\mathrm{0,1}\right]$ with initial condition $u\left(x\mathrm{,0}\right)=\sin \left(\pi x\right)$. Zero Dirichlet boundary conditions of $u\left(0\right)=0$ and $u\left(1\right)=0$ are applied. The exact solution is $u\left(x,t\right)=\exp \left(-{\pi }^{2}t\right)\sin \left(\pi x\right)$. The space derivative is discretized with the Chebyshev pseudospectral method [12] on the Chebyshev-Gauss-Lobbato (GCL) points which cluster densely around the boundary points and are less dense in the interior of the interval. The boundary conditions are enforced by setting the elements in the first and last rows of the derivative matrix to zero. This is a well-known “stiff” problem with large negative real eigenvalues. There are specifically designed numerical methods [7] that are appropriate for this type of problem and the Milne-Simpson method with its interval of absolute stability on the imaginary axis is certainly not one of them. However, with the use of filters it can nevertheless be applied.

The space interval is discretized with 20 CGL points and the $P_7^{\left(-3\right)}$ filtered MS method with a time step size of $k=0.00005$ is applied. The solution is advanced to time $t=0.4$, and the numerical and exact solution differ by 7.5e-14. The scaled eigenvalues and the absolute stability region of the MS method that is filtered every 6 time steps by $P_7^{\left(-3\right)}$ is shown in Figure 5.

Table 4. Errors from problem (35) with $k=0.001$ and $N_0=10$ (Columns 3 to 5) and $N_0=100$ (Columns 6 to 8).

Figure 5. Stability region of the MS method filtered every 6 steps by $P_7^{\left(-3\right)}$ and the scaled eigenvalues of the space discretization of Equation (36).

7. Conclusions

The Milne-Simpson method has fourth-order convergence which is the highest order theoretically possible for a two-step LMM. However, the method is only weakly stable, has only an interval of absolute convergence on the imaginary axis, and it is implicit. These factors lead to the MS method seldom being used in applications.

Seven-point filters which retain fourth-order accuracy can be used to expand the size of and to manipulate the location of the absolute stability region of the method. Systems with coefficient matrices or Jacobian matrices with eigenvalues with negative real parts require more frequent filtering in the range of every 5 to 20-time steps. Even if the eigenvalues are purely imaginary and lie within the stability region of the unfiltered MS method, infrequent filtering every 100 to 200 stepping can be beneficial to accuracy and stability over long time intervals. The enhanced stability properties of the filtered MS method make it applicable to a much larger class of problems than the unfiltered MS method is applicable to. The filtered MS method is in particular an attractive time integrator in a method of lines approach with linear PDEs where it is more efficient than fourth-order Runge-Kutta methods that are very popular for this type of problem.

All seven of the filters can be effective, but we are unable to label which filter is definitively the best. Filters $P_7^{\left(1\right)}$, $P_7^{\left(2\right)}$, and $P_7^{\left(3\right)}$ require 4, 5, and 6 extra time levels respectively to be computed and disregarded and show no advantage over the other four filters. Due to their computational expense, they can not be recommended for applications. Filters $P_7^{\left(0\right)}$ and $P_7^{\left(-1\right)}$ allow for the largest stability domain coverage in the left half-plain and are well suited to problems involving eigenvalues with large negative real parts. Filter $P_7^{\left(-3\right)}$ is the only one of the filters that does not require additional time levels, but its application does not expand the stability region as much as the other filters and as a result, will require a smaller time-step size for stability than some of the other filtered MS methods.

The filtering technique that has been used in this work was previously used with the leapfrog (midpoint) method [4]. Both the leapfrog and Milne-Simpson methods are derived using polynomial basis functions. In particular, the leapfrog method is derived by discretizing the time derivative in (1) by a three-point centered polynomial approximation. Another way to approximate the derivative in (1) is by a three-point centered approximation using Radial Basis Functions (RBFs) [13]. The resulting method based on RBFs is also weakly stable. Our future work involves applying the filtering technique to methods based on RBFs.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Dahlquist, G. (1985) 33 Years of Numerical Instability, Part I. BIT Numerical Mathematics, 25, 188-204. https://doi.org/10.1007/BF01934997
[2]	Iri, M. (1963) A Stabilizing Device for Unstable Numerical Solutions of Ordinary Differential Equations—Design Principle and Application of a Filter. Journal of the Information Processing Society of Japan, 4, 249-260.
[3]	Aluthge, A. (1985) Filtering and Extrapolation Techniques in the Numerical Solution of Ordinary Differential Equations. Master’s Thesis, University of Ottawa, Ottawa, 1-70.
[4]	Aluthge, A., Sarra, S.A. and Estep, R. (2016) Filtered Leapfrog Time Integration with Enhanced Stability Properties. Journal of Applied Mathematics and Physics, 4, 1354-1370. https://doi.org/10.4236/jamp.2016.47145
[5]	Butcher, J. (2003) Numerical Methods for Ordinary Differential Equations. Wiley, Hoboken. https://doi.org/10.1002/0470868279
[6]	Hairer, E., Norsett, S. and Wanner, G. (2000) Solving Ordinary Differential Equations I: Nonstiff Problems. Springer, Berlin.
[7]	Hairer, E. and Wanner, G. (2000) Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems. 2nd Edition, Springer, Berlin.
[8]	Henrici, P. (1962) Discrete Variable Methods in Ordinary Differential Equations. Wiley, Hoboken.
[9]	Iserles, A. (1973) A First Course in the Analysis of Differential Equations. Wiley, Hoboken.
[10]	Lambert, J. (1973) Computational Methods in Ordinary Differential Equations. Wiley, Hoboken.
[11]	Gragg, W.B. (1963) Re-Peated Extrapolation to the Limit in the Numerical Solution of Ordinary Differential Equations. Ph.D. Thesis, University of California, Los Angeles (UCLA), Los Angeles, 1-111.
[12]	Canuto, C., Hussaini, M.Y., Quarteroni, A. and Zang, T.A. (1988) Spectral Methods in Fluid Dynamics. Springer, Berlin. https://doi.org/10.1007/978-3-642-84108-8
[13]	Sarra, S.A. and Kansa, E.J. (2009) Multiquadric Radial Basis Function Approximation Methods for the Numerical Solution of Partial Differential Equations. Advances in Computational Mechanics, 2, 1-199.