Stability of Taylor-Based ODE Integrators

Why Stability Matters

A numerical solution to an initial-value problem is produced by repeatedly taking steps. Even when each step is computed correctly, small errors can enter through truncation, rounding, or imperfect starting data. Stability asks whether the numerical method controls those errors as the steps accumulate.

For a Taylor-based method, the step size affects both the approximation made at each step and the way errors are propagated. These are related but distinct questions. A step can satisfy a method's stability condition and still be too large for an accurate approximation.

A standard way to study this distinction is to apply the method to the scalar test equation $y'=\lambda y$, where $\lambda$ may be real or complex. The exact solution is $y(t)=y(0)e^{\lambda t}$. If $\operatorname{Re}(\lambda)<0$, its magnitude decays as time increases.

Taylor Methods and the Amplification Factor

For a positive integer $p$, the degree-$p$ Taylor method advances an approximation from $t_n$ to $t_{n+1}=t_n+h$ by retaining terms through degree $p$ in the Taylor expansion about $t_n$. On $y'=\lambda y$, repeated differentiation gives $y^{(k)}=\lambda^k y$. Consequently, the update is multiplication by a polynomial in $z=h\lambda$.

Definition: For the degree-$p$ Taylor method applied to $y'=\lambda y$, the amplification factor is $$ R_p(z)=\sum_{k=0}^{p}\frac{z^k}{k!},\qquad z=h\lambda. $$ One numerical step is $y_{n+1}=R_p(z)y_n$.

Here $h$ has units of time and $\lambda$ has units of inverse time, so $z$ is dimensionless. After $n$ steps with constant $h$ and $\lambda$, the method gives $y_n=R_p(z)^n y_0$. In contrast, the exact solution over the same elapsed time has amplification $e^{nz}$.

Definition: The absolute stability region of a method is the set of values $z$ for which $|R_p(z)|\leq 1$. Values with $|R_p(z)|<1$ damp perturbations from step to step; values with $|R_p(z)|>1$ amplify them.

The boundary case $|R_p(z)|=1$ is not damped: a perturbation keeps the same magnitude under repeated steps. The stability region is a property of the method as a function of $z$, not of $h$ alone. A step-size restriction arises after $\lambda$ is specified.

Why the Amplification Criterion Works

The stability test follows directly from how the method updates a perturbed value. Suppose the computed value at step $n$ differs from a reference numerical value by $\delta_n$. Since the update is multiplication by the same factor for both values, their difference is also multiplied by that factor.

Apply the update to both values.
If the reference value is $y_n$, the perturbed value is $y_n+\delta_n$. Their next-step values are $R_p(z)y_n$ and $R_p(z)(y_n+\delta_n)$.

Subtract the updates.
The new difference is $\delta_{n+1}=R_p(z)(y_n+\delta_n)-R_p(z)y_n=R_p(z)\delta_n$.

Repeat for $n$ steps.
Iteration gives $\delta_n=R_p(z)^n\delta_0$, and hence $|\delta_n|=|R_p(z)|^n|\delta_0|$.

Therefore, if $|R_p(z)|>1$, even a small initial perturbation grows geometrically in this test problem. If $|R_p(z)|<1$, it decays. This derivation explains the amplification-factor criterion; it does not claim that every source of error in every ODE behaves exactly like this scalar test equation.

Real Negative Values and Step-Size Restrictions

For a real negative $\lambda$, the exact solution decays. Write $z=h\lambda$, so a positive step $h$ gives $z<0$. The stability condition can then be solved along the negative real axis.

For the first-degree Taylor method, $R_1(z)=1+z$. The condition $|1+z|\leq 1$ is equivalent to $-1\leq 1+z\leq 1$, or $-2\leq z\leq 0$. If $\lambda<0$, this means $h\leq 2/|\lambda|$.

For the second-degree method, $R_2(z)=1+z+z^2/2$. For real $z$, this polynomial is positive because its discriminant is $1-2=-1<0$ and its leading coefficient is positive. Thus the stability condition reduces to $R_2(z)\leq 1$. Subtracting one gives $z+z^2/2=z(1+z/2)\leq 0$, which holds for $-2\leq z\leq 0$. This method has the same negative-real-axis step bound as the first-degree method.

The third-degree method has $R_3(z)=1+z+z^2/2+z^3/6$. Its negative-axis stability endpoint is not $-2$. Since $R_3'(z)=1+z+z^2/2>0$ for all real $z$, $R_3$ is strictly increasing. Its stable interval on the negative real axis runs from the unique solution of $R_3(z)=-1$ up to zero. That endpoint is approximately $-2.512745$.

Worked Example: A Step-Size Bound for the First-Degree Method

Consider $y'=-4y$, with $\lambda=-4\ \text{s}^{-1}$. For the first-degree Taylor method, stability requires $-2\leq h\lambda\leq 0$. Substituting $\lambda=-4$ gives $-2\leq -4h\leq 0$. Dividing by $-4$ reverses the inequalities, so $0\leq h\leq 0.5\ \text{s}$.

