info

Internal stability of Runge-Kutta methods

Some experiments in NodePy

Note: this post was generated from an iPython notebook. You can download the notebook from github and execute all the code yourself.

Internal stability deals with the growth of errors (such as roundoff) introduced at the Runge-Kutta stages during a single Runge-Kutta step. It is usually important only for methods with a large number of stages, since that is when the internal amplification factors can be large. An excellent explanation of internal stability is given in this paper. Here we demonstrate some tools for studying internal stability in NodePy.

First, let’s load a couple of RK methods:

from nodepy import rk
reload(rk)
rk4 = rk.loadRKM('RK44')
ssprk4 = rk.loadRKM('SSP104')
print rk4
print ssprk4

Classical RK4
The original four-stage, fourth-order method of Kutta
 0   |
 1/2 |  1/2
 1/2 |  0    1/2
 1   |  0    0    1
_____|____________________
     |  1/6  1/3  1/3  1/6
SSPRK(10,4)
The optimal ten-stage, fourth order SSP Runge-Kutta method
 0   |
 1/6 |  1/6
 1/3 |  1/6   1/6
 1/2 |  1/6   1/6   1/6
 2/3 |  1/6   1/6   1/6   1/6
 1/3 |  1/15  1/15  1/15  1/15  1/15
 1/2 |  1/15  1/15  1/15  1/15  1/15  1/6
 2/3 |  1/15  1/15  1/15  1/15  1/15  1/6   1/6
 5/6 |  1/15  1/15  1/15  1/15  1/15  1/6   1/6   1/6
 1   |  1/15  1/15  1/15  1/15  1/15  1/6   1/6   1/6   1/6
_____|____________________________________________________________
     |  1/10  1/10  1/10  1/10  1/10  1/10  1/10  1/10  1/10  1/10

Absolute stability regions

First we can use NodePy to plot the region of absolute stability for each method. The absolute stability region is the set

\(\\{ z \in C : |\phi (z)|\le 1 \\}\)

where \(\phi(z)\) is the stability function of the method:

\(1 + z b^T (I-zA)^{-1}\)

