david/A formula for the numerator of the stability function for SDIRKs

Let \(\psi(z)=P(z)/Q(z)\) be the stability function of a diagonally implicit Runge–Kutta method with coefficients \((A,b)\). Here we derive a formula for the numerator \(P(z)\) as a power series in \(1+z/r\) in order to analyze its radius of absolute monotonicity.

We start with the form of the stability function given in equation (4.16) of GKS:

\[\psi(z) = v_{m+1} + (\alpha_{m+1} + z\beta_{m+1})(I-\alpha-z\beta)^{-1} v.\]

Let \(r\le R(a,b)\); then we can take \(\beta = \alpha/r\) with \(v\) and \(\alpha\) non-negative entrywise:

\[\psi(z) = v_{m+1} + (1 + z/r)\alpha_{m+1}(I-(1+z/r)\alpha)^{-1} v.\]

Here \(\alpha\) and \(v\) are the upper parts and \(\alpha_{m+1},v_{m+1}\) are the last row/entry.

Since the method is singly diagonally implicit, we can write \(\alpha = \gamma I + N\) where \(\gamma\ge0\) and \(N\) is nilpotent of degree \(s\) (i.e., \(N^s=0\)). Then we can simplify the matrix inverse appearing above as follows:

\[\begin{align*} (I-(1+z/r)\alpha)^{-1} & = (I-(1+z/r)(\gamma I +N))^{-1} \\ & = ((1-\gamma(1+z/r))I - (1+z/r)N)^{-1} \\ & = \frac{1}{1-\gamma(1+z/r)}(I-\frac{1+z/r}{1-\gamma(1+z/r)} N)^{-1} \\ & = \sum_{j=0}^{s-1} \frac{(1+z/r)^j}{(1-\gamma(1+z/r))^{j+1}} N^j \\ & = \frac{1}{(1-\gamma(1+z/r))^s} \sum_{j=0}^{s-1} (1+z/r)^j (1-\gamma(1+z/r))^{s-j-1} N^j.\end{align*}\]

Thus \[P(z) = (1-\gamma(1+z/r))^s v_{m+1} + \sum_{j=0}^{s-1} (1+z/r)^{j+1} (1-\gamma(1+z/r))^{s-j-1} \alpha_{m+1} N^j v.\]

Next we use the binomial expansion for both terms involving \(\gamma\) to obtain

\[P(z) = v_{m+1} \sum_{l=0}^s {s\choose l}(-\gamma)^l(1+z/r)^l + \sum_{j=1}^s \alpha_{m+1} N^{j-1} v \sum_{l=0}^{s-j} {s-j\choose l} (-\gamma)^l (1+z/r)^{j+l}.\]

Now everything is written in terms of powers of \((1+z/r)^l\); we just need to collect powers. We first change variables in the second sum to obtain

\[P(z) = v_{m+1} + v_{m+1} \sum_{l=1}^s {s\choose l}(-\gamma)^l(1+z/r)^l + \sum_{j=1}^s \alpha_{m+1} N^{j-1} v \sum_{l=j}^{s} {s-j\choose l-j} (-\gamma)^{l-j} (1+z/r)^l.\]

Finally, reversing the order of summation in the double sum yields:

\[P(z) = v_{m+1} + v_{m+1} \sum_{l=1}^s {s\choose l}(-\gamma)^l(1+z/r)^l + \sum_{l=1}^s \sum_{j=1}^{l}\alpha_{m+1} N^{j-1} v {s-j\choose l-j} (-\gamma)^{l-j} (1+z/r)^l,\]

and collecting coefficients we find

\[P(z) = v_{m+1} + \sum_{l=1}^s (1+z/r)^l\left({s\choose l}(-\gamma)^l v_{m+1} + \sum_{j=1}^{l}\alpha_{m+1} N^{j-1} v {s-j\choose l-j} (-\gamma)^{l-j} \right).\]

This is the final form; now we would like to show that the quantity in large parentheses is non-negative (for each \(l\) when \(0\le r\le R(A,b)\).