Take $h=0.6\ \text{s}$, which exceeds the bound. Then $z=h\lambda=(0.6\ \text{s})(-4\ \text{s}^{-1})=-2.4$, and $R_1(z)=1-2.4=-1.4$. Starting from $y_0=1$, one step gives $y_1=-1.4$; a second gives $y_2=(-1.4)(-1.4)=1.96$. Indeed, the two-step amplification is $R_1(z)^2=(-1.4)^2=1.96>1$, so the perturbation grows. The exact solution's two-step amplification over $1.2\ \text{s}$ is $e^{(-4)(1.2)}=e^{-4.8}<1$. Substitution into the update confirms that this numerical growth is caused by the chosen step, not by growth in the exact solution.

Worked Example: A Stable but Inaccurate Second-Degree Step

Use the second-degree Taylor method on $y'=-3y$, starting with $y_0=2$, and take $h=0.5\ \text{s}$. Here $\lambda=-3\ \text{s}^{-1}$, so $z=h\lambda=-1.5$. The negative-axis stability interval is $[-2,0]$, and $z=-1.5$ lies inside it.

The amplification factor is $R_2(-1.5)=1-1.5+(-1.5)^2/2=1-1.5+2.25/2=0.625=5/8$. Thus $y_1=(5/8)(2)=1.25$. Checking the update directly gives $1+z+z^2/2=1-1.5+1.125=0.625$, and $0.625(2)=1.25$, as required. The exact value after this step is $y(0.5)=2e^{-1.5}$, approximately $0.4463$. The numerical value is far from the exact value, even though the step is stable. Stability controls amplification of perturbations; it does not certify that this step is accurate.

Worked Example: A Third-Degree Step Outside the Stable Interval

For $y'=-y$, take $h=3\ \text{s}$. Then $z=-3$, which lies below the third-degree method's negative-axis stability endpoint, approximately $-2.512745$. The factor is

R_3(-3)=1-3+\frac{9}{2}-\frac{27}{6} =1-3+4.5-4.5=-2.

With $y_0=1$, the first computed value is $y_1=-2$, and the second is $y_2=(-2)(-2)=4$. Direct substitution into the recurrence confirms the two-step factor $R_3(-3)^2=(-2)^2=4$. The exact solution at $t=6\ \text{s}$ is $e^{-6}$, which is positive and less than one. The growing numerical values therefore signal instability, not growth in the exact solution.

Complex Values of the Test Parameter

If $\lambda$ is complex, stability still depends on the modulus of the amplification factor. For example, a negative real part of $\lambda$ makes the exact solution decay in magnitude, but it does not by itself guarantee that a particular step is stable. The value $z=h\lambda$ must lie in the method's stability region.

Worked Example: Testing a Complex Second-Degree Step

Let $\lambda=(-1+2i)\ \text{s}^{-1}$ and $h=0.5\ \text{s}$. Then $z=-0.5+i$. For the second-degree Taylor method,

$$ z^2=(-0.5+i)^2=0.25-i-1=-0.75-i, $$

R_2(z)=1+z+\frac{z^2}{2} =1+(-0.5+i)+(-0.375-0.5i) =0.125+0.5i.

Its modulus satisfies $|R_2(z)|^2=0.125^2+0.5^2=0.015625+0.25=17/64$, so $|R_2(z)|=\sqrt{17}/8$, approximately $0.515$, which is less than one. The step is stable. The exact one-step amplification has modulus $|e^z|=e^{\operatorname{Re}(z)}=e^{-0.5}$, approximately $0.607$. The numerical update damps this mode, although its damping is not exactly the same as the exact damping.

Stability Is Not Accuracy

For the test equation, the exact one-step amplification is $e^z$, while a degree-$p$ Taylor method uses $R_p(z)$. The Taylor polynomial agrees with the exponential through degree $p$, so $R_p(z)-e^z$ is of order $z^{p+1}$ as $z$ approaches zero. This describes the local approximation when the scaled step is small; it does not imply that every stable value of $z$ gives a small error.

A practical distinction is:

Stability: Does the numerical update keep perturbations from growing uncontrollably?
Accuracy: Is the numerical value sufficiently close to the exact solution for the purpose at hand?

A step-size choice must account for both. A stable step may have a large truncation error, as the second-degree example showed. Conversely, a method may have a small local Taylor error for a step but still amplify perturbations if its amplification factor has modulus greater than one. Error estimates and adaptive step selection address accuracy; amplification analysis identifies stability restrictions.

Common pitfall: A decaying exact solution does not guarantee that the numerical solution decays. For $y'=-4y$, the first-degree method with $h=0.6\ \text{s}$ has $z=-2.4$ and $R_1(z)=-1.4$. Since $|R_1(z)|=1.4>1$, its computed perturbations grow, even though $\lambda<0$. The correction is to test $|R_p(h\lambda)|$, not to infer numerical stability from the sign of $\lambda$ alone.

The test equation gives a clear, computable diagnostic, particularly when a problem contains decaying modes with different rates. A mode with a large negative $\lambda$ can impose a more restrictive step bound than a slowly decaying mode. The method's stability region and the values $h\lambda$ together determine whether the chosen step controls those modes.

Check Your Understanding

Use the amplification-factor definition and show the relevant calculations in each case.

For the first-degree Taylor method applied to $y'=-5y$, what is the largest step size allowed by the negative-real-axis stability condition?
For the second-degree Taylor method with $z=-1$, calculate $R_2(z)$. Is the step stable?
For the first-degree method with $z=-2.2$, calculate the one-step amplification and decide whether perturbations grow or decay.
Explain why a stable step need not give an accurate approximation to the exact solution.
For a complex value of $z$, which quantity must be compared with one to test absolute stability?