If we solve \(u'(t) = \lambda u\) with a given method, then \(z=\lambda \Delta t\) must lie inside this region or the computation will be unstable.

p,q = rk4.stability_function()
print p
h1=rk4.plot_stability_region()

         4          3       2
0.04167 x + 0.1667 x + 0.5 x + 1 x + 1

p,q = ssprk4.stability_function()
print p
h2=ssprk4.plot_stability_region()

           10             9            8             7           6
3.969e-09 x  + 2.381e-07 x + 6.43e-06 x + 0.0001029 x + 0.00108 x
            5           4          3       2
 + 0.00787 x + 0.04167 x + 0.1667 x + 0.5 x + 1 x + 1

Internal stability

The stability function tells us by how much errors from one step are amplified in the next one. This is important since we introduce truncation errors at every step. However, we also introduce roundoff errors at the each stage within a step. Internal stability tells us about the growth of those. Internal stability is typically less important than (step-by-step) absolute stability for two reasons:

  • Roundoff errors are typically much smaller than truncation errors, so moderate amplification of them typically is not significant
  • Although the propagation of stage errors within a step is governed by internal stability functions, in later steps these errors are propagated according to the (principal) stability function

Nevertheless, in methods with many stages, internal stability can play a key role.

Questions: In the solution of PDEs, large spatial truncation errors enter at each stage. Does this mean internal stability becomes more significant? How does this relate to stiff accuracy analysis and order reduction?

Internal stability functions

We can write the equations of a Runge-Kutta method compactly as

\(y = u^n e + h A F(y)\)
\(u^{n+1} = u^n + h b^T F(y),\)

where \(y\) is the vector of stage values, \(u^n\) is the previous step solution, \(e\) is a vector with all entries equal to 1, \(h\) is the step size, \(A\) and \(b\) are the coefficients in the Butcher tableau, and \(F(y)\) is the vector of stage derivatives. In floating point arithmetic, roundoff errors will be made at each stage. Representing these errors by a vector \(r\), we have

\(y = u^n e + h A F(y) + r.\)

Considering the test problem \(F(y)=\lambda y\) and solving for \(y\) gives

\(y = u^n (I-zA)^{-1}e + (I-zA)^{-1}r,\)

where \(z=h\lambda\). Substituting this result in the equation for \(u^{n+1}\) gives

\(u^{n+1} = u^n (1 + zb^T(I-zA)^{-1}e) + zb^T(I-zA)^{-1}r = \psi(z) u^n + \theta(z)^T r.\)

Here \(\psi(z)\) is the stability function of the method, that we already encountered above. Meanwhile, the vector \(\theta(z)\) contains the internal stability functions that govern the amplification of roundoff errors \(r\) within a step:

\(\theta(z) = z b^T (I-zA)^{-1}.\)

Let’s compute \(\theta\) for the classical RK4 method:

theta=rk4.internal_stability_polynomials()
theta

    [poly1d([1/24, 1/12, 1/6, 1/6, 0], dtype=object),
     poly1d([1/12, 1/6, 1/3, 0], dtype=object),
     poly1d([1/6, 1/3, 0], dtype=object),
     poly1d([1/6, 0], dtype=object)]

for theta_j in theta:
    print theta_j

         4           3          2
0.04167 x + 0.08333 x + 0.1667 x + 0.1667 x
         3          2
0.08333 x + 0.1667 x + 0.3333 x
        2
0.1667 x + 0.3333 x
 
0.1667 x

Thus the roundoff errors in the first stage are amplified by a factor \(z^4/24 + z^3/12 + z^2/6 + z/6\), while the errors in the last stage are amplified by a factor \(z/6\).

Internal instability

Usually internal stability is unimportant since it relates to amplification of roundoff errors, which are very small. Let’s think about when things can go wrong in terms of internal instability. If \(|\theta(z)|\) is of the order \(1/\epsilon_{machine}\), then roundoff errors could be amplified so much that they destroy the accuracy of the computation. More specifically, we should be concerned if \(|\theta(z)|\) is of the order \(tol/\epsilon_{machine}\) where \(tol\) is our desired error tolerance. Of course, we only care about values of \(z\) that lie inside the absolute stability region \(S\), since internal stability won’t matter if the computation is not absolutely stable.

We can get some idea about the amplification of stage errors by plotting the curves \(|\theta(z)|=1\) along with the stability region. Ideally these curves will all lie outside the stability region, so that all stage errors are damped.

rk4.internal_stability_plot()


ssprk4.internal_stability_plot()

For both methods, we see that some of the curves intersect the absolute stability region, so some stage errors are amplified. But by how much? We’d really like to know the maximum amplification of the stage errors under the condition of absolute stability. We therefore define the maximum internal amplification factor \(M\):

\(M = \max_j \max_{z \in S} |\theta_j(z)|\) 
print rk4.maximum_internal_amplification()
print ssprk4.maximum_internal_amplification()

2.15239281554
4.04399941143

We see that both methods have small internal amplification factors, so internal stability is not a concern in either case. This is not surprising for the method with only four stages; it is a surprisingly good property of the method with ten stages.

Questions: Do SSP RK methods always (necessarily) have small amplification factors? Can we prove it?

Now let’s look at some methods with many stages.

Runge-Kutta Chebyshev methods

The paper of Verwer, Hundsdorfer, and Sommeijer deals with RKC methods, which can have very many stages. The construction of these methods is implemented in NodePy, so let’s take a look at them. The functions RKC1(s) and RKC2(s) construct RKC methods of order 1 and 2, respectively, with \(s\) stages.

s=4
rkc = rk.RKC1(s)
print rkc

RKC41

 0    |
 1/16 |  1/16
 1/4  |  1/8   1/8
 9/16 |  3/16  1/4   1/8
______|________________________
      |   1/4   3/8   1/4   1/8
rkc.internal_stability_plot()

It looks like there could be some significant internal amplification here. Let’s see:

rkc.maximum_internal_amplification()

    11.760869405962685

Nothing catastrophic. Let’s try a larger value of \(s\):

s=20
rkc = rk.RKC1(s)
rkc.maximum_internal_amplification()

    42.665327220219126

As promised, these methods seem to have good internal stability properties. What about the second-order methods?

s=20
rkc = rk.RKC2(s)
rkc.maximum_internal_amplification()

    106.69110992619214

Again, nothing catastrophic. We could take \(s\) much larger than 20, but the calculations get to be rather slow (in Python) and since we’re using floating point arithmetic, the accuracy deteriorates.

Remark: we could do the calculations in exact arithmetic using Sympy, but things would get even slower. Perhaps there are some optimizations that could be done to speed this up. Or perhaps we should use Mathematica if we need to do this kind of thing.

Remark 2: of course, for the RKC methods the internal stability polynomials are shifted Chebyshev polynomials. So we could evaluate them directly in a stable manner using the three-term recurrence (or perhaps scipy’s special functions library). This would also be a nice check on the calculations above.

Other methods with many stages

Three other classes of methods with many stages have been implemented in NodePy:

  • SSP families
  • Integral deferred correction (IDC) methods
  • Extrapolation methods

SSP Families

s=20
ssprk = rk.SSPRK2(s)
ssprk.internal_stability_plot()
ssprk.maximum_internal_amplification()

    2.0212921484995547


s=25 # # of stages
ssprk = rk.SSPRK3(s)
ssprk.internal_stability_plot()
ssprk.maximum_internal_amplification()

    3.8049237837215397

The SSP methods seem to have excellent internal stability properties.

IDC methods

p=6 #order
idc = rk.DC(p-1)
print len(idc)
idc.internal_stability_plot()
idc.maximum_internal_amplification()

26
    6.4140166271998815

IDC methods also seem to have excellent internal stability.

Extrapolation methods

p=6 #order
ex = rk.extrap(p)
print len(ex)
ex.internal_stability_plot()
ex.maximum_internal_amplification()

16
6

Not so good. Let’s try a method with even more stages (this next computation will take a while; go stretch your legs).

p=10 #order
ex = rk.extrap(p)
print len(ex)
ex.maximum_internal_amplification()

46
    28073.244376758907

Now we’re starting to see something that might cause trouble, especially since such high order extrapolation methods are usually used when extremely tight error tolerances are required. Internal amplification will cause a loss of about 5 digits of accuracy here, so the best we can hope for is about 10 digits of accuracy in double precision. Higher order extrapolation methods will make things even worse. How large are their amplification factors? (Really long calculation here…)

pmax = 12
ampfac = np.zeros(pmax+1)
for p in range(1,pmax+1):
    ex = rk.extrap(p)
    ampfac[p] = ex.maximum_internal_amplification()
    print p, ampfac[p]

1 1.99777378912
2 2.40329384375
3
 5.07204078733
4
 17.747335803
5
 69.62805786
6
 97.6097450835
7
 346.277441462
8
 1467.40356089
9
 6344.16303534
10
 28073.2443768
11
 126011.586473
12
 169897.662582
    []


semilogy(ampfac,linewidth=3)

    []

We see roughly geometric growth of the internal amplification factor as a function of the order \(p\). It seems clear that very high order extrapolation methods applied to problems with high accuracy requirements will fall victim to internal stability issues.

Posted on .

previous next


history



, , , ,