Internal stability of optimal 2nd order explicit SSP Runge-Kutta methods

The \(s\)-stage method is

\[\begin{align} y_1 & = u^n \\ y_j & = y_{j-1} + \frac{h}{s-1} f(y_{j-1}) & 2\le j\le s \\ u^{n+1} & = \frac{1}{s} u^n + \frac{s-1}{s} \left(y_s + \frac{h}{s-1} f(y_s) \right). \end{align}\]

It is convenient to define \[\nu = 1 + \frac{z}{s-1}.\] Then the stability function is \[ P(z) = \frac{1}{s} + \frac{s-1}{s} \nu^s \] and the internal stability functions are \[Q_{sj} = \frac{s-1}{s} \nu^{s-j+1}.\]

We would like to know the maximum amplification factor; i.e. the maximum of \(|Q_{sj}(z)|\) for \(z\) such that \(|P(z)|\le 1\). The latter condition is equivalent to \[\left| \nu^s - \frac{-1}{s-1}\right|\le \frac{s}{s-1}\] i.e., \(\nu^s\) lies in a disk of radius \(s/(s-1)\) centered at \(-1/(s-1)\). Thus \[|\nu^s|\le \frac{s+1}{s-1}.\] Thus (considering only the \(|\nu|>1\) case, since \(|\nu|\le 1\) is trivial) \[\begin{align} |Q_{sj}| & = \frac{s-1}{s} |\nu|^{s-j+1} \\ & \le \frac{s-1}{s}|\nu|^{s-1} \\ & \le \frac{s-1}{s} \left(\frac{s+1}{s-1}\right)^{\frac{s-1}{s}} \\ & < \frac{s+1}{s}. \end{align}\]

\(\left| \frac{(s-1) \nu^s}{s}+\frac{1}{s}\right| =1\) stability region boundaries for \(2\le s\le 7\)

I have checked this against results computed in NodePy and it appears to be correct. A special case is \(s=2\), when the bound above is actually smaller than 1, but the maximum amplification factor occurs for \(j=1\) and is equal to unity.

Internal stability of optimal 3rd order explicit SSP Runge-Kutta methods

The Shu–Osher coefficients for the optimal 3rd order SSPRK with \(s=n^2\) stages are given by \[\begin{align*} \alpha_{i+1,i} &= \begin{cases} \dfrac{n-1}{2n-1}, & \mbox{if } i = \frac{n(n+1)}{2} \\\\ 1,& \mbox{otherwise,} \end{cases} \\ \alpha_{\frac{n(n+1)}{2}+1,\frac{(n-1)(n-2)}{2}+1} &= \frac{n}{2n-1}, \\ \beta_{i+1,i} &= \frac{\alpha_{i+1,i}}{n^2-n}. \end{align*}\] Let \(k = \frac{n(n+1)}{2}+1\) and \(m= \frac{(n-1)(n-2)}{2}+1\), then the method can be written as follows: \[\begin{align*} y_1 &= u^n \\ y_j &= y_{j-1} + \frac{h}{n^2-n} f(y_{j-1}) \quad 2 \leq j\leq s, j \neq k \\ y_k &= \frac{n-1}{2n-1}y_{k-1} + \frac{n}{2n-1}y_m + h\frac{n-1}{(2n-1)(n^2-n)} f(y_{k-1}) \\ u^{n+1} &= y_{s+1} = y_s + \frac{h}{n^2-n} f(y_{s}). \end{align*}\]

