david/A heuristic method for splitting RK methods into upwind and downwind parts
Preliminaries
We start with a Runge–Kutta method
\[\begin{align*} \label{RK} y & = e u_n + h\beta_0 f.\end{align*}\]
When studying absolute monotonicity, it is convenient to write method in canonical Shu-Osher form
\[\begin{align*} \label{cso} y & = d_r u_n + \alpha_r \left( y + \frac{h}{r} f\right)\end{align*}\]
where
\[\begin{align*} \alpha_r & = r(I + r\beta_0)^{-1} \beta_0 \\ d_r & = (I + r\beta_0)^{-1} e.\end{align*}\]
We consider additive Runge–Kutta methods
\[\begin{align*} \label{additiveRK} y & = e u_n + h{\beta^\textup{up}}f + h{\beta^\textup{down}}{\tilde{f}}\end{align*}\]
characterized by the coefficients \({\beta^\textup{up}}, {\beta^\textup{down}}\). The underlying Runge–Kutta method corresponding to is obtained by setting \({\tilde{f}}=-f\). Conversely, we say that is an upwind/downwind splitting, or just a splitting, of if it reduces to when \({\tilde{f}}=-f\), i.e. if
\[\begin{align*} \label{betadiff} {\beta^\textup{up}}-{\beta^\textup{down}}=\beta_0.\end{align*}\]
We further define the SSP coefficient of as the largest value \(R\) such that the upwind/downwind method produces a solution satisfying \[\|u_{n+1}\|\le\|u_n\|\] whenever \(h\le R h_0\) and \(f,{\tilde{f}}\) satisfy
\[\begin{align*} \|u + hf(u)\| & \le \|u\| \mbox{ for all } u \in V, 0\le h\le h_0 \\ \|u + h{\tilde{f}}(u)\| & \le \|u\| \mbox{ for all } u \in V, 0\le h\le h_0.\end{align*}\]
It turns out that \(R\) depends only on \({\beta^\textup{up}}, {\beta^\textup{down}}\). We would like to find, for a given \(\beta_0\), the largest value of \(R\) over all \({\beta^\textup{up}}, {\beta^\textup{down}}\) satisfying . This question was posed by Higueras .
It is helpful to consider a form of similar to the form . The following appears as Lemma 3.1 in ; it follows from Theorems 3.5 and 3.6 of .
The SSP coefficient of the upwind/downwind Runge-Kutta method is the largest \(r\ge0\) such that the method can be written in the form
\[\begin{align*} \label{dwrk_cso} y & = \gamma_r u_n + \alpha^\textup{up}_r \left(y+ \frac{h}{r}f\right) + \alpha^\textup{down}_r \left(y+ \frac{h}{r}{\tilde{f}}\right)\end{align*}\]
where
\[\begin{align*} \label{poscoef} \gamma_r,\alpha^\textup{up}_r,\alpha^\textup{down}_r\ge 0.\end{align*}\]
The following proposition characterizes which methods of the form are splittings of a given method .
Method is a splitting of method if and only if
\[\begin{align*} \label{alpharel1} (I-2\alpha^\textup{down}_r)\alpha_r & = (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \\ (I-2\alpha^\textup{down}_r)d_r & = \gamma_r.\end{align*}\]
Take \({\tilde{f}}= -f\) in to obtain:
\[\begin{align*} y & = \gamma_r u_n + (\alpha^\textup{up}_r + \alpha^\textup{down}_r) y + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \frac{h}{r}f.\end{align*}\]
Subtract \(2\alpha^\textup{down}_r y\) from both sides to get
\[\begin{align*} (I-2\alpha^\textup{down}_r) y & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \left(y + \frac{h}{r}f\right).\end{align*}\]
Substituting in the above gives
\[\begin{align*} (I-2\alpha^\textup{down}_r) d_r u_n + (I-2\alpha^\textup{down}_r)\alpha_r \left(y + \frac{h}{r}f\right) & = \gamma_r u_n + (\alpha^\textup{up}_r - \alpha^\textup{down}_r) \left(y + \frac{h}{r}f\right),\end{align*}\]
Equating coefficients yields .
For most (including all explicit) methods, \(I-2\alpha^\textup{down}_r\) is invertible and can be written
\[\begin{align*} \label{alpharel} \alpha_r & = (I-2\alpha^\textup{down}_r)^{-1}(\alpha^\textup{up}_r - \alpha^\textup{down}_r) \\ d_r & = (I-2\alpha^\textup{down}_r)^{-1}\gamma_r.\end{align*}\]
From now on, we’ll just assume invertibility.
There is a small defect in the ’only if’ part of the above proposition. If \(\alpha_r\) has a zero row, then one of the stages is equal to \(u_n\) and it is possible to exchange values between \(\gamma_r\) and the corresponding column of \(\alpha^\textup{up}_r, \alpha^\textup{down}_r\) without changing the underlying method.
The question, then, is how to choose among all splittings satisfying for a given \(r\), one that yields positive coefficients. The following result shows that we only need to focus on a certain way of splitting \(\alpha_r\).
[signlemma] Given an explicit Runge–Kutta method , Let \(\alpha^\textup{up}_r\ge0, \alpha^\textup{down}_r\ge0\) denote coefficients of a splitting of . Then
[alphar]
\[\begin{align*} \alpha^\textup{up}_r & = (I-2\alpha^-_r)^{-1} \alpha^+_r \\ \alpha^\textup{down}_r & =-(I-2\alpha^-_r)^{-1} \alpha^-_r,\end{align*}\]
where \(\alpha^+_r\ge 0\), \(\alpha^-_r\le 0\), and \(\alpha_r=\alpha^+_r + \alpha^-_r\).
Define
[a1a2]
\[\begin{align*} \alpha^+_r & = (I-2\alpha^\textup{down})^{-1} \alpha^\textup{up}_r \\ \alpha^-_r & =-(I-2\alpha^\textup{down})^{-1} \alpha^\textup{down}_r.\end{align*}\]
Then \(\alpha_r=\alpha^+_r+\alpha^-_r\). Furthermore, since \(I-2\alpha^\textup{down}_r\) is an \(M\)-matrix, we have \(\alpha^+_r\ge0\) and \(\alpha^-_r\le0\). Solving for \(\alpha^\textup{up}_r,\alpha^\textup{down}_r\) gives .
Optimal splittings
The foregoing suggests an approach to obtain optimal splittings. Given a method and some \(r>0\), we first write the method in canonical Shu-Osher form and then take
[alpharchoice]
\[\begin{align*} \alpha^\textup{up}_r & = (I-2(\alpha_r)^-)^{-1} (\alpha_r)^+ \\ \alpha^\textup{down}_r & =-(I-2(\alpha_r)^-)^{-1} (\alpha_r)^-, \gamma_r & = (I-2\alpha^\textup{down}_r) d_r.\end{align*}\]
where
\[\begin{align*} ((x)^+)_{ij} = \begin{cases} x_{ij} & \mbox{ if } x_{ij}\ge0 \\ 0 & \mbox{ if } x_{ij}<0. \end{cases} \\ ((x)^-)_{ij} = \begin{cases} 0 & \mbox{ if } x_{ij}\ge0 \\ x_{ij} & \mbox{ if } x_{ij}<0. \end{cases}\end{align*}\]
These matrices satisfy the conditions in Lemma [signlemma]. If they are non-negative, then the split method has SSP coefficient at least \(r\). A natural approach is to find the largest \(r\) such that non-negativity of the coefficient matrices holds.
For methods (such as all explicit methods) with a zero row, a slight modification is useful. Suppose that \(y_1=u_n\). Then the underlying method is unchanged under the transformation
\[\begin{align*} (\gamma_r)_i & \to (\gamma_r)_i - \epsilon \\ (\alpha^\textup{up}_r)_{i1} & \to (\alpha^\textup{up}_r)_{i1} + \epsilon/2 \\ (\alpha^\textup{down}_r)_{i1} & \to (\alpha^\textup{down}_r)_{i1} + \epsilon/2.\end{align*}\]
For any \(i\) such that \((\gamma_r)_i >0\), it is clearly advantageous to apply this tranformation with \(\epsilon=(\gamma_r)_i\).
This approach seems to lead to optimal splittings in all tested cases. We don’t know that it’s optimal because one could use
\[\begin{align*} (\alpha_r)^+ + \delta \\ (\alpha_r)^- - \delta,\end{align*}\]
in place of \((\alpha_r)^+, (\alpha_r)^-\), where \(\delta\) is any non-negative matrix.
#Old Stuff below here
Given a method \(K\), and \(r\ge0\), define \[d = (I+rK)^{-1}e\] \[P = r(I+rK)^{-1}K.\] Furthermore, let \(P^+, P^-\) denote the positive and negative parts of \(P\), such that \(P^+\ge0\), \(P^-\le 0\), and \(p=P^+ + P^-\). Then set \[K^+ = (I-2P^-)^{-1}P^+\] \[K^- = (I-2P^-)^{-1}P^-\] \[\gamma = (I-2P^-)^{-1}d\]
Then the method can be written as \[Y = \gamma u^n + K^+ (Y+\frac{\Delta t}{r} F) + K^- (Y-\frac{\Delta t}{r} F)\] Replacing the second \(F\) with \(-\tilde{F}\) gives \[Y = \gamma u^n + K^+ (Y+\frac{\Delta t}{r} F) + K^- (Y+\frac{\Delta t}{r} \tilde{F})\] The method is SSP with coefficient \(r\) if \(\gamma, K^+, K^-\) are all non-negative.
(can state an if theorem; want to get an only if result)
For explicit methods, positive entries in \(d\) can advantageously be reallocated to the first column of \(K^\pm\), since both coefficients multiply \(u^n\). If corresponding entries of \(K^+\) and \(K^-\) are incremented by the same amount, then the corresponding \(f^+,f^-\) terms cancel out. So the optimal thing to do is allocate half of \(d\) to each of \(K^+,K^-\).
DERIVATION HERE