docu.tex

\documentclass{article}
\usepackage{amsmath,amssymb,amsthm}
\usepackage[utf8]{inputenc}
\usepackage[ruled,algo2e,linesnumbered,algonl]{algorithm2e}
\usepackage{color}

\newtheorem{remark}{Remark}
\newtheorem{lemma}{Lemma}
\newcommand{\ve}[1]{\mathbf{#1}}
\newcommand{\todo}[1]{\textcolor{red}{[TODO -- #1]}}

%\bibliographystyle{plain}

%\bibliography{refs}
\begin{document}
\section*{Multi-rate spectral deferred corrections}
Consider a single time step $[T_{n}, T_{n+1}]$.
Now denote as $t_m$, $m=1, \ldots, M$ with $T_{n} \leq t_1 < \ldots < t_{M} \leq T_{n+1}$ a set of quadrature nodes in this time step.
Denote as $l_m$, $m=1, \ldots, M$ the Lagrange polynomials satisfying 
\begin{equation}
	l_m(t_j) = \delta_{m,j} \ \text{for} \ m,j=1, \ldots, M.
\end{equation}
In each sub step $[t_{m-1}, t_{m}]$ with length $\Delta t_{m}=t_{m}-t_{m-1}$, place a set of embedded nodes $t_{m,p}$, $m=1,\ldots,M$, $p=1,\ldots,P$ with $t_{m-1} \leq t_{m,1} < \ldots < t_{m,P} \leq t_{m}$.
Denote as $l_{m,p}$ the Lagrange polynomials for an embedded set of quadrature nodes satisfying
\begin{equation}
	l_{m,p}(t_{m,q}) = \delta_{p,q} \ \text{for} \ p,q=1,\ldots,P.
\end{equation}
\begin{remark}
We refer to nodes and weights associated with $t_m$ as \textbf{standard} and use indices $m$ and $j$.
Nodes and weights associated with an embedded quadrature rule are called \textbf{embedded} and indexed with $p$ and $q$.
\end{remark}
Note that there are $M$ many standard nodes, polynomials and weights and $M \times P$ many embedded nodes, polynomials and weights.
\begin{remark}
We assume that for both standard and embedded quadrature rules, the right endpoint is a quadrature nodes, so that $t_{M} = T_{n+1}$ and $t_{m,P} = t_{m}$.
\end{remark}
\begin{remark}
However, to allow for Radau-type rules, we do \textbf{not} require that the left endpoint matches the first quadrature nodes and typically we will have $T_n < t_1$ as well as $t_{m-1} < t_{m,1}$.
To simplify notation, we use the following conventions
\begin{equation}
	t_0 := T_n \ \text{and} \ t_{m,0} := t_{m-1}.
\end{equation}
\end{remark}

\subsection*{Weights.}
We need three different sets of quadrature weights.
First the standard weights
\begin{equation}
	s_{m,j} := \int_{t_{m-1}}^{t_{m}} l_{j}(s)~ds, \ \text{for} \ m,j=1,\ldots,M.
\end{equation}
For the embedded weights we need
\begin{equation}
	s_{m,p,q} := \int_{t_{m,p-1}}^{t_{m,p}} l_{m,q}(s)~ds, \ \text{for} \ m=1,\ldots,M, \ p,q=1,\ldots,P.
\end{equation}
To approximate integration of a function given at the embedded nodes over a standard sub-step, the embedded weights need to be summed up
\begin{equation}
	\hat{s}_{m,q} := \int_{t_{m-1}}^{t_{m}} l_{m,q}(s)~ds = \sum_{p=1}^{P} \int_{t_{m,p-1}}^{t_{m,p}} l_{m,q}(s)~ds = \sum_{p=1}^{P} s_{m,p,q}
\end{equation}
for $m=1,\ldots,M$ and $q=1,\ldots,P$.
We also need the following \textbf{mixed} weights for integrating the Lagrange polynomials for the standard weights over sub steps associated with embedded quadrature rules
\begin{equation}
	\tilde{s}_{m,p,j} := \int_{t_{m,p-1}}^{t_{m,p}} l_{j}(s)~ds, \ \text{for} \ m,j=1,\ldots,M, \ p=1,\ldots,P.
\end{equation}
\begin{lemma}\label{lemma:weights_match}
The mixed weights are consistent with the standard weights in the sense that
\begin{equation}
	s_{m,j} = \int_{t_{m-1}}^{t_{m}} l_j(s)~ds = \sum_{p=1}^{P} \int_{t_{m,p-1}}^{t_{m,p}} l_j(s)~ds = \sum_{p=1}^{P} \tilde{s}_{m,p,j}.
\end{equation}
Verified by \texttt{test\_weights\_match}.
\end{lemma}

\begin{table}[h]
\centering
\begin{tabular}{|cc|cc|} \hline
\multicolumn{2}{|c|}{Function values} & \multicolumn{2}{c|}{Integral} \\ \hline
&          & Standard        & Embedded            \\ \hline
& Standard & $s_{m,j}$       & $\tilde{s}_{m,p,j}$ \\
& Embedded & $\hat{s}_{m,q}$ & $s_{m,p,q}$ \\ \hline
\end{tabular}
\caption{Quadrature weights for integrals over standard or embedded sub steps depending on whether function values are given at standard or embedded nodes.}
\end{table}

Furthermore, in case we want to use the ``zero-to-node'' update formulation, we also define the following weights
\begin{equation}
	q_{m,j} := \int_{T_{n}}^{t_m} l_j(s)~ds.
\end{equation}
Simple addition shows that
\begin{equation}
	q_{m,j} = \sum_{n=1}^{m} s_{n,j}.
\end{equation}
We also define the following embedded weights
\begin{equation}
	q_{m,p,q} := \int_{t_{m-1}}^{t_{m,p}} l_{m,q}(s)~ds = \int_{t_{m,0}}^{t_{m,p}} l_{m,q}(s)~ds
\end{equation}
and note that
\begin{equation}
	q_{m,p,q} = \sum_{r=1}^{p} s_{m,r,q}.
\end{equation}
%
% Quadrature nodes
%
\subsection*{Quadrature rules.}
Denote as $\ve{u}^{s}$ (``standard'') a vector with $M$ approximate solutions at the standard nodes and as $\ve{u}^{e}$ (``embedded'') a vector composed of $M$ vectors $\ve{u}^{e}_{m}$, $m=1,\ldots,M$ with each $\ve{u}^{e}_{m}$ containing $P$ approximate solutions at the embedded nodes $t_{m,p}$, $p=1, \ldots, P$.
We will need the following three integration operators:
\begin{equation}
	\int_{t_{m-1}}^{t_{m}} u(s)~ds \approx I_{m-1}^{m}(\ve{u}^s, \ve{u}^{e}_{m}) := \sum_{j=1}^{M} s_{m,j} \ve{u}^{s}_{j} + \sum_{p=1}^{P} \hat{s}_{m,p} \ve{u}^{e}_{m,p}.
\end{equation}
Here, $\ve{u}^{s}_{j}$ denotes the $j$\textsuperscript{th} entry in $\ve{u}^{s}$ while $\ve{u}^{e}_{m,p}$ is the $p$\textsuperscript{th} entry in $\ve{u}^{e}_{m}$.
\begin{remark}
Test \texttt{...} verifies this identity for linear functions.
\end{remark}
Furthermore, we can approximate the integral over an embedded sub step by
\begin{equation}
	\int_{t_{m,p-1}}^{t_{m,p}} u(s)~ds \approx I_{m,p-1}^{p}(\ve{u}^s, \ve{u}^{e}_{m}) := \sum_{j=1}^{M} \tilde{s}_{m,p,j} \ve{u}^{s}_{j} + \sum_{q=1}^{P} s_{m,p,q} \ve{u}^{e}_{m,q}.
\end{equation}
\begin{lemma}\label{lemma:quadrature_match}
The two integration operators are consistent in the sense that
\begin{equation}
	I_{m-1}^{m}(\ve{u}^s, \ve{u}^{e}_{m}) = \sum_{p=1}^{P} I_{m,p-1}^{p}(\ve{u}^s, \ve{u}^{e}_{m}).
\end{equation}
\begin{proof}
By definition of the operators and $\hat{s}_{m,q}$ and using Lemma~\ref{lemma:weights_match} it holds that
\begin{align*}
	\sum_{p=1}^{P} I_{m,p-1}^{p}(\ve{u}^s, \ve{u}^{e}_{m}) &= \sum_{p=1}^{P} \left( \sum_{j=1}^{M} \tilde{s}_{m,p,j} \ve{u}^{s}_j + \sum_{q=1}^{P} s_{m,p,q} \ve{u}^{e}_{m,q} \right) \\
		&= \sum_{j=1}^{M} \left( \sum_{p=1}^{P} \tilde{s}_{m,p,j} \right) \ve{u}^{s}_{j} + \sum_{q=1}^{P} \left( \sum_{p=1}^{P} s_{m,p,q} \right) \ve{u}^{e}_{m,q} \\
		&= \sum_{j=1}^{M} s_{m,j} \ve{u}^{s}_{j} + \sum_{q=1}^{P} \hat{s}_{m,q} \ve{u}^{e}_{m,q} = I_{m-1}^{m}(\ve{u}^s, \ve{u}^{e}_{m})
\end{align*}
\end{proof}
\end{lemma}
%
%
%
\paragraph{Matrix formulation.}
To represent the above quadrature rules in matrix form, define the following four matrices
\begin{align*}
	\ve{S} &:= \left( s_{m,j} \right) \in \mathbb{R}^{M,M} \ \text{for} \ m,j=1,\ldots,M \\ 
	\hat{\ve{S}} &:= \left( \hat{s}_{m,p} \right) \in \mathbb{R}^{M,P} \ \text{for} \ m=1,\ldots,M, \ p=1,\ldots,P \\
	\tilde{\ve{S}}_m &:= \left( \tilde{s}_{m,p,j} \right) \in \mathbb{R}^{P,M} \ \text{for} \ \ p=1,\ldots,P, \ j=1,\ldots,M \\
	\ve{S}_m &:= \left( s_{m,p,q} \right) \in \mathbb{R}^{P,P} \ \text{for} \ p,q=1,\ldots,P.
\end{align*}
We can write the integrals over standard sub steps as
\begin{equation}
	\ve{I}^s := 
	\begin{bmatrix}
		I_{0}^1 \\ I_1^2 \\ \vdots \\ I_{M-1}^{M}
	\end{bmatrix}
	=
	\ve{S} \ve{u}^s +
	\begin{bmatrix}
		\hat{\ve{S}}[1,:] & & \\ & \hat{\ve{S}}[2,:] \\ & & \ddots \\ & & & \hat{\ve{S}}[M,:] 
	\end{bmatrix}
	\begin{bmatrix}
		\ve{u}^e_1 \\ \ve{u}^e_2 \\ \vdots \\ \ve{u}^{e}_{M}
	\end{bmatrix}.
\end{equation}
Integrals over embedded steps can be written as
\begin{equation}
	\ve{I}^e_m :=
	\begin{bmatrix}
		I_{m,0}^{1} \\ I_{m,1}^{2} \\ \vdots \\ I_{m,P-1}^{P}
	\end{bmatrix}
	=
	\tilde{\ve{S}}[m,:,:] \ve{u}^s + \ve{S}[m,:,:] \ve{u}^e_m, \ \text{for} \ m=1,\ldots,M
\end{equation}
so that
\begin{equation}
	\ve{I}^e = 
	\begin{bmatrix}
		\ve{I}^e_1 \\ \ve{I}^e_2 \\ \vdots \\ \ve{I}^e_M
	\end{bmatrix} = 
	\begin{bmatrix}
		\tilde{\ve{S}}[1,:,:] \\ \tilde{\ve{S}}[2,:,:]  \\ \vdots \\ \tilde{\ve{S}}[M,:,:] 
	\end{bmatrix}
	\ve{u}^s
	+
	\begin{bmatrix}
		\ve{S}[1,:,:] \\ \ve{S}[2,:,:] \\ \vdots \\ \ve{S}[M,:,:] 
	\end{bmatrix}
	\begin{bmatrix}
		\ve{u}^e_1 \\ \ve{u}^e_2 \\ \vdots \\ \ve{u}^e_M.
	\end{bmatrix}
\end{equation}
%
%
%
\subsection*{Multi-rate SDC sweeps}
Consider an initial value problem
\begin{equation}
	\dot{u}(t) = f^I(u(t))+f^E(u(t)) + g(u(t))
\end{equation}
where we want to integrate $f=f^I+f^E$ with a large time step and $g$ explicitly with a small time step. Analogous to the IMEX-sdc method, we split the slow parts in an implicit part $f^I$ and an explicit part $f^E$.
As before, denote as $\ve{u}^s$ a vector with approximate solutions at the standard nodes and as $\ve{u}^{e}$ the vector with approximate solutions at the embedded nodes.
The predictor computes the initial set of values
\begin{itemize}
	\item Compute standard step
	\begin{equation*}
		u^*_{m} = u^0_{m-1} + \Delta t_m \left(f^I(u^*_{m})+f^E(u^0_{m-1})\right)
	\end{equation*}
	\item Compute $P$ embedded steps
	\begin{equation*}
		u^0_{m-1,p} = u^0_{m-1,p-1} + \Delta t_{m,p} \left(f^I(u^*_{m})+f^E(u_{m-1})\right) + \Delta t_{m,p} g(u^0_{m-1,p-1})
	\end{equation*}
	\item Set $u^0_{m} = u^0_{m-1,P}$ for use as initial value in next standard step.
\end{itemize}

Given values $f^I(u^k_m)$, $f^E(u^k_m)$ and $g(u^k_{m,p})$, the multi-rate SDC sweep then consists of
\begin{itemize}
\item Compute standard step
	\begin{equation}
	\begin{aligned}
		u^{*}_{m} = & u^{k+1}_{m-1} + \Delta t_m \left( f^I(u^{*}_{m}) - f^I(u^k_{m}) + f^E(u^{k+1}_{m-1})-f^E(u^{k}_{m-1}) \right) \cr + & I_{m-1}^{m} \left( f(\ve{u}^{s,k}) , g(\ve{u}^{e,k}_m) \right)
    \end{aligned}
	\end{equation}
\item Compute $P$ embedded steps
	\begin{equation}
	\begin{aligned}
		u_{m-1,p}^{k+1} = u_{m-1,p-1}^{k+1} &+ \Delta t_{m,p} \left( f^I(u^{*}_{m}) - f^I(u^k_{m}) + f^E(u^{k+1}_{m-1})-f^E(u^{k}_{m-1}) \right) \\
			& + \Delta t_{m,p} \left( g(u^{k+1}_{m-1,p-1}) - g(u^k_{m-1,p-1}) \right) \\
			& + I_{m-1,p-1}^{p}  \left( f(\ve{u}^{s,k}) , g(\ve{u}^{e,k}_{m-1}) \right) 
	\end{aligned}
	\end{equation}
	for $p=1, \ldots, P$.
\item Set $u_{m}^{k+1} = u_{m-1,P}^{k+1}$ for use as initial value in next standard step.
\end{itemize}
%
%
%
\begin{algorithm2e}
	\caption{Multi-rate SDC prediction step.}
        \SetKwInOut{Input}{input}
         \SetKwInOut{Output}{output}
         \Input{$u_0$}
         \Output{$f^I(u^0_m)$, $f^E(u^0_m)$ and $g(u^0_{m,p})$ for $m=1,\ldots,M$ and $P=1,\ldots,P$.}
         $u^0_0 \leftarrow u_0$\\
         \For{$m=1, M$}{
         	Solve $u^* = u^0_{m-1} + \Delta t_m f(u^*)$ \\
		$f^* \leftarrow f(u^*)$ \\
		$u^0_{m-1,0} \leftarrow u^0_{m-1}$\\
		\For{$p=1,P$}{
			$u^0_{m-1,p} = u^0_{m-1,p-1} + \Delta t_{m,p} f^*+ \Delta t_{m,p} g(u^0_{m-1,p-1})$ \\
			Evaluate $g(u^0_{m-1,p})$
		}
		$u^0_m \leftarrow u^0_{m-1,P}$ \\
		Evaluate $f^I(u^0_m)$ and $f^E(u^0_m)$
         }
\end{algorithm2e}
%
%
%
\begin{algorithm2e}
	\caption{Multi-rate SDC sweep.\todo{initial value for first embedded step is not correct... explicit terms cancel out}}
        \SetKwInOut{Input}{input}
         \SetKwInOut{Output}{output}
         \Input{$u_0$ and $f^I(u^k_m)$, $f^E(u^k_m)$, $g(u^k_{m,p})$ for $m=1,\ldots,M$ and $p=1,\ldots,P$.}
         \Output{updated values $f(u^{k+1}_m)$, $g(u^{k+1}_{m,p})$}
         Compute $I_{m-1}^{m}$, $m=1, \ldots, M$\\
         Compute $I_{m-1,p-1}^{p}$, $m=1,\ldots,M$; $p=1,\ldots,P$ \\
         $u^{k+1}_0 \leftarrow u_0$ \\
         \For{$m=1,M$}{
         	Solve $u^* = u^{k+1}_{m-1} + \Delta t_m \left( f^I(u^{*}_{m}) - f^I(u^k_{m})+f^E(u^{k+1}_{m-1})-f^E(u^k_{m-1}) \right) + I_{m-1}^{m}$ \\
		$f^* \leftarrow f^I(u^*)-f^I(u^k_m)+f^E(u^{k+1}_{m-1})-f^E(u^{k}_{m-1})$\\
		$u^{k+1}_{m-1,0} \leftarrow u^{k+1}_{m-1}$\\
		\For{$p=1,P$}{
			$ 
			 \begin{aligned}
			u_{m-1,p}^{k+1} =& u_{m-1,p-1}^{k+1} + \Delta t_{m,p} f^* \cr +&  \Delta t_{m,p} \left( g(u^{k+1}_{m-1,p-1}) - g(u^k_{m-1,p-1}) \right)  + I_{m,p-1}^{p}
			\end{aligned}$\\
			Evaluate $g(u^{k+1}_{m-1,p})$
		}
		$u^{k+1}_m \leftarrow u^{k+1}_{m-1,P}$\\
		Evaluate $f^I(u^{k+1}_m)$, $f^E(u^{k+1}_m)$
         }
\end{algorithm2e}
%
%
%
\begin{remark}
If the iteration converges, that is and $u^{k+1}_{m} - u^{k}_{m} \to 0$ and $u^{k+1}_{m,p} - u^{k}_{m,p} \to 0$, the sweeps reduce to the underlying collocation equations
\begin{equation}
	u_{m} = u_{m-1} + I_{m-1}^{m} \left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right)
\end{equation}
and
\begin{equation}
	u_{m,p} = u_{m,p-1} + I_{m,p-1}^{p} \left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right)
\end{equation}
Based on this, for given values of $\ve{u}^s$, $\ve{u}^3$, we can define the standard 
\begin{equation}
 r^s := \max_{m=1,\ldots,M} \left\| \ve{u}^s_{m} - \ve{u}^s_{m-1} - I_{m-1}^{m}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right) \right\|
\end{equation}
and embedded residual
\begin{equation}
	r^e := \max_{m=1,\ldots,M} \max_{p=1,\ldots,P} \left\| \ve{u}^e_{m,p} - \ve{u}^e_{m,p-1} - I_{m,p-1}^{p}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right) \right\|.
\end{equation}
\end{remark}
%
%
%
\begin{remark}
In case a set of function values $f(\ve{u}^{s})$, $g(\ve{u}^e)$ and an initial value $u_0$ is given, the solution can be reconstructed iteratively from
\begin{equation}
	u_{m,p} = u_{m,p-1} + I_{0,p-1}^{p}
\end{equation}
for $p=1, \ldots, P$ with $u_{0,0} = u_0$ in the first standard step and $u_{m,0} = u_{m-1}$ for $m > 1$.
After the $u_{m,p}$ are computed, the value for the next standard step is obtained from
\begin{equation}
	u_{m} := u_{m,P-1}.
\end{equation}
\todo{hm... how does this help?}
\end{remark}
%
%
%
\begin{remark}
 In case of vanishing fast part, i.e. $g=0$, we obtain the classical IMEX-SDC.
\end{remark}
%
%
%
\begin{lemma}
The collocation solution is invariant under SDC sweeps. Let $\ve{u}^s$ and $\ve{u}^e$ satisfy the collocation condition.
Then,
\begin{align*}
	u^{k+1}_{m} &= u_{m-1}  + \Delta t_m \left( f(u^{k+1}_m) - f(u_m) \right) + I_{m-1}^{m}\left( f(\ve{u}^s, g(\ve{u}^e_m) \right) \\
		&= u_{m} + \Delta t_{m} \left( f(u^{k+1}_m) - f(u_m) \right) \\
	\Rightarrow u^{k+1}_m - \Delta t_m f(u^{k+1}_{m}) &= u_{m} - \Delta t_{m} f(u_m)
\end{align*}
so that $u^{k+1}_m = u_m$.
The embedded sweep will then give
\begin{align*}
	u^{k+1}_{m,p} &= u_{m,p-1} + \Delta t_{m,p} \left( g(u^{k+1}_{m-1,p-1}) - g(u_{m-1,p-1} \right) + I_{m,p-1}^{p} \left( f(\ve{u}^s, g(\ve{u}^e_m) \right) \\
				&= u_{m,p} - \Delta t_{m,p} \left( g(u^{k+1}_{m-1,p-1}) - g(u_{m-1,p-1}) \right)
\end{align*}
so that $u^{k+1}_{m,p} = u_{m,p}$.
Finally,
\begin{equation}
	u_{m} = u_{m-1,P} = u_{m-1,0} + \sum_{p=1}^{P} I_{m,p-1}^{p}\left( f(\ve{u}^s, g(\ve{u}^e_m) \right) = u_{m-1} + I_{m-1}^{m}\left( f(\ve{u}^s, g(\ve{u}^e_m) \right) 
\end{equation}
so that overwriting $u_{m}$ with $u_{m-1,P}$ does not modify the value.
\end{lemma}
%
%
%
\begin{lemma}
The collocation solutions are consistent in the sense that
\begin{align*}
	u_{m,P} &= u_{m-1,P-1} + I_{m,P-1}^{P}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right) \\
			&= \ldots = u_{m-1,0} + \sum_{p=1}^{P} I_{m,p-1}^{p}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right) \\
			&= u_{m-1} + I_{m-1}^{m}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right) = u_{m}
\end{align*}
using Lemma~\ref{lemma:quadrature_match}.
\end{lemma}
%
%
\begin{remark}
For approximations of the collocation solution, this consistency criterion is not satisfied and typically
\begin{equation}
	u_{m}^{k} \neq u_{m,P}^{k}.
\end{equation}
\end{remark}

\paragraph{A comment on updates.}
Instead of simply overwriting $u^{k+1}_{m}$ with $u^{k+1}_{m-1,P}$ in the multi-rate SDC sweep, it is possible to compute the proper collocation update
\begin{equation}
	u^{k+1}_{m} = u^{k+1}_{m-1} + \sum_{p=1}^{P} q_{m, p} f(u^{k+1}_{m,p})
\end{equation}
with weights
\begin{equation}
	q_{m,p} := \int_{t_{m-1}}^{t_{m}} l_{m,p}(s)~ds.
\end{equation}
\todo{effect??}

\subsection*{Zero-to-node formulation}
Recursively applying the node-to-node collocation equations provides the zero-to-node formulation for both standard and embedded quadrature rules.
For the embedded step, this gives us
\begin{subequations}
\begin{align}
	u^{k+1}_{m-1,p} &= u^{k+1}_{m-1,0} + \left( f^{I}(u^{*}_{m}) - f^{I}(u^k_m)  \right)  \sum_{r=1}^{p} \Delta t_{m,p} \\
				& \quad +  \left( f^{E}(u^{k+1}_{m-1}) - f^{E}(u^k_{m-1}) \right) \sum_{r=1}^{p} \Delta t_{m,p} \\
				& \quad + \sum_{r=1}^{p} \Delta t_{m,p} \left( g(u^{k+1}_{m-1,r-1}) - g(u^k_{m-1,r-1}) \right) \\
				& \quad + \sum_{r=1}^{p} I_{m-1,r-1}^{r}\left( f(\ve{u}^{s}) , g(\ve{u}^{e}_m) \right)
\end{align}
\end{subequations}
Note that since the $f^I$ and $f^E$ terms do not depend on $p$, we can simply move them in front of the sum.
To use St. Martin's trick, replace the $\Delta t_{m,p}$ with the entries from the LU factorisation.
\todo{Because we set $u^{k+1}_m$ as the final value from the embedded sweep, I don't see how we can rewrite this as zero-to-node? Might have to settle for using St. Martin's trick for embedded sweep only.}

\todo{Have we tried using Picard (i.e $\theta=0$) in standard sweep but $\theta = 1.0$ in embedded sweeps?}

\subsection*{Sweep convergence}
Analogous to Martin Weiser, we apply the MRSDC method on Dahlquists equation. For our method, we have to use the modified equation
\begin{align}
 \dot{u}=&\lambda_1 u + \lambda_2 u
\end{align}
with $f^I(u)=\lambda_1 u$, $f^E(u)=0$ and $g(u)=\lambda_2 u$. (\todo{what about $f^E(u)=\lambda_3u$?}) 
For given $u^k_{m}$ and $u^{k}_{m,p}$, we obtain the expressions
\begin{align*}
 u^*_m=&\frac{1}{1-\lambda_1\Delta t_m }\left(u^{k+1}_{m-1}-\lambda_1\Delta t_m u^{k}_m+ I_{m-1}^m\right) \cr
 f^*_m =&\frac{\lambda_1}{1-\lambda_1\Delta t_m }\left(u^{k+1}_{m-1}-\lambda_1\Delta t_m u^{k}_m+ I_{m-1}^m\right) \cr
 I_{m-1}^m =& \lambda_1 \sum_{j=1}^M s_{m,j} u^{k}_j + \lambda_2 \sum_{p=1}^{P} \hat{s}_{m,p} u^{k}_{m,p}
\end{align*}
for the provisional value $u^*_m$ and 
\begin{align*}
 u^{k+1}_{m-1,p}=&u^{k}_{m-1,p-1} +  \frac{\lambda_1\Delta t_{m,p}}{1-\lambda_1\Delta t_m }\left(u^{k+1}_{m-1}-\lambda_1\Delta t_m u^{k}_m+ I_{m-1}^m\right) \cr
   +&\lambda_2 \Delta t_{m,p}\left(u^{k+1}_{m-1,p-1}-u^{k}_{m-1,p-1}\right) + I_{m,p-1}^p \cr
 I_{m,p-1}^p =& \lambda_1 \sum_{j=1}^M \tilde{s}_{m,p,j} u^k_j + \lambda_2 \sum_{q=1}^P s_{m,p,q} u^k_{m,q} %\\
 %\sum_{p=1}^P I_{m,p-1}^p =&\lambda_1 \sum_j^M s_{m,j} u^k_j + \lambda_2 \sum_p^P \hat{s}_{m,p} u^k_{m,p}
\end{align*}
for the embedded values $u_{m,p}$. Now, we collect the embedded values by iteration 
\begin{align*}
 u^{k+1}_{m-1,p}-\frac{\lambda_1 \Delta t_{m,p}}{1-\lambda_1 \Delta t_{m}}u^{k+1}_{m-1}-\lambda_2\Delta t_{m,p} u^{k+1}_{m-1,p-1}=&u^{k}_{m-1,p-1} +  \frac{\lambda_1\Delta t_{m,p}}{1-\lambda_1\Delta t_m }\left(-\lambda_1\Delta t_m u^{k}_m+ I_{m-1}^m\right) \cr
   +&\lambda_2 \Delta t_{m,p}\left(-u^{k}_{m-1,p-1}\right) + I_{m,p-1}^p
\end{align*}
Due to the relation $u^{k+1}_{m}=u^{k+1}_{m-1,P}$, we replace all occurances of $u^{k}_m$ by $u^{k}_{m+1,P}$ and obtain an iterative scheme for the embedded values
\begin{align*}
 u^{k+1}_{m-1,p}-\frac{\lambda_1 \Delta t_{m,p}}{1-\lambda_1 \Delta t_{m}}u^{k+1}_{m-2,P}-\lambda_2\Delta t_{m,p} u^{k+1}_{m-1,p-1}=&u^{k}_{m-1,p-1} -  \frac{\lambda_1\Delta t_{m,p}}{1-\lambda_1\Delta t_m }\lambda_1\Delta t_m u^{k}_{m-1,P} \cr
   +&\frac{\lambda_1\Delta t_{m,p}}{1-\lambda_1\Delta t_m }I_{m-1}^m\cr
   +&\lambda_2 \Delta t_{m,p}\left(-u^{k}_{m-1,p-1}\right) + I_{m,p-1}^p
\end{align*}
\todo{Build vector, write as $\mathbf{u}^{k+1}=G\mathbf{u}^k+d$}
\section*{Two component test problem}
We consider the following simplified test problem to study stability and accuracy
\begin{equation}
	\begin{pmatrix} \dot{y}_1 \\ \dot{y}_2 \end{pmatrix} = \nu \begin{pmatrix} y_1 \\ y_2 \end{pmatrix} + \begin{pmatrix} a \left( y_1 - y_2 \right) \\ -a (y_1 - y_2) \end{pmatrix}
		= \begin{pmatrix}  \nu + a & -a \\ -a & \nu  +a \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}
\end{equation}
assuming that both machine parts have the same heat transport coefficient $\nu$.
The eigenvalues are the roots of
\begin{align*}
	\left| \begin{pmatrix}  \nu + a & -a \\ -a & \nu  +a \end{pmatrix} - \lambda \ve{I} \right| &= \left( \nu + a - \lambda \right)^2 - a^2 \\
	&= \left( \nu + a \right)^2 - 2 \lambda \left( \nu + a \right) + \lambda^2 - a^2 = 0
\end{align*}
so that
\begin{equation}
	\lambda_{1,2} = \nu + a \pm \sqrt{ \left(\nu+a\right)^2 + a^2 - \left( \nu + a \right)^2 } = \nu + a \pm \left| a \right|.
\end{equation}
Assume that $a > 0$ so that
\begin{equation}
	\lambda_1 = \nu \ \text{and} \ \lambda_2 = \nu + 2 a.
\end{equation}
The matrix can be diagonalised as
\begin{equation}
	\begin{pmatrix} \nu+a & -a \\ -a & \nu + a \end{pmatrix} = \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} \nu & 0 \\ 0 & \nu + 2 a \end{pmatrix} \begin{pmatrix} \frac{1}{2} & \frac{1}{2} \\ -\frac{1}{2} & \frac{1}{2} \end{pmatrix} =: S D S^{-1}.
\end{equation}
For some initial value $\ve{y}_0$, the solution $\ve{y} := (y_1, y_2)$ can be found from
\begin{equation}
	\dot{\ve{z}} = D \ve{z}, \ \ve{z}(0) = S^{-1} \ve{y}_0 = \begin{pmatrix} \frac{1}{2} y_1(0) + \frac{1}{2} y_2(0) \\ -\frac{1}{2} y_1(0) + \frac{1}{2} y_2(0) \end{pmatrix}
\end{equation}
by computing $\ve{y} = S \ve{z}$.
It is straightforward to compute
\begin{equation}
	\ve{z} = \begin{pmatrix} z_1(t) \\ z_2(t) \end{pmatrix} = \begin{pmatrix} z_1(0) e^{\nu t} \\ z_2(0) e^{\left(\nu + 2a\right)t} \end{pmatrix} = \begin{pmatrix} \frac{1}{2} \left(y_1(0) + y_2(0) \right) e^{\nu t} \\ \frac{1}{2} \left( y_2(0) - y_1(0) \right) e^{(\nu + 2 a) t} \end{pmatrix}
\end{equation}
and then
\begin{equation}
	\ve{y} = S \ve{z} = \begin{pmatrix} z_{1}(t) - z_{2}(t) \\ z_{1}(t) + z_{2}(t) \end{pmatrix} = \frac{1}{2}  \begin{pmatrix}  y_1(0) + y_2(0)  \\ y_1(0) + y_2(0) \end{pmatrix} e^{\nu t} + \frac{1}{2} \begin{pmatrix} y_1(0) - y_2(0)   \\ y_2(0) - y_1(0) \end{pmatrix} e^{(\nu+2 a)t}
\end{equation}

\paragraph{Stability.}
To assess stability, we need to find a matrix $\ve{R}$ such that
\begin{equation}
	\ve{y}^1 = \ve{R} \ve{y}_0
\end{equation}
for a single time step of length $1$.
If $\sigma(\ve{R}) \leq 1$, the series generated from repeatedly applying $\ve{R}$ won't diverge and the method is stable.

\section{Underlying continuous collocation formula}
The initial value problem over a time step $[T_n, T_{n+1}]$ is
\begin{equation}
	\dot{u}(t) = f(u(t)) + g(u(t), t), \ u(T_n) = u_0.
\end{equation}
The first interpolation polynomial, with respect to the standard nodes, reads
\begin{equation}
	p^s(t) = \sum_{m=1}^{M} u^s_m l_m(t)
\end{equation}
with the Lagrange polynomials $l_m$ defined by $l_m(t_j) = \delta_{mj}$.
Furthermore, we have $M$ interpolation polynomials with respect to the embedded nodes, reading
\begin{equation}
	p^e_m(t) = \sum_{p=1}^{P} u_{m,p} l_{m,p}(t)
\end{equation}
with $l_{m,p}(t_{m,q}) = \delta_{p q}$.
From the ODE in integral form, we have the conditions that
\begin{equation}
	p^s(t_m) = u_0 + \int_{T_n}^{t_m} f(p^e(s))~ds + \int_{T_n}^{t_m} g(p^e(s))~ds
\end{equation}
\todo{is that so? what exactly are the continuous collocation conditions here???}
\end{document}