When applying the method to the linear problem \(u'(t) = \lambda u(t)\), it is convenient to define \[\nu = 1 + \frac{z}{n^2-n},\] where \(z=\lambda h\). Assuming that \(r_j\) are perturbations of the \(j\)-stage, \(2 \leq j \leq s\), then the stage’s approximations become \[\begin{align*} y_j &= \nu y_{j-1} + r_j \quad 2 \leq j\leq s, j \neq k \\ y_k &= \nu \frac{n-1}{2n-1}y_{k-1} + \frac{n}{2n-1}y_m + r_k. \end{align*}\] Then the solution at the \(n+1\) step is \[\begin{align*} u^{n+1} &= \nu^{s-k+1}y_k+ \sum_{j=1}^{s-k} \nu^j r_{s-j+1} \\ &= \frac{n-1}{2n-1} \nu^{s-k+2} y_{k-1} + \frac{n}{2n-1} \nu^{s-k+1} y_m + \sum_{j=1}^{s-k+1} \nu^j r_{s-j+1}. \end{align*}\] For \(n \geq 2\) we have \(k \geq m\), thus \[\begin{align*} y_m = \nu^{m-1}u^n + \sum_{j=1}^{m-1} \nu^{j-1} r_{m-j+1} \end{align*}\] and \[\begin{align*} y_{k-1} = \nu^{k-2}u^n + \sum_{j=1}^{k-2} \nu^{j-1} r_{k-j}. \end{align*}\]

Substituting \(u_m\) and \(u_{k-1}\) to the solution equation yields \[\begin{align*} u^{n+1} &= \Bigl(\frac{n-1}{2n-1} \nu^s +\frac{n}{2n-1} \nu^{s-k+m}\Bigr) u^n \\ & \quad + \frac{n-1}{2n-1}\sum_{j=1}^{k-2} \nu^{s-k+j+1} r_{k-j} + \frac{n}{2n-1}\sum_{j=1}^{m-1} \nu^{s-k+j} r_{m-j+1} + \sum_{j=1}^{s-k+1} \nu^j r_{s-j+1}. \end{align*}\] It is convenient to change the indices in the summations and after some algebra we rewrite the above equation as \[\begin{align} u^{n+1} &= \Bigl(\frac{n-1}{2n-1} \nu^{n^2} +\frac{n}{2n-1} \nu^{(n-1)^2}\Bigr) u^n + \\ &\quad \sum_{j=2}^{m}\Bigl(\frac{n-1}{2n-1} \nu^{n^2-j+1} +\frac{n}{2n-1} \nu^{(n-1)^2-j+1}\Bigr) r_j + \frac{n-1}{2n-1}\sum_{j=m+1}^{k-1} \nu^{n^2-j+1} r_j + \sum_{j=k}^{n^2} \nu^{n^2-j+1} r_j. \end{align}\]

Then the stability function is \[ P(\nu) = \frac{n-1}{2n-1} \nu^{n^2} +\frac{n}{2n-1} \nu^{(n-1)^2} \] and the internal stability functions are \[\begin{align*} Q_{sj}(\nu) = \begin{cases} \dfrac{n-1}{2n-1} \nu^{n^2-j+1} +\dfrac{n}{2n-1} \nu^{(n-1)^2-j+1}, & 2 \leq j \leq m \\\\ \dfrac{n-1}{2n-1} \nu^{n^2-j+1}, & m+1 \leq j \leq k-1 \\\\ \nu^{n^2-j+1}, & k \leq j \leq n^2. \end{cases} \end{align*}\]

We would like to know the maximum amplification factor; i.e. the maximum of \(|Q_{sj}(\nu)|\) for \(\nu \in \mathbb{C}\) such that \(|P(\nu)|\le 1\).

\(\left|\frac{n-1}{2n-1} \nu^{n^2} +\frac{n}{2n-1} \nu^{(n-1)^2}\right| =1\) stability region boundaries for \(2\le n\le 8\)

For a treatment on this case (just the steps of a possible proof), see this link as well.

Since the latter condition always gives a bounded region, in order to get a bound for \(|\nu|\) we seek
a circle that contains the stability region for any value of \(n\). Because \(\nu\) can be in absolute value bigger than one and still to satisfy the condition \(|P(\nu)|\le 1\), we want a circle that is slightly bigger than \[|z + n^2 - n| < n^2 - n.\] It can be shown that \[|z + n^2 - n| < n^2 - \frac{n^2}{n+4} = \frac{n^2(n+3)}{n+4}\] always includes the stability region (tested for values of \(n\) up to \(100\)). Then \[|\nu| < \frac{n(n+3)}{(n-1)(n+4)}\] and since all internal stability functions are less than \(|\nu|^{n^2-1}\), then \[|Q_{sj}| < \Bigl(\frac{n(n+3)}{(n-1)(n+4)}\Bigr)^{n^2-1}.\] This seems to be a rather poor estimate for values of \(n \leq 5\) for which a bound of \[|Q_{sj}| < \Bigl(\frac{n(n+3)}{(n-1)(n+4)}\Bigr)^{\frac{n^2-n}{2}}\] gives a better estimate.

% number of stages: n^2
n = 8;

phi = zeros(1000,1000);
phi2 = zeros(1000,1000);
v = zeros(1000,1000);

x = linspace(-500,500,1000);
y = linspace(-500,500,1000);
[x y] = meshgrid(x,y);
z = x + y*1i;

for i =1 :1000
    for j=1:1000
        v(i,j) = 1 + z(i,j)/(n^2-n);
    end
end

for i =1 :1000
    for j=1:1000
        phi(i,j) = (n-1)/(2*n-1)*v(i,j)^(n^2) + n/(2*n-1)*v(i,j)^((n-1)^2);
        phi2(i,j) = z(i,j)+n^2-n;
    end
end

contour(x,y,abs(phi), [1,1], 'r')
hold on
contour(x,y,abs(phi2), [n^2-n^2/(n+4),n^2-n^2/(n+4)], 'b')
xlabel('Real z'), ylabel('Imag z');
legend('stability region', 'Location', 'best');
hold all;
axis tight
grid on;
axis square;