msdc_heat.tex

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

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

\bibliographystyle{plain}

\begin{document}
 \begin{enumerate}
  \item $\dot u=f(u)+g(u)$, f: schnell, g: langsam
  \item Zeitbereich: $[t_0,t_e]$
  \item Unterteile $[t_0,t_e]$ in N Intervalle: $[t_0,t_1,\dots,t_{N+1}$
  \item Unterteile jedes Teilintervall $[t_m,t_{m+1}]$ in $N_f$ Teilintervalle ($N_{fast}$)
 \end{enumerate}
 Idea stems basicly from \cite{bourlioux2003high}. There are three scales, where both fast scales are threated implicitly.
 Prediktor: (strukturell ist das (ähnlich) eines M[IF]S $\rightarrow$ andere Vorschrift)
 \begin{itemize}
  \item Berechne $u_{n,s}^F$ mit $u_{n,s+1}^F=u_{n,s}^F+\frac{\Delta t}{N_f} (f(u_{n,s}^F)+g(u_n^F))$ (expliziter Euler, mit altem Wert für langsamen Prozess)
  \item Berechne $u_{n+1}$ mit $u_{n+1}=u_{n}+\Delta t(f(u_{n+1})+g(u_{n+1}^F))$ mit $u_{n+1}^F=u_{n,N_{f}}^F$ (TODO: extrapolieren?)
 \end{itemize}

 Korrektor:
 \begin{itemize}
  \item $\delta_{m,s+1}^{F,k}=\delta^{F,k}_{m,s}+\Delta_{s}(f(u_{m,s}^{F,k+1})-f(u_{m,s}^{F,k})+g(u_{m}^{k+1})-g(u_m^k))+E_{m,s+1}(u^k)-E_{m,s}(u^k)$
  \item $\delta_{m+1}^{k+1}=\delta_{m}^{k+1}+\Delta_m (f(u_{m+1}^{F,k+1})-f(u_{m+1}^{F,k})+g(u_{m+1}^{k+1})-g(u_{m+1}^k))+E_{m+1}(u^k)-E_{m}(u^k)$
  \item $\delta_{m,0}^{F,k}=\delta_{m}^k$(???)
 \end{itemize}
 bzw. iterationen:
  \begin{itemize}
   \item $u_{m,0}^{F,k+1}=u_{m}+\Delta_{m,0} (f(u_m))$
   \item $u_{m,s+1}^{F,k+1}=u_{m,s}^{F,k+1}+\Delta_{m,s}(f(u_{m,s}^{F,k+1})-f(u_{m,s}^{F,k})+g(u_{m}^{k+1})-g(u_m^k))+I_{s}^{s+1}(f(u^k)+g(u^k))$
   \item $u_{m+1}^{k+1}=u_{m}^{k+1}+\Delta_{m} (f(\underbrace{u^{F,k+1}_{m+1}}_{u^{F,k+1}_{m,N_f}})-f(\underbrace{u^{F,k}_{m+1}}_{u^{F,k}_{m,N_f}})+g(u_{m+1}^{k+1})-g(u_{m+1}^k))+I_{m}^{m+1}(f(u^k)+g(u^k))$
   \item $I_s^{s+1}(\psi)=\int_{t_{m,s}}^{t_{m,s+1}} \psi(\tau) d\tau\approx \sum_l w_l \psi(c_l)$:
    \begin{itemize}
     \item $u^k$ von grober Skala: $$I_{s}^{s+1}(f(u^k)+g(u^k))\approx \sum_l w_l( f(u^k(c_l))+g(u^k(c_l)))$$, $u^k(t)=\sum_m u^k_m P_m^M(t)$, $P_m^M(t)$: Lagrange polynomial through $c_m$ of degree $M-1$
     \item $u^k$ von feiner Skala: $$I_s^{s+1} (f(u^k)+g(u^k))\approx \sum_l w_l (f(u^{F,k}(c_l))+g(u^{F,k}(c_l)))$$, $u^{F,k}(t)=\sum_m u^{F,k}_m P_m^P(t)$
     \item $u^k$ mixed: A more detailed explaination can be found in \cite{bouzarth2010multirate}, using explicit steps for every part
      \begin{align*}
        f(u^k(\tau))=& \sum_{l=1}^P f(u^{F,k}_{l}) L_l^P(\tau) \cr
         g(u^k(\tau))=&\sum_{l=1}^M g(u^{k}_l) L_l^M(\tau) \cr
 	I_{t_{m,s}}^{t_{m,s+1}}(f(u^k)+g(u^k))\approx & \sum_l^P w_l^P f(u^{F,k}_l) + \sum_{l}^M w_l^{M,s} g(u^k_l) \cr
 	w_l^{M,s}=&\int_{t_{m,s}}^{t_{m,s+1}} L_l^{M}(\tau) d\tau &\text{non classic}\cr
 	w_l^{P}=&\int_{t_{m,s}}^{t_{m,s+1}} L_l^P(\tau) d\tau & \text{classic}
      \end{align*}
    \end{itemize}
  \end{itemize}
  teste an $\dot u=\lambda_fu + \lambda_s u$, $u(0)=1.0$, $t_e=1.0$
 Stabilität, Testproblem: $M\dot u = Lu+B(t)u$, L: symmetrisch \& pos. Def., B(t): symmetrisch $\rightarrow (\lambda_f, \lambda_s)\in \mathbb{R}^2$
 
\paragraph{Some notes.}
Single SDC step for slow process $F_A$, which we treat implicitly even though it is the slow process:
\begin{equation*}
	u_{A,m+1}^{k+1} = u_{A,m}^{k+1} + \Delta t_m \left( F_{A}(u^{k+1}_{A,m+1}) -F_{A}(u^{k}_{m+1}) \right) + I_m^{m+1}\left( F_{A}(u^k) + F_{D}(u^k) \right).
\end{equation*}
Note that $u^k_{m+1}$ is equal to the \emph{fine} level solution at $t_{m+1}$ and is constructed from one of the sub collocation rules.
Single SDC step for fast process $F_D$, which we treat explicitly with small steps:
\begin{align*}
	u_{p+1}^{k+1} = u_{p}^{k+1} &+ \Delta t_p \left( F_{D}(u^{k+1}_{p}) - F_{D}(u^k_p) \right) \\
							   &+ \Delta t_p \left( F_{A}(u^{k+1}_{A,m+1}) - F_{A}(u^k_m) \right) \\
							   &+ I_p^{p+1}\left( F_A(u^k) + F_D(u^k) \right)
\end{align*}
The solution on the sub collocation nodes is the ``real" solution, that is $u^k$ contains all fine level solutions at $k$.
As $u^{k+1} - u^k \to 0$ and also $u^{k+1}_A - u^k \to 0$, we have
\begin{align*}
	u_{m+1} &= u_{m} + I_m^{m+1} \left( F_A(u) + F_D(u) \right) \\
	u_{p+1} &= u_{p} + I_{p}^{p+1} \left( F_A(u) + F_D(u) \right)
\end{align*}
so that we get convergence to the coarse and fine level collocation solutions.

\subsection{Multi-rate SDC.}
\todo{The indices do not really match -- should maybe consider $t_{m-1}$ to $t_m$ and $t_{m,p-1}$ to $t_{m,p}$ instead}
The Lagrange polynomials on $[T_n, T_{n+1}]$
\begin{equation}
	l_m(t), m = 1, \ldots, M \quad \text{with} \quad l_m(t_j) = \delta_{mj}.
\end{equation}
For each sub-step $[t_m, t_{m+1}]$, there are $P$ many Lagrange polynomials for the nodes $t_{m} \leq t_{m,1} < \ldots < t_{m,P} \leq t_{m+1}$:
\begin{equation}
	l_{m,p}(t), p=1, \ldots, P \quad \text{with} \quad l_{m,p}(t_{m,j}) = \delta_{pj}.
\end{equation}
Note that there are $M$ many polynomials $l_m(s)$ and $M \times P$ many $l_{m,p}$.
Also, there are $M$ many nodes $t_m$ and $M \times P$ many nodes $t_{m,p}$.
The ``standard" quadrature weights are
\begin{equation}
	s_{m,n} = \int_{t_m}^{t_{m+1}} l_n(s)~ds, \ m=1, \ldots, M; \ n=1, \ldots, M
\end{equation}
and
\begin{equation}
	s_{m,p,j} = \int_{t_{m,p}}^{t_{m,p+1}}  l_{m,j}(s)~ds, \ m=1, \ldots, M; \ p=1, \ldots, P; \ j=1, \ldots, P.
\end{equation}
In addition, we will need the ``mixed" weights
\begin{equation}
	\tilde{s}_{m,n,p} := \int_{t_{m,p}}^{t_{m,p+1}} l_{n}(s)~ds, m=1, \ldots, M; \ n=1, \ldots, M; p=1, \ldots, P
\end{equation}
This leads to the following split approximations of integrals
\begin{equation}
	\int_{t_{m,p}}^{t_{m,p+1}} u_1(s) + u_2(s)~ds \approx \sum_{n=1}^{M} \tilde{s}_{m,n,p} u_1(t_{n}) + \sum_{j=1}^{P} s_{m,p,j} u_2(t_{m,j}) =: I_{m,p}^{p+1}(u_1,u_{m,2})
\end{equation}
and
\begin{align*}
	\int_{t_m}^{t_{m+1}} u_1(s) + u_2(s)~ds &\approx \sum_{n=1}^{M} s_{m,n} u_1(t_n) + \sum_{p=1}^{P} \sum_{j=1}^{P} s_{m,p,j} u_2(t_{m,j}) \\
					&=  \sum_{n=1}^{M} s_{m,n} u_1(t_n) + \sum_{j=1}^{P} u_2(t_{m,j}) \left( \sum_{p=1}^{P} s_{m,p,j} \right) \\
					&=: \sum_{n=1}^{M} s_{m,n} u_1(t_n) + \sum_{j=1}^{P} \hat{s}_{m,j} u_2(t_{m,j}) =: I_m^{m+1}(u_1, u_{m,2})
\end{align*}
with $u_1 \in \mathbb{R}^M$ and $u_2 \in \mathbb{R}^{MP}$.
\begin{remark}
The class \texttt{MultirateCollocation.py} contains the weights $s_{m,n}$, $s_{m,p,j}$ and $\tilde{s}_{m,n,p}$ and functions to evaluate the corresponding quadrature rules.
\end{remark}
Now consider an IVP 
\begin{equation}
	\dot{u}(t) = f_1(u(t)) + f_2(u(t))
\end{equation}
where we want to integrate $f_1$ implicitly with large and $f_2$ explicitly with small time steps.
Denote as $\ve{u}_m \in \mathbb{R}^P$ the vector containing solutions at all sub-steps between two nodes of the coarse level collocation rule, that is $t_m \leq t_{m,1} < \ldots < t_{m,P} \leq t_{m+1}$ for $m=1, \ldots, M$.
Furthermore, denote $\ve{u} \in \mathbb{R}^M$ the vector containing solutions at the quadrature nodes $T_{n} \leq t_{1} < \ldots < t_{M} \leq T_{n+1}$.
Given $\ve{u}^k_m$ and $\ve{u}^k$, we first compute the entries of $\ve{u}^{k+1}$ from
\begin{equation}
	\label{eq:coarse_step}
	u^{k+1}_{m+1} = u^{k+1}_m + \Delta t_m \left( f_1(u^{k+1}_{m+1}) - f_1(u^k_{m+1}) \right) + I_{m}^{m+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right).
\end{equation}
Then, we can perform the SDC sweep on the sub-step
\begin{align}
	\label{eq:mr_sdc_sweep}
	u_{m,p+1}^{k+1} = u_{m,p}^{k+1} &+ \Delta t_{m,p} \left( f_1(u^{k+1}_{m+1}) - f_1(u^k_{m+1}) \right) \\
					&+ \Delta t_{m,p} \left( f_2(u^{k+1}_{p}) - f_2(u^k_{p}) \right) \\
					&+ I_{m,p}^{p+1}\left( f_1(\ve{u}^k) + f_2(\ve{u}_m^k) \right)
\end{align}
for $p=1, \ldots, P$.
\begin{algorithm2e}[h]
	\caption{Multirate SDC sweeps}
	Run predictor to generate $\ve{u}^0_m$ and $\ve{u}^0$\\
	\For{$k=1, K$}{
		Update $I_{m,p}^{p+1}(f_1(\ve{u}^k) + f_2(\ve{u}^k_m))$\\
		\For{$m=1, M$}{
			Compute $u^{k+1}_{m+1}$ using~\eqref{eq:coarse_step}\\
			Compute $f_1(u^{k+1}_{m+1})$\\
			\For{$p=1, P$}{
				Update $u^{k+1}_{m,p+1}$ using~\eqref{eq:mr_sdc_sweep}
			}
		}
	}
\end{algorithm2e}
Upon convergence, as $u^{k+1}_{m} - u^k_m \to 0$ and $u^{k+1}_{m,p} - u^k_{m,p} \to 0$, the solutions approach
\begin{equation}
	u_{m+1} = u_{m} + I_m^{m+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right)
\end{equation}
and
\begin{equation}
	u_{m,p+1} = u_{m,p} + I_{m,p}^{p+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right).
\end{equation}
\begin{remark}
Is there any reason to expect that if $t_{m,P} = t_{m+1}$ we will get $u_{m,P} = u_{m+1}$? Consider that
\begin{equation}
	u_{m,P} = u_{m} + \sum_{p=1}^{P} I_{m,p}^{p+1}
\end{equation}
(using the sub-step update recursively) and
\begin{equation}
	u_{m+1} = u_m + I_{m}^{m+1}
\end{equation}
so that
\begin{equation}
	u_{m+1} - u_{m,P} = I_{m}^{m+1} - \sum_{p=1}^{P} I_{m,p}^{p+1}
\end{equation}
The question becomes whether $\sum_{p=1}^{P} I_{m,p}^{p+1} = I_{m}^{m+1}$.
Obviously, this is true in the continuous case.
Also, because of the definition of the $\hat{s}_{m,j}$, it holds true for the $u_2$ part in both operators.
For the $u_1$ part, if $t_{m,1} = t_m$ and $t_{m,P+1} = t_{m+1}$, we have
\begin{equation}
	s_{m,n} = \int_{t_m}^{t_{m+1}} l_n(s)~ds = \sum_{p=1}^{P} \int_{t_{m,p}}^{t_{m,p+1}} l_n(s)~ds = \sum_{p=1}^{P} \tilde{s}_{m,n,p}
\end{equation}
so that the condition is also satisfied for the $u_1$ term.
\end{remark}

\subsubsection{A note on testing.}

\paragraph{Coarse and fine sweep match.}
Setting $f_2 \equiv 0$ reduces the coarse level sweep to
\begin{equation}
	u^{k+1}_{m+1} = u^{k}_m + \Delta t_m \left( f_1(u^{k+1}_{m+1}) - f_1(u^k_{m+1}) \right) + I_m^{m+1}(f_1(\ve{u}^k))
\end{equation}
and thus a standard implicit Euler SDC.
The coarse level values should therefore show the expected convergence behaviour in this case.
The sub-steps become
\begin{align}
	u^{k+1}_{m,p+1} &= u^{k+1}_{m,p} + \Delta t_{m,p} \left( f_1(u^{k+1}_{m+1}) - f_1(u^{k}_{m+1}) \right) + I_{m, p}^{p+1} (f_1(\ve{u}^k)) \\
					&= u^{k+1}_m + \left( f_1(u^{k+1}_{m+1}) - f_1(u^{k}_{m+1}) \right) \sum_{q=1}^{p} \Delta t_{m,p} + \sum_{q=1}^{p} I_{m, p}^{p+1} (f_1(\ve{u}^k))
\end{align}
so that the final sub-step value should be
\begin{equation}
	u^{k+1}_{m,P} = u^{k+1}_{m} +  \left( f_1(u^{k+1}_{m+1}) - f_1(u^{k}_{m+1}) \right) \Delta t_{m} + I_m^{m+1} (f_1(\ve{u}^k))
\end{equation}
using that $\sum_{p=1}^{P} \Delta t_{m,p} = \Delta t_m$ (for Radau nodes at least) and that
\begin{equation}
	\sum_{p=1}^{P} I_{m,p}^{p+1} = I_{m}^{m+1}.
\end{equation}
Therefore, if $f_2 \equiv 0$, the sub-stepping should provide the same values at the coarse quadrature nodes as the coarse sweep.

\paragraph{Coarse sweep invariance.}
If the collocation solution is used for $\ve{u}^k$ in the coarse level sweep, it remains unchanged.
The collocation solution satisfies
\begin{equation}
	u_m = u_{0} + \sum_{j=1}^{M} q_{m,j} f(u_j) = u_{m-1} + \sum_{j=1}^{M} s_{m,j} f(u_j).		
\end{equation}
Inserting it into the implicit Euler SDC sweep gives, assuming $u^{k+1}_{m-1} = u_{m-1}$ (induction)
\begin{align}
	u^{k+1}_{m} &= u^{k+1}_{m-1} + \Delta t_m \left( f_1(u^{k+1}_m) - f_1(u_m) \right) + \sum_{j=1}^{M} s_{m,j} f_1(u_j) \\
	u^{k+1}_m - \Delta t_m f_1(u^{k+1}_{m}) &= u_m - \Delta t_m f_1(u_m)
\end{align}
from which it follows that $u^{k+1}_m = u_m$.

\paragraph{Fine sweep invariance.}
Assume that in each sub-step the solution is the collocation solution given by
\begin{equation}
	u_{m,p} = u_{m,p-1} + \sum_{j=1}^{P} s_{m,p,j} f_2(u_{m,j})
\end{equation}
with some given initial value $u_{m,0}$.
If we assume that $f_1 \equiv 0$, the sub-step sweep becomes
\begin{align*}
	u^{k+1}_{m,p} &= u^{k+1}_{m,p-1} + \Delta t_{m,p} \left( f_2(u^{k+1}_{m,p-1}) - f_2(u_{m,p-1}) \right) + \sum_{j=1}^{P} s_{m,p,j} f_2(u_{m,j}) \\
				  &= u_{m,p}  + \Delta t_{m,p} \left( f_2(u^{k+1}_{m,p-1}) - f_2(u_{m,p-1}) \right) \\
				  &= u_{m,p}
\end{align*}
assuming that $u^{k+1}_{m,p-1} = u_{m,p-1}$ (induction).

\paragraph{Coarse sweep residual.}

\paragraph{Fine sweep residual.}
Assume the solution at sub-steps is given by
\begin{equation}
	u_{m,p} = u_{m,p-1} + \sum_{j=1}^{P} s_{m,p,j} f_2(u_{m,j}) = \ldots = u_{m,0} +  \sum_{j=1}^{P} q_{m,p,j} f_2(u_{m,j})
\end{equation}
with given initial value $u_{m,0}$.
Then, the residual with entries
\begin{equation}
	r_{m,p} := u_{m,p} - u_{m,p-1} - \sum_{j=1}^{P} s_{m,p,j} f_2(u_{m,j})
\end{equation}
has to be zero.

\paragraph{Multi-rate collocation solution.}
Assume that we have solutions $\ve{u} \in \mathbb{R}^{M}$ and $\ve{u}_m \in \mathbb{R}^{P}$, $m=1, \ldots, M$ that satisfy
\begin{align*}
	u_{m,p} &= u_{m,p-1} + I_{m,p}^{p+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right) \\
			&= u_{m,0} + \sum_{q=1}^{p} I_{m,q}^{q+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right).
\end{align*}
Further, assume that the coarse quadrature entries satisfy
\begin{align*}
	u_{m} &= u_{m-1} + I_{m}^{m+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right) \\
		  &= u_0 + \sum_{j=1}^{m} I_{j}^{j+1} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right)
\end{align*}
As above, the coarse level sweep results in
\begin{align*}
	u_{m}^{k+1} &= u^{k+1}_{m-1} + \Delta t_m \left( f_1(u^{k+1}_m) - f_1(u_m) \right) + I_m^{m+1}\left( f_1(\ve{u}) + f_2(\ve{u}_m) \right) \\
				&= u_{m} + \Delta t_m \left( f_1(u^{k+1}_m) - f_1(u_m) \right) \\
\Rightarrow u^{k+1}_m - \Delta t_m f_1(u^{k+1}_{m}) &= u_m - \Delta t_m f_1(u_m) \\
\Rightarrow u^{k+1}_m &= u_m
\end{align*}
assuming $u^{k+1}_{m-1} = u_{m-1}$ (induction).
Therefore, for the sub-step sweep, the $f_1$ implicit Euler terms cancel out and we get
\begin{align*}
	u_{m,p}^{k+1} &= u_{m,p-1}^{k+1} + \Delta t_m \left( f_2(u_{m,p-1}^{k+1}) - f_2(u_{m,p-1}) \right) + I_{m,p-1}^{p} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right).
\end{align*}
Assuming $u^{k+1}_{m,p-1} = u_{m,p-1}$ (induction), we get
\begin{equation}
	u_{m,p}^{k+1} = u_{m,p-1} + I_{m,p-1}^{p} \left( f_1(\ve{u}) + f_2(\ve{u}_m) \right) = u_{m,p}
\end{equation}
\bibliography{refs}
\end{document}